kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

0.1.4-dev-4 #86

Merged
altavir merged 38 commits from dev into master 2020-04-30 11:36:59 +03:00
Conversation 0 Commits 38 Files Changed 72 +2730 -152
72 changed files with 2730 additions and 152 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
README.md
View File

@ -1,9 +1,12 @@
[![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)
[![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
# 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.

4
build.gradle.kts
View File

@ -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")

18
examples/build.gradle.kts
View File

@ -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<AllOpenExtension> {
@ -13,10 +13,9 @@ configure<AllOpenExtension> {
}
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)
}

71
examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt Normal file
View File

@ -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()
}
}
}

8
examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt Normal file
View File

@ -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")
}

30
examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt
View File

@ -4,7 +4,13 @@ import kotlinx.coroutines.GlobalScope
import scientifik.kmath.operations.RealField
import kotlin.system.measureTimeMillis
fun main(args: Array<String>) {
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<String>) {
//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<String>) {
}
}
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<Double> = one
repeat(n) {
@ -45,10 +46,7 @@ fun main(args: Array<String>) {
}
}
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<String>) {
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<Double> = one
@ -72,6 +68,4 @@ fun main(args: Array<String>) {
}
}
println("Generic addition completed in $genericTime millis")
}

35
examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt Normal file
View File

@ -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<Double, RealField>.simple() {
val m1 = produce<D2, D3> { i, j -> (i + j).toDouble() }
val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
//Dimension-safe addition
m1.transpose() + m2
}
object D5 : Dimension {
override val dim: UInt = 5u
}
fun DMatrixContext<Double, RealField>.custom() {
val m1 = produce<D2, D5> { i, j -> (i + j).toDouble() }
val m2 = produce<D5, D2> { i, j -> (i - j).toDouble() }
val m3 = produce<D2, D2> { i, j -> (i - j).toDouble() }
(m1 dot m2) + m3
}
fun main() {
DMatrixContext.real.run {
simple()
custom()
}
}

BIN
gradle/wrapper/gradle-wrapper.jar vendored
View File

Binary file not shown.

2
gradle/wrapper/gradle-wrapper.properties vendored
View File

@ -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

3
gradlew.bat vendored
View File

@ -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"

3
kmath-commons/build.gradle.kts
View File

@ -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")
}

2
kmath-commons/src/test/kotlin/scientifik/kmath/commons/expressions/AutoDiffTest.kt
View File

@ -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 <R> diff(order: Int, vararg parameters: Pair<String, Double>, block: DerivativeStructureField.() -> R) =

10
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt
View File

@ -41,16 +41,6 @@ fun <T : Any> Structure2D.Companion.square(vararg elements: T): FeaturedMatrix<T
return BufferMatrix(size, size, buffer)
}
fun <T : Any> Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)
class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
operator fun invoke(vararg elements: T): FeaturedMatrix<T> {
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()
/**

0
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt → kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt
View File

14
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt Normal file
View File

@ -0,0 +1,14 @@
package scientifik.kmath.linear
import scientifik.kmath.structures.Structure2D
import scientifik.kmath.structures.asBuffer
class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
operator fun invoke(vararg elements: T): FeaturedMatrix<T> {
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 <T : Any> Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)

37
kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/functions.kt
View File

@ -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<T, C : SpaceOperations<T>> {
suspend operator fun C.invoke(arg: T): T
}
suspend fun UFunction<Double, RealField>.invoke(arg: Double) = RealField.invoke(arg)
/**
* A suspendable multivariate (N->1) function defined on algebraic context
*/
interface MFunction<T, C : SpaceOperations<T>> {
/**
* The input dimension of the function
*/
val dimension: UInt
suspend operator fun C.invoke(vararg args: T): T
}
suspend fun MFunction<Double, RealField>.invoke(args: DoubleArray) = RealField.invoke(*args.toTypedArray())
@JvmName("varargInvoke")
suspend fun MFunction<Double, RealField>.invoke(vararg args: Double) = RealField.invoke(*args.toTypedArray())
/**
* A suspendable univariate function with parameter
*/
interface ParametricUFunction<T, P, C : SpaceOperations<T>> {
suspend operator fun C.invoke(arg: T, parameter: P): T
}

6
kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
View File

@ -6,14 +6,14 @@ annotation class KMathContext
/**
* Marker interface for any algebra
*/
interface Algebra
interface Algebra<T>
inline operator fun <T : Algebra, R> T.invoke(block: T.() -> R): R = run(block)
inline operator fun <T : Algebra<*>, R> T.invoke(block: T.() -> R): R = run(block)
/**
* Space-like operations without neutral element
*/
interface SpaceOperations<T> : Algebra {
interface SpaceOperations<T> : Algebra<T> {
/**
* Addition operation for two context elements
*/

11
kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt
View File

@ -2,3 +2,14 @@ package scientifik.kmath.operations
fun <T> Space<T>.sum(data: Iterable<T>): T = data.fold(zero) { left, right -> add(left, right) }
fun <T> Space<T>.sum(data: Sequence<T>): T = data.fold(zero) { left, right -> add(left, right) }
fun <T : Any, S : Space<T>> Iterable<T>.sumWith(space: S): T = space.sum(this)
//TODO optimized power operation
fun <T> RingOperations<T>.power(arg: T, power: Int): T {
var res = arg
repeat(power - 1) {
res *= arg
}
return res
}

484
kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt Normal file
View File

@ -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<BigInt> {
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<BigInt> {
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<BigInt, BigInt> {
// 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<UInt, String> = 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
}

4
kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt
View File

@ -29,7 +29,7 @@ fun <T : MathElement<out TrigonometricOperations<T>>> ctg(arg: T): T = arg.conte
/**
* A context extension to include power operations like square roots, etc
*/
interface PowerOperations<T> {
interface PowerOperations<T> : Algebra<T> {
fun power(arg: T, pow: Number): T
fun sqrt(arg: T) = power(arg, 0.5)
@ -42,7 +42,7 @@ fun <T : MathElement<out PowerOperations<T>>> sqr(arg: T): T = arg pow 2.0
/* Exponential */
interface ExponentialOperations<T> {
interface ExponentialOperations<T>: Algebra<T> {
fun exp(arg: T): T
fun ln(arg: T): T
}

2
kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
View File

@ -73,6 +73,8 @@ fun <T> Buffer<T>.asSequence(): Sequence<T> = Sequence(::iterator)
fun <T> Buffer<T>.asIterable(): Iterable<T> = asSequence().asIterable()
val Buffer<*>.indices: IntRange get() = IntRange(0, size - 1)
interface MutableBuffer<T> : Buffer<T> {
operator fun set(index: Int, value: T)

1
kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
View File

@ -2,7 +2,6 @@ package scientifik.kmath.structures
import scientifik.kmath.operations.*
interface ExtendedNDField<T : Any, F, N : NDStructure<T>> :
NDField<T, F, N>,
TrigonometricOperations<N>,

48
kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
View File

@ -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<Double>(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])
}
}
}

50
kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntAlgebraTest.kt Normal file
View File

@ -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")
}
}
}

26
kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConstructorTest.kt Normal file
View File

@ -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)
}
}

43
kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConversionsTest.kt Normal file
View File

@ -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())
}
}

381
kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntOperationsTest.kt Normal file
View File

@ -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)
}
}

0
kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumers.kt → kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt
View File

28
kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt
View File

@ -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<out R> {
interface Chain<out R>: Flow<R> {
/**
* Generate next value, changing state if needed
*/
@ -36,14 +38,15 @@ interface Chain<out R> {
*/
fun fork(): Chain<R>
@InternalCoroutinesApi
override suspend fun collect(collector: FlowCollector<R>) {
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 <R> Chain<R>.flow(): Flow<R> = kotlinx.coroutines.flow.flow { while (true) emit(next()) }
fun <T> Iterator<T>.asChain(): Chain<T> = SimpleChain { next() }
fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain()
@ -127,6 +130,21 @@ fun <T, R> Chain<T>.map(func: suspend (T) -> R): Chain<R> = object : Chain<R> {
override fun fork(): Chain<R> = 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 <T> Chain<T>.filter(block: (T) -> Boolean): Chain<T> = object : Chain<T> {
override suspend fun next(): T {
var next: T
do {
next = this@filter.next()
} while (!block(next))
return next
}
override fun fork(): Chain<T> = this@filter.fork().filter(block)
}
/**
* Map the whole chain
*/

9
kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt
View File

@ -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) {

2
kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/RingBufferTest.kt
View File

@ -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 {

19
kmath-dimensions/build.gradle.kts Normal file
View File

@ -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"))
}
}
}

35
kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt Normal file
View File

@ -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 <D : Dimension> KClass<D>.dim(): UInt = Dimension.resolve(this).dim
expect fun <D : Dimension> Dimension.Companion.resolve(type: KClass<D>): D
expect fun Dimension.Companion.of(dim: UInt): Dimension
inline fun <reified D : Dimension> 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
}

161
kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt Normal file
View File

@ -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<T, R : Dimension, C : Dimension> : Structure2D<T> {
companion object {
/**
* Coerces a regular matrix to a matrix with type-safe dimensions and throws a error if coercion failed
*/
inline fun <T, reified R : Dimension, reified C : Dimension> coerce(structure: Structure2D<T>): DMatrix<T, R, C> {
if (structure.rowNum != Dimension.dim<R>().toInt()) {
error("Row number mismatch: expected ${Dimension.dim<R>()} but found ${structure.rowNum}")
}
if (structure.colNum != Dimension.dim<C>().toInt()) {
error("Column number mismatch: expected ${Dimension.dim<C>()} but found ${structure.colNum}")
}
return DMatrixWrapper(structure)
}
/**
* The same as [coerce] but without dimension checks. Use with caution
*/
fun <T, R : Dimension, C : Dimension> coerceUnsafe(structure: Structure2D<T>): DMatrix<T, R, C> {
return DMatrixWrapper(structure)
}
}
}
/**
* An inline wrapper for a Matrix
*/
inline class DMatrixWrapper<T, R : Dimension, C : Dimension>(
val structure: Structure2D<T>
) : DMatrix<T, R, C> {
override val shape: IntArray get() = structure.shape
override fun get(i: Int, j: Int): T = structure[i, j]
}
/**
* Dimension-safe point
*/
interface DPoint<T, D : Dimension> : Point<T> {
companion object {
inline fun <T, reified D : Dimension> coerce(point: Point<T>): DPoint<T, D> {
if (point.size != Dimension.dim<D>().toInt()) {
error("Vector dimension mismatch: expected ${Dimension.dim<D>()}, but found ${point.size}")
}
return DPointWrapper(point)
}
fun <T, D : Dimension> coerceUnsafe(point: Point<T>): DPoint<T, D> {
return DPointWrapper(point)
}
}
}
/**
* Dimension-safe point wrapper
*/
inline class DPointWrapper<T, D : Dimension>(val point: Point<T>) :
DPoint<T, D> {
override val size: Int get() = point.size
override fun get(index: Int): T = point[index]
override fun iterator(): Iterator<T> = point.iterator()
}
/**
* Basic operations on dimension-safe matrices. Operates on [Matrix]
*/
inline class DMatrixContext<T : Any, Ri : Ring<T>>(val context: GenericMatrixContext<T, Ri>) {
inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
if (rowNum != Dimension.dim<R>().toInt()) {
error("Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum")
}
if (colNum != Dimension.dim<C>().toInt()) {
error("Column number mismatch: expected ${Dimension.dim<C>()} but found $colNum")
}
return DMatrix.coerceUnsafe(this)
}
/**
* Produce a matrix with this context and given dimensions
*/
inline fun <reified R : Dimension, reified C : Dimension> produce(noinline initializer: (i: Int, j: Int) -> T): DMatrix<T, R, C> {
val rows = Dimension.dim<R>()
val cols = Dimension.dim<C>()
return context.produce(rows.toInt(), cols.toInt(), initializer).coerce<R,C>()
}
inline fun <reified D : Dimension> point(noinline initializer: (Int) -> T): DPoint<T, D> {
val size = Dimension.dim<D>()
return DPoint.coerceUnsafe(
context.point(
size.toInt(),
initializer
)
)
}
inline infix fun <reified R1 : Dimension, reified C1 : Dimension, reified C2 : Dimension> DMatrix<T, R1, C1>.dot(
other: DMatrix<T, C1, C2>
): DMatrix<T, R1, C2> {
return context.run { this@dot dot other }.coerce()
}
inline infix fun <reified R : Dimension, reified C : Dimension> DMatrix<T, R, C>.dot(vector: DPoint<T, C>): DPoint<T, R> {
return DPoint.coerceUnsafe(context.run { this@dot dot vector })
}
inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, R, C>.times(value: T): DMatrix<T, R, C> {
return context.run { this@times.times(value) }.coerce()
}
inline operator fun <reified R : Dimension, reified C : Dimension> T.times(m: DMatrix<T, R, C>): DMatrix<T, R, C> =
m * this
inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.plus(other: DMatrix<T, C, R>): DMatrix<T, C, R> {
return context.run { this@plus + other }.coerce()
}
inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.minus(other: DMatrix<T, C, R>): DMatrix<T, C, R> {
return context.run { this@minus + other }.coerce()
}
inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.unaryMinus(): DMatrix<T, C, R> {
return context.run { this@unaryMinus.unaryMinus() }.coerce()
}
inline fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.transpose(): DMatrix<T, R, C> {
return context.run { (this@transpose as Matrix<T>).transpose() }.coerce()
}
/**
* A square unit matrix
*/
inline fun <reified D : Dimension> one(): DMatrix<T, D, D> = produce { i, j ->
if (i == j) context.elementContext.one else context.elementContext.zero
}
inline fun <reified R : Dimension, reified C : Dimension> zero(): DMatrix<T, R, C> = produce { _, _ ->
context.elementContext.zero
}
companion object {
val real = DMatrixContext(MatrixContext.real)
}
}

30
kmath-dimensions/src/commonTest/kotlin/scientifik/dimensions/DMatrixContextTest.kt Normal file
View File

@ -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<D2, D2> { 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<D2, D3> { i, j -> (i + j).toDouble() }
val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
//Dimension-safe addition
m1.transpose() + m2
}
}
}

22
kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt Normal file
View File

@ -0,0 +1,22 @@
package scientifik.kmath.dimensions
import kotlin.reflect.KClass
private val dimensionMap = hashMapOf<UInt, Dimension>(
1u to D1,
2u to D2,
3u to D3
)
@Suppress("UNCHECKED_CAST")
actual fun <D : Dimension> Dimension.Companion.resolve(type: KClass<D>): 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
}
}
}

18
kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt Normal file
View File

@ -0,0 +1,18 @@
package scientifik.kmath.dimensions
import kotlin.reflect.KClass
actual fun <D:Dimension> Dimension.Companion.resolve(type: KClass<D>): 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
}
}
}

11
kmath-for-real/build.gradle.kts Normal file
View File

@ -0,0 +1,11 @@
plugins {
id("scientifik.mpp")
}
kotlin.sourceSets {
commonMain {
dependencies {
api(project(":kmath-core"))
}
}
}

0
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt → kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt
View File

146
kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt Normal file
View File

@ -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<DoubleArray>.toMatrix() = toList().let {
MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] }
}
fun Matrix<Double>.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<Double>.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row, col] * double
}
operator fun Matrix<Double>.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row, col] + double
}
operator fun Matrix<Double>.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row, col] - double
}
operator fun Matrix<Double>.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row, col] / double
}
operator fun Double.times(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
row, col -> this * matrix[row, col]
}
operator fun Double.plus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
row, col -> this * matrix[row, col]
}
operator fun Double.minus(matrix: Matrix<Double>) = 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<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
// row, col -> matrix[row, col] / this
//}
/*
* Per-element (!) square and power operations
*/
fun Matrix<Double>.square() = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row, col].pow(2)
}
fun Matrix<Double>.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) {
i, j -> this[i, j].pow(n)
}
/*
* Operations on two matrices (per-element!)
*/
operator fun Matrix<Double>.times(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row, col] * other[row, col]
}
operator fun Matrix<Double>.plus(other: Matrix<Double>) = MatrixContext.real.add(this, other)
operator fun Matrix<Double>.minus(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) {
row, col -> this[row,col] - other[row,col]
}
/*
* Operations on columns
*/
inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double>) -> Double) =
MatrixContext.real.produce(rowNum,colNum+1) {
row, col ->
if (col < colNum)
this[row, col]
else
mapper(rows[row])
}
fun Matrix<Double>.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) {
row, col -> this[row, columnRange.first + col]
}
fun Matrix<Double>.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex)
fun Matrix<Double>.sumByColumn() = MatrixContext.real.produce(1, colNum) { _, j ->
val column = columns[j]
with(elementContext) {
sum(column.asSequence())
}
}
fun Matrix<Double>.minByColumn() = MatrixContext.real.produce(1, colNum) {
_, j -> columns[j].asSequence().min() ?: throw Exception("Cannot produce min on empty column")
}
fun Matrix<Double>.maxByColumn() = MatrixContext.real.produce(1, colNum) {
_, j -> columns[j].asSequence().max() ?: throw Exception("Cannot produce min on empty column")
}
fun Matrix<Double>.averageByColumn() = MatrixContext.real.produce(1, colNum) {
_, j -> columns[j].asSequence().average()
}
/*
* Operations processing all elements
*/
fun Matrix<Double>.sum() = elements().map { (_, value) -> value }.sum()
fun Matrix<Double>.min() = elements().map { (_, value) -> value }.min()
fun Matrix<Double>.max() = elements().map { (_, value) -> value }.max()
fun Matrix<Double>.average() = elements().map { (_, value) -> value }.average()

16
kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt Normal file
View File

@ -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)
}
}

36
kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt Normal file
View File

@ -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])
}
}

0
kmath-for-real/src/jvmMain/kotlin/.gitkeep Normal file
View File

11
kmath-functions/build.gradle.kts Normal file
View File

@ -0,0 +1,11 @@
plugins {
id("scientifik.mpp")
}
kotlin.sourceSets {
commonMain {
dependencies {
api(project(":kmath-core"))
}
}
}

57
kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt Normal file
View File

@ -0,0 +1,57 @@
package scientifik.kmath.functions
import scientifik.kmath.operations.Ring
interface Piecewise<T, R> {
fun findPiece(arg: T): R?
}
interface PiecewisePolynomial<T : Any> :
Piecewise<T, Polynomial<T>>
/**
* Ordered list of pieces in piecewise function
*/
class OrderedPiecewisePolynomial<T : Comparable<T>>(delimeter: T) :
PiecewisePolynomial<T> {
private val delimiters: ArrayList<T> = arrayListOf(delimeter)
private val pieces: ArrayList<Polynomial<T>> = 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<T>) {
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<T>) {
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<T>? {
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 <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
findPiece(arg)?.value(ring, arg)
fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }

73
kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt Normal file
View File

@ -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<T : Any>(val coefficients: List<T>) {
constructor(vararg coefficients: T) : this(coefficients.toList())
}
fun Polynomial<Double>.value() =
coefficients.reduceIndexed { index: Int, acc: Double, d: Double -> acc + d.pow(index) }
fun <T : Any, C : Ring<T>> Polynomial<T>.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 <T : Any, C : Ring<T>> Polynomial<T>.asMathFunction(): MathFunction<T, out C, T> = object :
MathFunction<T, C, T> {
override fun C.invoke(arg: T): T = value(this, arg)
}
/**
* Represent the polynomial as a regular context-less function
*/
fun <T : Any, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
/**
* An algebra for polynomials
*/
class PolynomialSpace<T : Any, C : Ring<T>>(val ring: C) : Space<Polynomial<T>> {
override fun add(a: Polynomial<T>, b: Polynomial<T>): Polynomial<T> {
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<T>, k: Number): Polynomial<T> {
ring.run {
return Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * k })
}
}
override val zero: Polynomial<T> =
Polynomial(emptyList())
operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
}
fun <T : Any, C : Ring<T>, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R {
return PolynomialSpace(this).run(block)
}

33
kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt Normal file
View File

@ -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<T, C : Algebra<T>, R> {
operator fun C.invoke(arg: T): R
}
fun <R> MathFunction<Double, RealField, R>.invoke(arg: Double): R = RealField.invoke(arg)
/**
* A suspendable function defined in algebraic context
*/
interface SuspendableMathFunction<T, C : Algebra<T>, R> {
suspend operator fun C.invoke(arg: T): R
}
suspend fun <R> SuspendableMathFunction<Double, RealField, R>.invoke(arg: Double) = RealField.invoke(arg)
/**
* A parametric function with parameter
*/
interface ParametricFunction<T, P, C : Algebra<T>> {
operator fun C.invoke(arg: T, parameter: P): T
}

45
kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt Normal file
View File

@ -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<X, Y> {
fun interpolate(points: XYPointSet<X, Y>): (X) -> Y
}
interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T> {
val algebra: Ring<T>
fun getDefaultValue(): T = error("Out of bounds")
fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T>
override fun interpolate(points: XYPointSet<T, T>): (T) -> T = { x ->
interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue()
}
}
fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
x: Buffer<T>,
y: Buffer<T>
): PiecewisePolynomial<T> {
val pointSet = BufferXYPointSet(x, y)
return interpolatePolynomials(pointSet)
}
fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
data: Map<T, T>
): PiecewisePolynomial<T> {
val pointSet = BufferXYPointSet(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
return interpolatePolynomials(pointSet)
}
fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
data: List<Pair<T, T>>
): PiecewisePolynomial<T> {
val pointSet = BufferXYPointSet(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
return interpolatePolynomials(pointSet)
}

26
kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt Normal file
View File

@ -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<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
override fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T> = 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)
}
}
}
}

296
kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt Normal file
View File

@ -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<T : Comparable<T>>(override val algebra: Ring<T>) : PolynomialInterpolator<T> {
// /**
// * 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<Double, Double, Double>.smooth(): XYPointSet<Double, Double> {
// 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<Double>, weights: Buffer<Double>,
// 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<Double>, 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|<sup>3</sup>)<sup>3</sup>` for |x| &lt; 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<Pair<T, T>>): PiecewisePolynomial<T> {
// 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
// }
//}

58
kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt Normal file
View File

@ -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<T : Comparable<T>>(
override val algebra: Field<T>,
val bufferFactory: MutableBufferFactory<T>
) : PolynomialInterpolator<T> {
//TODO possibly optimize zeroed buffers
override fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T> = 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<T>(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)
}
}
}
}

54
kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt Normal file
View File

@ -0,0 +1,54 @@
package scientifik.kmath.interpolation
import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.Structure2D
interface XYPointSet<X, Y> {
val size: Int
val x: Buffer<X>
val y: Buffer<Y>
}
interface XYZPointSet<X, Y, Z> : XYPointSet<X, Y> {
val z: Buffer<Z>
}
internal fun <T : Comparable<T>> insureSorted(points: XYPointSet<T, *>) {
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<T>(val structure: Structure2D<T>, val column: Int) : Buffer<T> {
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<T> = sequence {
repeat(size) {
yield(get(it))
}
}.iterator()
}
class BufferXYPointSet<X, Y>(override val x: Buffer<X>, override val y: Buffer<Y>) : XYPointSet<X, Y> {
init {
require(x.size == y.size) { "Sizes of x and y buffers should be the same" }
}
override val size: Int
get() = x.size
}
fun <T> Structure2D<T>.asXYPointSet(): XYPointSet<T, T> {
require(shape[1] == 2) { "Structure second dimension should be of size 2" }
return object : XYPointSet<T, T> {
override val size: Int get() = this@asXYPointSet.shape[0]
override val x: Buffer<T> get() = NDStructureColumn(this@asXYPointSet, 0)
override val y: Buffer<T> get() = NDStructureColumn(this@asXYPointSet, 1)
}
}

27
kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt Normal file
View File

@ -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<Double> = 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))
}
}

9
kmath-geometry/build.gradle.kts Normal file
View File

@ -0,0 +1,9 @@
plugins {
id("scientifik.mpp")
}
kotlin.sourceSets.commonMain {
dependencies {
api(project(":kmath-core"))
}
}

58
kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean2DSpace.kt Normal file
View File

@ -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<Double>, Vector, SpaceElement<Vector2D, Vector2D, Euclidean2DSpace> {
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<Double> = 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<Vector2D> {
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
}

57
kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean3DSpace.kt Normal file
View File

@ -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<Double>, Vector, SpaceElement<Vector3D, Vector3D, Euclidean3DSpace> {
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<Double> = 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<Vector3D> {
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
}

17
kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/GeometrySpace.kt Normal file
View File

@ -0,0 +1,17 @@
package scientifik.kmath.geometry
import scientifik.kmath.operations.Space
interface Vector
interface GeometrySpace<V: Vector>: Space<V> {
/**
* L2 distance
*/
fun V.distanceTo(other: V): Double
/**
* Scalar product
*/
infix fun V.dot(other: V): Double
}

6
kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Line.kt Normal file
View File

@ -0,0 +1,6 @@
package scientifik.kmath.geometry
data class Line<V: Vector>(val base: V, val direction: V)
typealias Line2D = Line<Vector2D>
typealias Line3D = Line<Vector3D>

4
kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/ReferenceFrame.kt Normal file
View File

@ -0,0 +1,4 @@
package scientifik.kmath.geometry
interface ReferenceFrame {
}

1
kmath-histograms/build.gradle.kts
View File

@ -5,5 +5,6 @@ plugins {
kotlin.sourceSets.commonMain {
dependencies {
api(project(":kmath-core"))
api(project(":kmath-for-real"))
}
}

23
kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt
View File

@ -45,8 +45,7 @@ class RealHistogram(
private val values: NDStructure<LongCounter> = NDStructure.auto(strides) { LongCounter() }
//private val weight: NDStructure<DoubleCounter?> = ndStructure(strides){null}
private val weights: NDStructure<DoubleCounter> = 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<out Double>){
// val index = getIndex(point)
// values[index].increment()
// }
override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
if (weight != 1.0) TODO("Implement weighting")
val index = getIndex(point)
values[index].increment()
weights[index].add(weight)
}
override fun iterator(): Iterator<MultivariateBin<Double>> = values.elements().map { (index, value) ->
override fun iterator(): Iterator<MultivariateBin<Double>> = 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<Number> {
return NDStructure.auto(this.values.shape) { values[it].sum() }
fun values(): NDStructure<Number> {
return NDStructure.auto(values.shape) { values[it].sum() }
}
/**
* Sum of weights
*/
fun weights():NDStructure<Double>{
return NDStructure.auto(weights.shape) { weights[it].sum() }
}
companion object {

15
kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt
View File

@ -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()
@ -91,3 +95,12 @@ class ByteBufferMemory(
override fun writer(): MemoryWriter = writer
}
/**
* Use direct memory-mapped buffer from file to read something and close it afterwards.
*/
fun <R> 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()
}
}

8
kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt
View File

@ -58,9 +58,13 @@ fun <T : Any> Sampler<T>.sampleBuffer(
//creating temporary storage once
val tmp = ArrayList<T>(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] }
}
}

1
kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt
View File

@ -12,4 +12,3 @@ class RandomChain<out R>(val generator: RandomGenerator, private val gen: suspen
}
fun <R> RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain<R> = RandomChain(this, gen)
fun <R> RandomGenerator.flow(gen: suspend RandomGenerator.() -> R) = chain(gen).fork()

5
kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt
View File

@ -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<T>(val space: Space<T>) : ComposableStatistic<T, Pair<T, Int>, T> {
/**
* Non-composable median
*/
class Median<T>(comparator: Comparator<T>) : Statistic<T, T> {
class Median<T>(private val comparator: Comparator<T>) : Statistic<T, T> {
override suspend fun invoke(data: Buffer<T>): 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 {

0
kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/distributions.kt → kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/UniformDistribution.kt
View File

17
kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/SamplerTest.kt Normal file
View File

@ -0,0 +1,17 @@
package scientifik.kmath.prob
import kotlinx.coroutines.runBlocking
import kotlin.test.Test
class SamplerTest {
@Test
fun bufferSamplerTest(){
val sampler: Sampler<Double> =
BasicSampler { it.chain { nextDouble() } }
val data = sampler.sampleBuffer(RandomGenerator.default, 100)
runBlocking {
println(data.next())
}
}
}

4
kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt
View File

@ -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() {

10
kmath-viktor/build.gradle.kts Normal file
View File

@ -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")
}

20
kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt Normal file
View File

@ -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<Double> {
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<Double> {
return ViktorBuffer(flatArray.copy().flatten())
}
override fun iterator(): Iterator<Double> = flatArray.data.iterator()
}

86
kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt Normal file
View File

@ -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<Double> {
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<Pair<IntArray, Double>> {
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<Double, RealField, ViktorNDStructure> {
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()
}

16
settings.gradle.kts
View File

@ -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"
)