Fix tests #91

Merged
pklimai merged 5 commits from more-tests into dev 2020-05-06 20:03:00 +03:00
72 changed files with 2730 additions and 152 deletions
Showing only changes of commit f04eeac3b4 - Show all commits

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: [ ![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) 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 # KMath
Could be pronounced as `key-math`. 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. 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. * **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/) * **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 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. to submit a feature request if you want something to be done first.

View File

@ -1,8 +1,8 @@
plugins { 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 bintrayRepo by extra("scientifik")
val githubProject by extra("kmath") val githubProject by extra("kmath")

View File

@ -4,8 +4,8 @@ import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
plugins { plugins {
java java
kotlin("jvm") kotlin("jvm")
kotlin("plugin.allopen") version "1.3.60" kotlin("plugin.allopen") version "1.3.71"
id("kotlinx.benchmark") version "0.2.0-dev-5" id("kotlinx.benchmark") version "0.2.0-dev-7"
} }
configure<AllOpenExtension> { configure<AllOpenExtension> {
@ -13,10 +13,9 @@ configure<AllOpenExtension> {
} }
repositories { repositories {
maven("https://dl.bintray.com/kotlin/kotlin-eap")
maven("http://dl.bintray.com/kyonifer/maven") 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() mavenCentral()
} }
@ -29,12 +28,11 @@ dependencies {
implementation(project(":kmath-coroutines")) implementation(project(":kmath-coroutines"))
implementation(project(":kmath-commons")) implementation(project(":kmath-commons"))
implementation(project(":kmath-koma")) implementation(project(":kmath-koma"))
implementation(project(":kmath-viktor"))
implementation(project(":kmath-dimensions"))
implementation("com.kyonifer:koma-core-ejml:0.12") implementation("com.kyonifer:koma-core-ejml:0.12")
implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:${Scientifik.ioVersion}") 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")
implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-2")
"benchmarksCompile"(sourceSets.main.get().compileClasspath) "benchmarksCompile"(sourceSets.main.get().compileClasspath)
} }

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

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

View File

@ -4,7 +4,13 @@ import kotlinx.coroutines.GlobalScope
import scientifik.kmath.operations.RealField import scientifik.kmath.operations.RealField
import kotlin.system.measureTimeMillis 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 dim = 1000
val n = 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. //A generic boxing field. It should be used for objects, not primitives.
val genericField = NDField.boxing(RealField, dim, dim) val genericField = NDField.boxing(RealField, dim, dim)
measureAndPrint("Automatic field addition") {
val autoTime = measureTimeMillis {
autoField.run { autoField.run {
var res = one var res = one
repeat(n) { repeat(n) {
@ -25,18 +30,14 @@ fun main(args: Array<String>) {
} }
} }
println("Automatic field addition completed in $autoTime millis") measureAndPrint("Element addition"){
val elementTime = measureTimeMillis {
var res = genericField.one var res = genericField.one
repeat(n) { repeat(n) {
res += 1.0 res += 1.0
} }
} }
println("Element addition completed in $elementTime millis") measureAndPrint("Specialized addition") {
val specializedTime = measureTimeMillis {
specializedField.run { specializedField.run {
var res: NDBuffer<Double> = one var res: NDBuffer<Double> = one
repeat(n) { repeat(n) {
@ -45,10 +46,7 @@ fun main(args: Array<String>) {
} }
} }
println("Specialized addition completed in $specializedTime millis") measureAndPrint("Lazy addition") {
val lazyTime = measureTimeMillis {
val res = specializedField.one.mapAsync(GlobalScope) { val res = specializedField.one.mapAsync(GlobalScope) {
var c = 0.0 var c = 0.0
repeat(n) { repeat(n) {
@ -60,9 +58,7 @@ fun main(args: Array<String>) {
res.elements().forEach { it.second } res.elements().forEach { it.second }
} }
println("Lazy addition completed in $lazyTime millis") measureAndPrint("Generic addition") {
val genericTime = measureTimeMillis {
//genericField.run(action) //genericField.run(action)
genericField.run { genericField.run {
var res: NDBuffer<Double> = one var res: NDBuffer<Double> = one
@ -72,6 +68,4 @@ fun main(args: Array<String>) {
} }
} }
println("Generic addition completed in $genericTime millis")
} }

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

Binary file not shown.

View File

@ -1,5 +1,5 @@
distributionBase=GRADLE_USER_HOME distributionBase=GRADLE_USER_HOME
distributionPath=wrapper/dists 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 zipStoreBase=GRADLE_USER_HOME
zipStorePath=wrapper/dists zipStorePath=wrapper/dists

3
gradlew.bat vendored
View File

@ -29,6 +29,9 @@ if "%DIRNAME%" == "" set DIRNAME=.
set APP_BASE_NAME=%~n0 set APP_BASE_NAME=%~n0
set APP_HOME=%DIRNAME% 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. @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" set DEFAULT_JVM_OPTS="-Xmx64m" "-Xms64m"

View File

@ -8,7 +8,6 @@ dependencies {
api(project(":kmath-core")) api(project(":kmath-core"))
api(project(":kmath-coroutines")) api(project(":kmath-coroutines"))
api(project(":kmath-prob")) api(project(":kmath-prob"))
api(project(":kmath-functions"))
api("org.apache.commons:commons-math3:3.6.1") api("org.apache.commons:commons-math3:3.6.1")
testImplementation("org.jetbrains.kotlin:kotlin-test")
testImplementation("org.jetbrains.kotlin:kotlin-test-junit")
} }

View File

@ -1,7 +1,7 @@
package scientifik.kmath.commons.expressions package scientifik.kmath.commons.expressions
import org.junit.Test
import scientifik.kmath.expressions.invoke import scientifik.kmath.expressions.invoke
import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
inline fun <R> diff(order: Int, vararg parameters: Pair<String, Double>, block: DerivativeStructureField.() -> R) = inline fun <R> diff(order: Int, vararg parameters: Pair<String, Double>, block: DerivativeStructureField.() -> R) =

View File

@ -41,16 +41,6 @@ fun <T : Any> Structure2D.Companion.square(vararg elements: T): FeaturedMatrix<T
return BufferMatrix(size, size, buffer) 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() val Matrix<*>.features get() = (this as? FeaturedMatrix)?.features?: emptySet()
/** /**

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)

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
}

View File

@ -6,14 +6,14 @@ annotation class KMathContext
/** /**
* Marker interface for any algebra * 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 * Space-like operations without neutral element
*/ */
interface SpaceOperations<T> : Algebra { interface SpaceOperations<T> : Algebra<T> {
/** /**
* Addition operation for two context elements * Addition operation for two context elements
*/ */

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: 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> 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
}

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
}

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 * 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 power(arg: T, pow: Number): T
fun sqrt(arg: T) = power(arg, 0.5) 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 */ /* Exponential */
interface ExponentialOperations<T> { interface ExponentialOperations<T>: Algebra<T> {
fun exp(arg: T): T fun exp(arg: T): T
fun ln(arg: T): T fun ln(arg: T): T
} }

View File

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

View File

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

View File

@ -1,26 +1,13 @@
package scientifik.kmath.linear package scientifik.kmath.linear
import scientifik.kmath.structures.Matrix import scientifik.kmath.structures.Matrix
import scientifik.kmath.structures.NDStructure
import scientifik.kmath.structures.as2D
import kotlin.test.Test import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
class MatrixTest { 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 @Test
fun testTranspose() { fun testTranspose() {
val matrix = MatrixContext.real.one(3, 3) val matrix = MatrixContext.real.one(3, 3)
@ -28,21 +15,6 @@ class MatrixTest {
assertEquals(matrix, transposed) 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 @Test
fun testBuilder() { fun testBuilder() {
val matrix = Matrix.build<Double>(2, 3)( val matrix = Matrix.build<Double>(2, 3)(
@ -74,4 +46,20 @@ class MatrixTest {
val toTenthPower = transitionMatrix pow 10 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])
}
}
} }

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

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

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

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

View File

@ -16,7 +16,9 @@
package scientifik.kmath.chains package scientifik.kmath.chains
import kotlinx.coroutines.InternalCoroutinesApi
import kotlinx.coroutines.flow.Flow import kotlinx.coroutines.flow.Flow
import kotlinx.coroutines.flow.FlowCollector
import kotlinx.coroutines.sync.Mutex import kotlinx.coroutines.sync.Mutex
import kotlinx.coroutines.sync.withLock import kotlinx.coroutines.sync.withLock
@ -25,7 +27,7 @@ import kotlinx.coroutines.sync.withLock
* A not-necessary-Markov chain of some type * A not-necessary-Markov chain of some type
* @param R - the chain element type * @param R - the chain element type
*/ */
interface Chain<out R> { interface Chain<out R>: Flow<R> {
/** /**
* Generate next value, changing state if needed * Generate next value, changing state if needed
*/ */
@ -36,14 +38,15 @@ interface Chain<out R> {
*/ */
fun fork(): Chain<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 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> Iterator<T>.asChain(): Chain<T> = SimpleChain { next() }
fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain() 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) 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 * Map the whole chain
*/ */

View File

@ -3,11 +3,12 @@ package scientifik.kmath.streaming
import kotlinx.coroutines.* import kotlinx.coroutines.*
import kotlinx.coroutines.flow.asFlow import kotlinx.coroutines.flow.asFlow
import kotlinx.coroutines.flow.collect import kotlinx.coroutines.flow.collect
import org.junit.Test import org.junit.jupiter.api.Timeout
import scientifik.kmath.coroutines.async import scientifik.kmath.coroutines.async
import scientifik.kmath.coroutines.collect import scientifik.kmath.coroutines.collect
import scientifik.kmath.coroutines.mapParallel import scientifik.kmath.coroutines.mapParallel
import java.util.concurrent.Executors import java.util.concurrent.Executors
import kotlin.test.Test
@ExperimentalCoroutinesApi @ExperimentalCoroutinesApi
@ -17,7 +18,8 @@ class BufferFlowTest {
val dispatcher = Executors.newFixedThreadPool(4).asCoroutineDispatcher() val dispatcher = Executors.newFixedThreadPool(4).asCoroutineDispatcher()
@Test(timeout = 2000) @Test
@Timeout(2000)
fun map() { fun map() {
runBlocking { runBlocking {
(1..20).asFlow().mapParallel( dispatcher) { (1..20).asFlow().mapParallel( dispatcher) {
@ -31,7 +33,8 @@ class BufferFlowTest {
} }
} }
@Test(timeout = 2000) @Test
@Timeout(2000)
fun async() { fun async() {
runBlocking { runBlocking {
(1..20).asFlow().async(dispatcher) { (1..20).asFlow().async(dispatcher) {

View File

@ -2,8 +2,8 @@ package scientifik.kmath.streaming
import kotlinx.coroutines.flow.* import kotlinx.coroutines.flow.*
import kotlinx.coroutines.runBlocking import kotlinx.coroutines.runBlocking
import org.junit.Test
import scientifik.kmath.structures.asSequence import scientifik.kmath.structures.asSequence
import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
class RingBufferTest { class RingBufferTest {

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

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
}

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

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

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

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

View File

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

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

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

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

View File

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

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

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

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
}

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

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

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
// }
//}

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

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

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

View File

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

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
}

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
}

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
}

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>

View File

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

View File

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

View File

@ -45,8 +45,7 @@ class RealHistogram(
private val values: NDStructure<LongCounter> = NDStructure.auto(strides) { LongCounter() } 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 override val dimension: Int get() = lower.size
@ -102,21 +101,33 @@ class RealHistogram(
return MultivariateBin(getDef(index), getValue(index)) 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) { override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
if (weight != 1.0) TODO("Implement weighting")
val index = getIndex(point) val index = getIndex(point)
values[index].increment() 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()) MultivariateBin(getDef(index), value.sum())
}.iterator() }.iterator()
/** /**
* Convert this histogram into NDStructure containing bin values but not bin descriptions * Convert this histogram into NDStructure containing bin values but not bin descriptions
*/ */
fun asNDStructure(): NDStructure<Number> { fun values(): NDStructure<Number> {
return NDStructure.auto(this.values.shape) { values[it].sum() } 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 { companion object {

View File

@ -1,6 +1,10 @@
package scientifik.memory package scientifik.memory
import java.nio.ByteBuffer 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) return ByteBufferMemory(buffer)
} }
class ByteBufferMemory( private class ByteBufferMemory(
val buffer: ByteBuffer, val buffer: ByteBuffer,
val startOffset: Int = 0, val startOffset: Int = 0,
override val size: Int = buffer.limit() override val size: Int = buffer.limit()
@ -91,3 +95,12 @@ class ByteBufferMemory(
override fun writer(): MemoryWriter = writer 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()
}
}

View File

@ -58,9 +58,13 @@ fun <T : Any> Sampler<T>.sampleBuffer(
//creating temporary storage once //creating temporary storage once
val tmp = ArrayList<T>(size) val tmp = ArrayList<T>(size)
return sample(generator).collect { chain -> return sample(generator).collect { chain ->
for (i in tmp.indices) { //clear list from previous run
tmp[i] = chain.next() tmp.clear()
//Fill list
repeat(size){
tmp.add(chain.next())
} }
//return new buffer with elements from tmp
bufferFactory(size) { tmp[it] } bufferFactory(size) { tmp[it] }
} }
} }

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.chain(gen: suspend RandomGenerator.() -> R): RandomChain<R> = RandomChain(this, gen)
fun <R> RandomGenerator.flow(gen: suspend RandomGenerator.() -> R) = chain(gen).fork()

View File

@ -11,6 +11,7 @@ import scientifik.kmath.coroutines.mapParallel
import scientifik.kmath.operations.* import scientifik.kmath.operations.*
import scientifik.kmath.structures.Buffer import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.asIterable import scientifik.kmath.structures.asIterable
import scientifik.kmath.structures.asSequence
/** /**
* A function, that transforms a buffer of random quantities to some resulting value * 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 * 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 { 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 { companion object {

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

View File

@ -3,7 +3,7 @@ package scientifik.kmath.prob
import kotlinx.coroutines.flow.drop import kotlinx.coroutines.flow.drop
import kotlinx.coroutines.flow.first import kotlinx.coroutines.flow.first
import kotlinx.coroutines.runBlocking import kotlinx.coroutines.runBlocking
import scientifik.kmath.chains.flow
import scientifik.kmath.streaming.chunked import scientifik.kmath.streaming.chunked
import kotlin.test.Test import kotlin.test.Test
@ -13,7 +13,7 @@ class StatisticTest {
//Create a stateless chain from generator. //Create a stateless chain from generator.
val data = generator.chain { nextDouble() } val data = generator.chain { nextDouble() }
//Convert a chaint to Flow and break it into chunks. //Convert a chaint to Flow and break it into chunks.
val chunked = data.flow().chunked(1000) val chunked = data.chunked(1000)
@Test @Test
fun testParallelMean() { fun testParallelMean() {

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

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

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

View File

@ -1,10 +1,10 @@
pluginManagement { pluginManagement {
plugins { plugins {
id("scientifik.mpp") version "0.2.5" id("scientifik.mpp") version "0.4.1"
id("scientifik.jvm") version "0.2.5" id("scientifik.jvm") version "0.4.1"
id("scientifik.atomic") version "0.2.5" id("scientifik.atomic") version "0.4.1"
id("scientifik.publish") version "0.2.5" id("scientifik.publish") version "0.4.1"
} }
repositories { repositories {
@ -13,13 +13,14 @@ pluginManagement {
gradlePluginPortal() gradlePluginPortal()
maven("https://dl.bintray.com/kotlin/kotlin-eap") maven("https://dl.bintray.com/kotlin/kotlin-eap")
maven("https://dl.bintray.com/mipt-npm/scientifik") maven("https://dl.bintray.com/mipt-npm/scientifik")
maven("https://dl.bintray.com/mipt-npm/dev")
maven("https://dl.bintray.com/kotlin/kotlinx") maven("https://dl.bintray.com/kotlin/kotlinx")
} }
resolutionStrategy { resolutionStrategy {
eachPlugin { eachPlugin {
when (requested.id.id) { 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( include(
":kmath-memory", ":kmath-memory",
":kmath-core", ":kmath-core",
":kmath-functions",
// ":kmath-io", // ":kmath-io",
":kmath-coroutines", ":kmath-coroutines",
":kmath-histograms", ":kmath-histograms",
":kmath-commons", ":kmath-commons",
":kmath-viktor",
":kmath-koma", ":kmath-koma",
":kmath-prob", ":kmath-prob",
":kmath-io", ":kmath-io",
":kmath-dimensions",
":kmath-for-real",
":kmath-geometry",
":examples" ":examples"
) )