diff --git a/README.md b/README.md index bcdda2030..9dd62a21e 100644 --- a/README.md +++ b/README.md @@ -1,53 +1,54 @@ # KMath -Kotlin MATHematics library is intended as a kotlin based analog of numpy python library. Contrary 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. ## Features * **Algebra** - * Mathematical operation entities like rings, spaces and fields with (**TODO** add example to wiki) - * Basic linear algebra operations (summs products, etc) backed by `Space` API. - * Complex numbers backed by `Field` API (meaning that they will be useable in any structures like vectors and NDArrays). - * [In progress] advanced linear algebra operations like matrix inversions. -* **Array-like structures** Full support of numpy-like ndarray including mixed arithmetic operations and function operations -on arrays and numbers just like it works in python (with benefit of static type checking). + * Algebraic structures like rings, spaces and field (**TODO** add example to wiki) + * Basic linear algebra operations (sums, products, etc.), backed by the `Space` API. + * Complex numbers backed by the `Field` API (meaning that they will be usable in any structure like vectors and N-dimensional arrays). + * [In progress] advanced linear algebra operations like matrix inversion and LU decomposition. +* **Array-like structures** Full support of [numpy-like ndarrays](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ndarray.html) including mixed arithmetic operations and function operations over arrays and numbers just like in Python (with the added benefit of static type checking). -* **Expressions** Expressions are one of the ultimate goals of kmath. It is planned to be able to write some mathematical -expression once an then apply it to different types of objects by providing different context. Expressions could be used -for a wide variety of purposes from high performance calculations to code generation. +* **Expressions** Expressions are one of the ultimate goals of KMath. By writing a single mathematical expression +once, users will be able to apply different types of objects to the expression by providing a context. Exceptions +can be used for a wide variety of purposes from high performance calculations to code generation. ## Planned features * **Common mathematics** 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 suite kotlin programming paradigm. There is no fixed priority list for that. Feel free -to submit a future request if you want something to be done first. +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. * **Messaging** A mathematical notation to support multi-language and multi-node communication for mathematical tasks. ## Multi-platform support -KMath is developed as a multi-platform library, which means that most of interfaces are declared in common module. -Implementation is also done in common module wherever it is possible. In some cases features are delegated to -platform even if they could be done in common module because of platform performance optimization. -Currently the main focus of development is the JVM platform, contribution of implementations for Kotlin - Native and -Kotlin - JS is welcome. + +KMath is developed as a multi-platform library, which means that most of interfaces are declared in the [common module](kmath-core/src/commonMain). +Implementation is also done in the common module wherever possible. In some cases, features are delegated to +platform-specific implementations even if they could be done in the common module for performance reasons. +Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome. ## Performance -The calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve -both performance and flexibility. We expect to firstly focus on creating convenient universal API and then work on -increasing performance for specific cases. We expect the worst KMath performance still be better than natural python, -but worse than optimized native/scipy (mostly due to boxing operations on primitive numbers). The best performance -of optimized parts should be better than scipy. + +Calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve +both performance and flexibility. We expect to focus on creating convenient universal API first and then work on +increasing performance for specific cases. We expect the worst KMath benchmarks will perform better than native Python, +but worse than optimized native/SciPy (mostly due to boxing operations on primitive numbers). The best performance +of optimized parts should be better than SciPy. ## Releases -The project is currently in pre-release stage. Nightly builds could be used by adding additional repository to (groovy) gradle config: +The project is currently in pre-release stage. Nightly builds can be used by adding an additional repository to the Gradle config like so: + ```groovy repositories { maven { url = "http://npm.mipt.ru:8081/artifactory/gradle-dev" } mavenCentral() } ``` -or for kotlin gradle dsl: + +or for the Gradle Kotlin DSL: ```kotlin repositories { @@ -56,16 +57,20 @@ repositories { } ``` -Then use regular dependency like +Then use a regular dependency like so: + ```groovy compile(group: 'scientifik', name: 'kmath-core', version: '0.0.1-SNAPSHOT') ``` -or in kotlin + +or in the Gradle Kotlin DSL: + ```kotlin compile(group = "scientifik", name = "kmath-core", version = "0.0.1-SNAPSHOT") ``` -Work builds could be obtained with [![](https://jitpack.io/v/altavir/kmath.svg)](https://jitpack.io/#altavir/kmath). +Working builds can be obtained here: [![](https://jitpack.io/v/altavir/kmath.svg)](https://jitpack.io/#altavir/kmath). ## Contributing -The project requires a lot of additional work. Please fill free to contribute in any way and propose new features. + +The project requires a lot of additional work. Please fill free to contribute in any way and propose new features. \ No newline at end of file diff --git a/benchmarks/build.gradle b/benchmarks/build.gradle new file mode 100644 index 000000000..989590397 --- /dev/null +++ b/benchmarks/build.gradle @@ -0,0 +1,39 @@ +plugins { + id "java" + id "me.champeau.gradle.jmh" version "0.4.8" + id 'org.jetbrains.kotlin.jvm' +} + +repositories { + maven { url 'https://dl.bintray.com/kotlin/kotlin-eap' } + maven{ url "http://dl.bintray.com/kyonifer/maven"} + mavenCentral() +} + +dependencies { + implementation project(":kmath-core") + implementation project(":kmath-coroutines") + implementation project(":kmath-commons") + implementation project(":kmath-koma") + implementation group: "com.kyonifer", name:"koma-core-ejml", version: "0.12" + implementation "org.jetbrains.kotlinx:kotlinx-io-jvm:0.1.5" + //compile "org.jetbrains.kotlin:kotlin-stdlib-jdk8" + //jmh project(':kmath-core') +} + +jmh{ + warmupIterations = 1 +} + +jmhClasses.dependsOn(compileKotlin) + +compileKotlin { + kotlinOptions { + jvmTarget = "1.8" + } +} +compileTestKotlin { + kotlinOptions { + jvmTarget = "1.8" + } +} \ No newline at end of file diff --git a/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt b/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt new file mode 100644 index 000000000..393d7b06a --- /dev/null +++ b/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt @@ -0,0 +1,48 @@ +package scientifik.kmath.structures + +import org.openjdk.jmh.annotations.Benchmark +import org.openjdk.jmh.annotations.Scope +import org.openjdk.jmh.annotations.State +import java.nio.IntBuffer + + +@State(Scope.Benchmark) +open class ArrayBenchmark { + + @Benchmark + fun benchmarkArrayRead() { + var res = 0 + for (i in 1..size) { + res += array[size - i] + } + } + + @Benchmark + fun benchmarkBufferRead() { + var res = 0 + for (i in 1..size) { + res += arrayBuffer.get(size - i) + } + } + + @Benchmark + fun nativeBufferRead() { + var res = 0 + for (i in 1..size) { + res += nativeBuffer.get(size - i) + } + } + + companion object { + val size = 1000 + + val array = IntArray(size) { it } + val arrayBuffer = IntBuffer.wrap(array) + val nativeBuffer = IntBuffer.allocate(size).also { + for (i in 0 until size) { + it.put(i, i) + } + + } + } +} \ No newline at end of file diff --git a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt b/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt similarity index 67% rename from kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt rename to benchmarks/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt index fd27ae3e0..5ca05d451 100644 --- a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt +++ b/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt @@ -1,10 +1,10 @@ package scientifik.kmath.structures -import org.openjdk.jmh.annotations.* +import org.openjdk.jmh.annotations.Benchmark +import org.openjdk.jmh.annotations.Scope +import org.openjdk.jmh.annotations.State import scientifik.kmath.operations.Complex -@Warmup(iterations = 1) -@Measurement(iterations = 5) @State(Scope.Benchmark) open class BufferBenchmark { @@ -22,17 +22,17 @@ open class BufferBenchmark { @Benchmark fun complexBufferReadWrite() { - val buffer = Complex.createBuffer(size/2) - (0 until size/2).forEach { + val buffer = MutableBuffer.complex(size / 2) + (0 until size / 2).forEach { buffer[it] = Complex(it.toDouble(), -it.toDouble()) } - (0 until size/2).forEach { + (0 until size / 2).forEach { buffer[it] } } companion object { - const val size = 1000 + const val size = 100 } } \ No newline at end of file diff --git a/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt b/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt new file mode 100644 index 000000000..0cef4bdb4 --- /dev/null +++ b/benchmarks/src/jmh/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt @@ -0,0 +1,69 @@ +package scientifik.kmath.structures + +import org.openjdk.jmh.annotations.Benchmark +import scientifik.kmath.operations.RealField + +open class NDFieldBenchmark { + + @Benchmark + fun autoFieldAdd() { + bufferedField.run { + var res: NDBuffer = one + repeat(n) { + res += one + } + } + } + + @Benchmark + fun autoElementAdd() { + var res = genericField.one + repeat(n) { + res += 1.0 + } + } + + @Benchmark + fun specializedFieldAdd() { + specializedField.run { + var res: NDBuffer = one + repeat(n) { + res += 1.0 + } + } + } + + + @Benchmark + fun lazyFieldAdd() { + lazyNDField.run { + var res = one + repeat(n) { + res += one + } + + res.elements().sumByDouble { it.second } + } + } + + + @Benchmark + fun boxingFieldAdd() { + genericField.run { + var res: NDBuffer = one + repeat(n) { + res += one + } + } + } + + companion object { + val dim = 1000 + val n = 100 + + val bufferedField = NDField.auto(RealField, intArrayOf(dim, dim)) + val specializedField = NDField.real(intArrayOf(dim, dim)) + val genericField = NDField.buffered(intArrayOf(dim, dim), RealField) + val lazyNDField = NDField.lazy(intArrayOf(dim, dim), RealField) + } +} \ No newline at end of file diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt new file mode 100644 index 000000000..d57631eba --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt @@ -0,0 +1,61 @@ +package scientifik.kmath.linear + +import koma.matrix.ejml.EJMLMatrixFactory +import kotlin.random.Random +import kotlin.system.measureTimeMillis + +fun main() { + val random = Random(12224) + val dim = 100 + //creating invertible matrix + val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } + val matrix = l dot u + + val n = 500 // iterations + + val solver = LUSolver.real + + repeat(50) { + val res = solver.inverse(matrix) + } + + val inverseTime = measureTimeMillis { + repeat(n) { + val res = solver.inverse(matrix) + } + } + + println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis") + + //commons-math + + val cmContext = CMMatrixContext + + val commonsTime = measureTimeMillis { + cmContext.run { + val cm = matrix.toCM() //avoid overhead on conversion + repeat(n) { + val res = inverse(cm) + } + } + } + + + println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis") + + //koma-ejml + + val komaContext = KomaMatrixContext(EJMLMatrixFactory()) + + val komaTime = measureTimeMillis { + komaContext.run { + val km = matrix.toKoma() //avoid overhead on conversion + repeat(n) { + val res = inverse(km) + } + } + } + + println("[koma-ejml] Inversion of $n matrices $dim x $dim finished in $komaTime millis") +} \ No newline at end of file diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt new file mode 100644 index 000000000..0b3b7c275 --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt @@ -0,0 +1,45 @@ +package scientifik.kmath.linear + +import koma.matrix.ejml.EJMLMatrixFactory +import kotlin.random.Random +import kotlin.system.measureTimeMillis + +fun main() { + val random = Random(12224) + val dim = 1000 + //creating invertible matrix + val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + +// //warmup +// matrix1 dot matrix2 + + CMMatrixContext.run { + val cmMatrix1 = matrix1.toCM() + val cmMatrix2 = matrix2.toCM() + + val cmTime = measureTimeMillis { + cmMatrix1 dot cmMatrix2 + } + + println("CM implementation time: $cmTime") + } + + + KomaMatrixContext(EJMLMatrixFactory()).run { + val komaMatrix1 = matrix1.toKoma() + val komaMatrix2 = matrix2.toKoma() + + val komaTime = measureTimeMillis { + komaMatrix1 dot komaMatrix2 + } + + println("Koma-ejml implementation time: $komaTime") + } + + val genericTime = measureTimeMillis { + val res = matrix1 dot matrix2 + } + + println("Generic implementation time: $genericTime") +} \ No newline at end of file diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt new file mode 100644 index 000000000..46b1b56dd --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt @@ -0,0 +1,34 @@ +package scientifik.kmath.structures + +import kotlin.system.measureTimeMillis + +fun main() { + val dim = 1000 + val n = 1000 + + val realField = NDField.real(dim, dim) + val complexField = NDField.complex(dim, dim) + + + val realTime = measureTimeMillis { + realField.run { + var res: NDBuffer = one + repeat(n) { + res += 1.0 + } + } + } + + println("Real addition completed in $realTime millis") + + val complexTime = measureTimeMillis { + complexField.run { + var res: ComplexNDElement = one + repeat(n) { + res += 1.0 + } + } + } + + println("Complex addition completed in $complexTime millis") +} \ No newline at end of file diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt new file mode 100644 index 000000000..53a23c2bc --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt @@ -0,0 +1,77 @@ +package scientifik.kmath.structures + +import kotlinx.coroutines.GlobalScope +import scientifik.kmath.operations.RealField +import kotlin.system.measureTimeMillis + +fun main(args: Array) { + val dim = 1000 + val n = 1000 + + // automatically build context most suited for given type. + val autoField = NDField.auto(RealField, dim, dim) + // specialized nd-field for Double. It works as generic Double field as well + val specializedField = NDField.real(dim, dim) + //A generic boxing field. It should be used for objects, not primitives. + val genericField = NDField.buffered(intArrayOf(dim, dim), RealField) + + + val autoTime = measureTimeMillis { + autoField.run { + var res = one + repeat(n) { + res += 1.0 + } + } + } + + println("Automatic field addition completed in $autoTime millis") + + val elementTime = measureTimeMillis { + var res = genericField.one + repeat(n) { + res += 1.0 + } + } + + println("Element addition completed in $elementTime millis") + + val specializedTime = measureTimeMillis { + specializedField.run { + var res: NDBuffer = one + repeat(n) { + res += 1.0 + } + } + } + + println("Specialized addition completed in $specializedTime millis") + + + val lazyTime = measureTimeMillis { + val res = specializedField.one.mapAsync(GlobalScope) { + var c = 0.0 + repeat(n) { + c += 1.0 + } + c + } + + res.elements().forEach { it.second } + } + + println("Lazy addition completed in $lazyTime millis") + + val genericTime = measureTimeMillis { + //genericField.run(action) + genericField.run { + var res: NDBuffer = one + repeat(n) { + res += 1.0 + } + } + } + + println("Generic addition completed in $genericTime millis") + +} \ No newline at end of file diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureReadBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureReadBenchmark.kt new file mode 100644 index 000000000..ecfb4ab20 --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureReadBenchmark.kt @@ -0,0 +1,36 @@ +package scientifik.kmath.structures + +import kotlin.system.measureTimeMillis + +fun main(args: Array) { + val n = 6000 + + val array = DoubleArray(n * n) { 1.0 } + val buffer = DoubleBuffer(array) + val strides = DefaultStrides(intArrayOf(n, n)) + + val structure = BufferNDStructure(strides, buffer) + + measureTimeMillis { + var res: Double = 0.0 + strides.indices().forEach { res = structure[it] } + } // warmup + + val time1 = measureTimeMillis { + var res: Double = 0.0 + strides.indices().forEach { res = structure[it] } + } + println("Structure reading finished in $time1 millis") + + val time2 = measureTimeMillis { + var res: Double = 0.0 + strides.indices().forEach { res = buffer[strides.offset(it)] } + } + println("Buffer reading finished in $time2 millis") + + val time3 = measureTimeMillis { + var res: Double = 0.0 + strides.indices().forEach { res = array[strides.offset(it)] } + } + println("Array reading finished in $time3 millis") +} \ No newline at end of file diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt new file mode 100644 index 000000000..e5147b941 --- /dev/null +++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt @@ -0,0 +1,39 @@ +package scientifik.kmath.structures + +import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory +import kotlin.system.measureTimeMillis + + +fun main(args: Array) { + + val n = 6000 + + val structure = NDStructure.build(intArrayOf(n, n), DoubleBufferFactory) { 1.0 } + + structure.mapToBuffer { it + 1 } // warm-up + + val time1 = measureTimeMillis { + val res = structure.mapToBuffer { it + 1 } + } + println("Structure mapping finished in $time1 millis") + + val array = DoubleArray(n * n) { 1.0 } + + val time2 = measureTimeMillis { + val target = DoubleArray(n * n) + val res = array.forEachIndexed { index, value -> + target[index] = value + 1 + } + } + println("Array mapping finished in $time2 millis") + + val buffer = DoubleBuffer(DoubleArray(n * n) { 1.0 }) + + val time3 = measureTimeMillis { + val target = DoubleBuffer(DoubleArray(n * n)) + val res = array.forEachIndexed { index, value -> + target[index] = value + 1 + } + } + println("Buffer mapping finished in $time3 millis") +} \ No newline at end of file diff --git a/build.gradle.kts b/build.gradle.kts index 5b95be71a..dadb99422 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,13 +1,13 @@ -buildscript { - extra["kotlinVersion"] = "1.3.11" - extra["ioVersion"] = "0.1.2-dev-2" - extra["coroutinesVersion"] = "1.0.1" +import org.jetbrains.kotlin.gradle.dsl.KotlinMultiplatformExtension - val kotlinVersion: String by extra - val ioVersion: String by extra - val coroutinesVersion: String by extra +buildscript { + val kotlinVersion: String by rootProject.extra("1.3.21") + val ioVersion: String by rootProject.extra("0.1.5") + val coroutinesVersion: String by rootProject.extra("1.1.1") + val atomicfuVersion: String by rootProject.extra("0.12.1") repositories { + //maven("https://dl.bintray.com/kotlin/kotlin-eap") jcenter() } @@ -18,18 +18,39 @@ buildscript { } plugins { - id("com.jfrog.artifactory") version "4.8.1" apply false -// id("org.jetbrains.kotlin.multiplatform") apply false + id("com.jfrog.artifactory") version "4.9.1" apply false } allprojects { - apply(plugin = "maven-publish") - apply(plugin = "com.jfrog.artifactory") + if (project.name.startsWith("kmath")) { + apply(plugin = "maven-publish") + apply(plugin = "com.jfrog.artifactory") + } group = "scientifik" - version = "0.0.2-dev-1" + version = "0.0.3-dev" + + repositories { + //maven("https://dl.bintray.com/kotlin/kotlin-eap") + jcenter() + } + + extensions.findByType()?.apply { + jvm { + compilations.all { + kotlinOptions { + jvmTarget = "1.8" + } + } + } + targets.all { + sourceSets.all { + languageSettings.progressiveMode = true + } + } + } } -if(file("artifactory.gradle").exists()){ +if (file("artifactory.gradle").exists()) { apply(from = "artifactory.gradle") -} \ No newline at end of file +} diff --git a/doc/algebra.md b/doc/algebra.md new file mode 100644 index 000000000..2888f9484 --- /dev/null +++ b/doc/algebra.md @@ -0,0 +1,70 @@ +# Algebra and algebra elements + +The mathematical operations in `kmath` are generally separated from mathematical objects. +This means that in order to perform an operation, say `+`, one needs two objects of a type `T` and +and algebra context which defines appropriate operation, say `Space`. Next one needs to run actual operation +in the context: + +```kotlin +val a: T +val b: T +val space: Space + +val c = space.run{a + b} +``` + +From the first glance, this distinction seems to be a needless complication, but in fact one needs +to remember that in mathematics, one could define different operations on the same objects. For example, +one could use different types of geometry for vectors. + +## Algebra hierarchy + +Mathematical contexts have the following hierarchy: + +**Space** <- **Ring** <- **Field** + +All classes follow abstract mathematical constructs. +[Space](http://mathworld.wolfram.com/Space.html) defines `zero` element, addition operation and multiplication by constant, +[Ring](http://mathworld.wolfram.com/Ring.html) adds multiplication and unit `one` element, +[Field](http://mathworld.wolfram.com/Field.html) adds division operation. + +Typical case of `Field` is the `RealField` which works on doubles. And typical case of `Space` is a `VectorSpace`. + +In some cases algebra context could hold additional operation like `exp` or `sin`, in this case it inherits appropriate +interface. Also a context could have an operation which produces an element outside of its context. For example +`Matrix` `dot` operation produces a matrix with new dimensions which could not be compatible with initial matrix in +terms of linear operations. + +## Algebra element + +In order to achieve more familiar behavior (where you apply operations directly to mathematica objects), without involving contexts +`kmath` introduces special type objects called `MathElement`. A `MathElement` is basically some object coupled to +a mathematical context. For example `Complex` is the pair of real numbers representing real and imaginary parts, +but it also holds reference to the `ComplexField` singleton which allows to perform direct operations on `Complex` +numbers without explicit involving the context like: + +```kotlin + val c1 = Complex(1.0, 1.0) + val c2 = Complex(1.0, -1.0) + val c3 = c1 + c2 + 3.0.toComplex() + //or with field notation: + val c4 = ComplexField.run{c1 + i - 2.0} +``` + +Both notations have their pros and cons. + +The hierarchy for algebra elements follows the hierarchy for the corresponding algebra. + +**MathElement** <- **SpaceElement** <- **RingElement** <- **FieldElement** + +**MathElement** is the generic common ancestor of the class with context. + +One important distinction between algebra elements and algebra contexts is that algebra element has three type parameters: + +1. The type of elements, field operates on. +2. The self-type of the element returned from operation (must be algebra element). +3. The type of the algebra over first type-parameter. + +The middle type is needed in case algebra members do not store context. For example, it is not possible to add +a context to regular `Double`. The element performs automatic conversions from context types and back. +One should used context operations in all important places. The performance of element operations is not guaranteed. diff --git a/doc/buffers.md b/doc/buffers.md new file mode 100644 index 000000000..6208deaff --- /dev/null +++ b/doc/buffers.md @@ -0,0 +1 @@ +**TODO** \ No newline at end of file diff --git a/doc/contexts.md b/doc/contexts.md new file mode 100644 index 000000000..58b198046 --- /dev/null +++ b/doc/contexts.md @@ -0,0 +1,73 @@ +# Context-oriented mathematics + +## The problem + +A known problem for implementing mathematics in statically-typed languages (but not only in them) is that different +sets of mathematical operators can be defined on the same mathematical objects. Sometimes there is no single way to +treat some operations, including basic arithmetic operations, on a Java/Kotlin `Number`. Sometimes there are different ways to +define the same structure, such as Euclidean and elliptic geometry vector spaces over real vectors. Another problem arises when +one wants to add some kind of behavior to an existing entity. In dynamic languages those problems are usually solved +by adding dynamic context-specific behaviors at runtime, but this solution has a lot of drawbacks. + +## Context-oriented approach + +One possible solution to these problems is to divorce numerical representations from behaviors. +For example in Kotlin one can define a separate class which represents some entity without any operations, +ex. a complex number: + +```kotlin +data class Complex(val re: Double, val im: Double) +``` + +And then to define a separate class or singleton, representing an operation on those complex numbers: + +```kotlin +object ComplexOperations { + operator fun Complex.plus(other: Complex) = Complex(re + other.re, im + other.im) + operator fun Complex.minus(other: Complex) = Complex(re - other.re, im - other.im) +} +``` + +In Java, applying such external operations could be very cumbersome, but Kotlin has a unique feature which allows us +implement this naturally: [extensions with receivers](https://kotlinlang.org/docs/reference/extensions.html#extension-functions). +In Kotlin, an operation on complex number could be implemented as: + +```kotlin +with(ComplexOperations) { c1 + c2 - c3 } +``` + +Kotlin also allows the creation of functions with receivers: + +```kotlin +fun ComplexOperations.doSomethingWithComplex(c1: Complex, c2: Complex, c3: Complex) = c1 + c2 - c3 + +ComplexOperations.doComethingWithComplex(c1, c2, c3) +``` + +In fact, whole parts of a program may be run within a mathematical context or even multiple nested contexts. + +In KMath, contexts are not only responsible for operations, but also for raw object creation and advanced features. + +## Other possibilities + +### Type classes + +An obvious candidate to get more or less the same functionality is the type class, which allows one to bind a behavior to +a specific type without modifying the type itself. On the plus side, type classes do not require explicit context +declaration, so the code looks cleaner. On the minus side, if there are different sets of behaviors for the same types, +it is impossible to combine them into one module. Also, unlike type classes, context can have parameters or even +state. For example in KMath, sizes and strides for `NDElement` or `Matrix` could be moved to context to optimize +performance in case of a large amount of structures. + +### Wildcard imports and importing-on-demand + +Sometimes, one may wish to use a single context throughout a file. In this case, is possible to import all members +from a package or file, via `import context.complex.*`. Effectively, this is the same as enclosing an entire file +with a single context. However when using multiple contexts, this technique can introduce operator ambiguity, due to +namespace pollution. If there are multiple scoped contexts which define the same operation, it is still possible to +to import specific operations as needed, without using an explicit context with extension functions, for example: + +``` +import context.complex.op1 +import context.quaternion.op2 +``` diff --git a/doc/expressions.md b/doc/expressions.md new file mode 100644 index 000000000..1e05e5340 --- /dev/null +++ b/doc/expressions.md @@ -0,0 +1,26 @@ +# Expressions + +**Experimental: this API is in early stage and could change any time** + +Expressions is an experimental feature which allows to construct lazily or immediately calculated parametric mathematical +expressions. + +The potential use-cases for it (so far) are following: + +* Lazy evaluation (in general simple lambda is better, but there are some border cases) + +* Automatic differentiation in single-dimension and in multiple dimensions + +* Generation of mathematical syntax trees with subsequent code generation for other languages + +* Maybe symbolic computations (needs additional research) + +The workhorse of this API is `Expression` interface which exposes single `operator fun invoke(arguments: Map): T` +method. `ExpressionContext` is used to generate expressions and introduce variables. + +Currently there are two implementations: + +* Generic `ExpressionField` in `kmath-core` which allows construction of custom lazy expressions + +* Auto-differentiation expression in `kmath-commons` module allows to use full power of `DerivativeStructure` +from commons-math. **TODO: add example** diff --git a/doc/features.md b/doc/features.md new file mode 100644 index 000000000..e69de29bb diff --git a/doc/histograms.md b/doc/histograms.md new file mode 100644 index 000000000..6208deaff --- /dev/null +++ b/doc/histograms.md @@ -0,0 +1 @@ +**TODO** \ No newline at end of file diff --git a/doc/linear.md b/doc/linear.md new file mode 100644 index 000000000..6208deaff --- /dev/null +++ b/doc/linear.md @@ -0,0 +1 @@ +**TODO** \ No newline at end of file diff --git a/doc/nd-performance.md b/doc/nd-performance.md new file mode 100644 index 000000000..47653e3e8 --- /dev/null +++ b/doc/nd-performance.md @@ -0,0 +1,127 @@ +# Performance for n-dimensional structures operations + +One of the most sought after features of mathematical libraries is the high-performance operations on n-dimensional +structures. In `kmath` performance depends on which particular context was used for operation. + +Let us consider following contexts: +```kotlin + // specialized nd-field for Double. It works as generic Double field as well + val specializedField = NDField.real(intArrayOf(dim, dim)) + + // automatically build context most suited for given type. + val autoField = NDField.auto(intArrayOf(dim, dim), RealField) + + //A field implementing lazy computations. All elements are computed on-demand + val lazyField = NDField.lazy(intArrayOf(dim, dim), RealField) + + //A generic boxing field. It should be used for objects, not primitives. + val genericField = NDField.buffered(intArrayOf(dim, dim), RealField) +``` +Now let us perform several tests and see which implementation is best suited for each case: + +## Test case + +In order to test performance we will take 2d-structures with `dim = 1000` and add a structure filled with `1.0` +to it `n = 1000` times. + +## Specialized +The code to run this looks like: +```kotlin + specializedField.run { + var res = one + repeat(n) { + res += 1.0 + } + } +``` +The performance of this code is the best of all tests since it inlines all operations and is specialized for operation +with doubles. We will measure everything else relative to this one, so time for this test will be `1x` (real time +on my computer is about 4.5 seconds). The only problem with this approach is that it requires to specify type +from the beginning. Everyone do so anyway, so it is the recommended approach. + +## Automatic +Let's do the same with automatic field inference: +```kotlin + autoField.run { + var res = one + repeat(n) { + res += 1.0 + } + } +``` +Ths speed of this operation is approximately the same as for specialized case since `NDField.auto` just +returns the same `RealNDField` in this case. Of course it is usually better to use specialized method to be sure. + +## Lazy +Lazy field does not produce a structure when asked, instead it generates an empty structure and fills it on-demand +using coroutines to parallelize computations. +When one calls +```kotlin + lazyField.run { + var res = one + repeat(n) { + res += 1.0 + } + } +``` +The result will be calculated almost immediately but the result will be empty. In order to get the full result +structure one needs to call all its elements. In this case computation overhead will be huge. So this field never +should be used if one expects to use the full result structure. Though if one wants only small fraction, it could +save a lot of time. + +This field still could be used with reasonable performance if call code is changed: +```kotlin + lazyField.run { + val res = one.map { + var c = 0.0 + repeat(n) { + c += 1.0 + } + c + } + + res.elements().forEach { it.second } + } +``` +In this case it completes in about `4x-5x` time due to boxing. + +## Boxing +The boxing field produced by +```kotlin + genericField.run { + var res = one + repeat(n) { + res += 1.0 + } + } +``` +obviously is the slowest one, because it requires to box and unbox the `double` on each operation. It takes about +`15x` time (**TODO: there seems to be a problem here, it should be slow, but not that slow**). This field should +never be used for primitives. + +## Element operation +Let us also check the speed for direct operations on elements: +```kotlin + var res = genericField.one + repeat(n) { + res += 1.0 + } +``` +One would expect to be at least as slow as field operation, but in fact, this one takes only `2x` time to complete. +It happens, because in this particular case it does not use actual `NDField` but instead calculated directly +via extension function. + +## What about python? + +Usually it is bad idea to compare the direct numerical operation performance in different languages, but it hard to +work completely without frame of reference. In this case, simple numpy code: +```python +res = np.ones((1000,1000)) +for i in range(1000): + res = res + 1.0 +``` +gives the completion time of about `1.1x`, which means that specialized kotlin code in fact is working faster (I think it is +because better memory management). Of course if one writes `res += 1.0`, the performance will be different, +but it would be differenc case, because numpy overrides `+=` with in-place operations. In-place operations are +available in `kmath` with `MutableNDStructure` but there is no field for it (one can still work with mapping +functions). \ No newline at end of file diff --git a/doc/operations.md b/doc/operations.md new file mode 100644 index 000000000..425d791e8 --- /dev/null +++ b/doc/operations.md @@ -0,0 +1,40 @@ +## Spaces and fields + +An obvious first choice of mathematical objects to implement in a context-oriented style are algebraic elements like spaces, +rings and fields. Those are located in the `scientifik.kmath.operations.Algebra.kt` file. Alongside common contexts, the file includes definitions for algebra elements like `FieldElement`. A `FieldElement` object +stores a reference to the `Field` which contains additive and multiplicative operations, meaning +it has one fixed context attached and does not require explicit external context. So those `MathElements` can be operated without context: + +```kotlin +val c1 = Complex(1.0, 2.0) +val c2 = ComplexField.i +val c3 = c1 + c2 +``` + +`ComplexField` also features special operations to mix complex and real numbers, for example: + +```kotlin +val c1 = Complex(1.0, 2.0) +val c2 = ComplexField.run{ c1 - 1.0} // Returns: [re:0.0, im: 2.0] +val c3 = ComplexField.run{ c1 - i*2.0} +``` + +**Note**: In theory it is possible to add behaviors directly to the context, but currently kotlin syntax does not support +that. Watch [KT-10468](https://youtrack.jetbrains.com/issue/KT-10468) and [KEEP-176](https://github.com/Kotlin/KEEP/pull/176) for updates. + +## Nested fields + +Contexts allow one to build more complex structures. For example, it is possible to create a `Matrix` from complex elements like so: + +```kotlin +val element = NDElement.complex(shape = intArrayOf(2,2)){ index: IntArray -> + Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble()) +} +``` + +The `element` in this example is a member of the `Field` of 2-d structures, each element of which is a member of its own +`ComplexField`. The important thing is one does not need to create a special n-d class to hold complex +numbers and implement operations on it, one just needs to provide a field for its elements. + +**Note**: Fields themselves do not solve the problem of JVM boxing, but it is possible to solve with special contexts like +`BufferSpec`. This feature is in development phase. \ No newline at end of file diff --git a/gradle/wrapper/gradle-wrapper.jar b/gradle/wrapper/gradle-wrapper.jar new file mode 100644 index 000000000..457aad0d9 Binary files /dev/null and b/gradle/wrapper/gradle-wrapper.jar differ diff --git a/gradle/wrapper/gradle-wrapper.properties b/gradle/wrapper/gradle-wrapper.properties new file mode 100644 index 000000000..75b8c7c8c --- /dev/null +++ b/gradle/wrapper/gradle-wrapper.properties @@ -0,0 +1,5 @@ +distributionBase=GRADLE_USER_HOME +distributionPath=wrapper/dists +distributionUrl=https\://services.gradle.org/distributions/gradle-5.0-bin.zip +zipStoreBase=GRADLE_USER_HOME +zipStorePath=wrapper/dists diff --git a/gradlew b/gradlew new file mode 100755 index 000000000..af6708ff2 --- /dev/null +++ b/gradlew @@ -0,0 +1,172 @@ +#!/usr/bin/env sh + +############################################################################## +## +## Gradle start up script for UN*X +## +############################################################################## + +# Attempt to set APP_HOME +# Resolve links: $0 may be a link +PRG="$0" +# Need this for relative symlinks. +while [ -h "$PRG" ] ; do + ls=`ls -ld "$PRG"` + link=`expr "$ls" : '.*-> \(.*\)$'` + if expr "$link" : '/.*' > /dev/null; then + PRG="$link" + else + PRG=`dirname "$PRG"`"/$link" + fi +done +SAVED="`pwd`" +cd "`dirname \"$PRG\"`/" >/dev/null +APP_HOME="`pwd -P`" +cd "$SAVED" >/dev/null + +APP_NAME="Gradle" +APP_BASE_NAME=`basename "$0"` + +# Add default JVM options here. You can also use JAVA_OPTS and GRADLE_OPTS to pass JVM options to this script. +DEFAULT_JVM_OPTS='"-Xmx64m"' + +# Use the maximum available, or set MAX_FD != -1 to use that value. +MAX_FD="maximum" + +warn () { + echo "$*" +} + +die () { + echo + echo "$*" + echo + exit 1 +} + +# OS specific support (must be 'true' or 'false'). +cygwin=false +msys=false +darwin=false +nonstop=false +case "`uname`" in + CYGWIN* ) + cygwin=true + ;; + Darwin* ) + darwin=true + ;; + MINGW* ) + msys=true + ;; + NONSTOP* ) + nonstop=true + ;; +esac + +CLASSPATH=$APP_HOME/gradle/wrapper/gradle-wrapper.jar + +# Determine the Java command to use to start the JVM. +if [ -n "$JAVA_HOME" ] ; then + if [ -x "$JAVA_HOME/jre/sh/java" ] ; then + # IBM's JDK on AIX uses strange locations for the executables + JAVACMD="$JAVA_HOME/jre/sh/java" + else + JAVACMD="$JAVA_HOME/bin/java" + fi + if [ ! -x "$JAVACMD" ] ; then + die "ERROR: JAVA_HOME is set to an invalid directory: $JAVA_HOME + +Please set the JAVA_HOME variable in your environment to match the +location of your Java installation." + fi +else + JAVACMD="java" + which java >/dev/null 2>&1 || die "ERROR: JAVA_HOME is not set and no 'java' command could be found in your PATH. + +Please set the JAVA_HOME variable in your environment to match the +location of your Java installation." +fi + +# Increase the maximum file descriptors if we can. +if [ "$cygwin" = "false" -a "$darwin" = "false" -a "$nonstop" = "false" ] ; then + MAX_FD_LIMIT=`ulimit -H -n` + if [ $? -eq 0 ] ; then + if [ "$MAX_FD" = "maximum" -o "$MAX_FD" = "max" ] ; then + MAX_FD="$MAX_FD_LIMIT" + fi + ulimit -n $MAX_FD + if [ $? -ne 0 ] ; then + warn "Could not set maximum file descriptor limit: $MAX_FD" + fi + else + warn "Could not query maximum file descriptor limit: $MAX_FD_LIMIT" + fi +fi + +# For Darwin, add options to specify how the application appears in the dock +if $darwin; then + GRADLE_OPTS="$GRADLE_OPTS \"-Xdock:name=$APP_NAME\" \"-Xdock:icon=$APP_HOME/media/gradle.icns\"" +fi + +# For Cygwin, switch paths to Windows format before running java +if $cygwin ; then + APP_HOME=`cygpath --path --mixed "$APP_HOME"` + CLASSPATH=`cygpath --path --mixed "$CLASSPATH"` + JAVACMD=`cygpath --unix "$JAVACMD"` + + # We build the pattern for arguments to be converted via cygpath + ROOTDIRSRAW=`find -L / -maxdepth 1 -mindepth 1 -type d 2>/dev/null` + SEP="" + for dir in $ROOTDIRSRAW ; do + ROOTDIRS="$ROOTDIRS$SEP$dir" + SEP="|" + done + OURCYGPATTERN="(^($ROOTDIRS))" + # Add a user-defined pattern to the cygpath arguments + if [ "$GRADLE_CYGPATTERN" != "" ] ; then + OURCYGPATTERN="$OURCYGPATTERN|($GRADLE_CYGPATTERN)" + fi + # Now convert the arguments - kludge to limit ourselves to /bin/sh + i=0 + for arg in "$@" ; do + CHECK=`echo "$arg"|egrep -c "$OURCYGPATTERN" -` + CHECK2=`echo "$arg"|egrep -c "^-"` ### Determine if an option + + if [ $CHECK -ne 0 ] && [ $CHECK2 -eq 0 ] ; then ### Added a condition + eval `echo args$i`=`cygpath --path --ignore --mixed "$arg"` + else + eval `echo args$i`="\"$arg\"" + fi + i=$((i+1)) + done + case $i in + (0) set -- ;; + (1) set -- "$args0" ;; + (2) set -- "$args0" "$args1" ;; + (3) set -- "$args0" "$args1" "$args2" ;; + (4) set -- "$args0" "$args1" "$args2" "$args3" ;; + (5) set -- "$args0" "$args1" "$args2" "$args3" "$args4" ;; + (6) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" ;; + (7) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" ;; + (8) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" "$args7" ;; + (9) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" "$args7" "$args8" ;; + esac +fi + +# Escape application args +save () { + for i do printf %s\\n "$i" | sed "s/'/'\\\\''/g;1s/^/'/;\$s/\$/' \\\\/" ; done + echo " " +} +APP_ARGS=$(save "$@") + +# Collect all arguments for the java command, following the shell quoting and substitution rules +eval set -- $DEFAULT_JVM_OPTS $JAVA_OPTS $GRADLE_OPTS "\"-Dorg.gradle.appname=$APP_BASE_NAME\"" -classpath "\"$CLASSPATH\"" org.gradle.wrapper.GradleWrapperMain "$APP_ARGS" + +# by default we should be in the correct project dir, but when run from Finder on Mac, the cwd is wrong +if [ "$(uname)" = "Darwin" ] && [ "$HOME" = "$PWD" ]; then + cd "$(dirname "$0")" +fi + +exec "$JAVACMD" "$@" diff --git a/gradlew.bat b/gradlew.bat new file mode 100644 index 000000000..6d57edc70 --- /dev/null +++ b/gradlew.bat @@ -0,0 +1,84 @@ +@if "%DEBUG%" == "" @echo off +@rem ########################################################################## +@rem +@rem Gradle startup script for Windows +@rem +@rem ########################################################################## + +@rem Set local scope for the variables with windows NT shell +if "%OS%"=="Windows_NT" setlocal + +set DIRNAME=%~dp0 +if "%DIRNAME%" == "" set DIRNAME=. +set APP_BASE_NAME=%~n0 +set APP_HOME=%DIRNAME% + +@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" + +@rem Find java.exe +if defined JAVA_HOME goto findJavaFromJavaHome + +set JAVA_EXE=java.exe +%JAVA_EXE% -version >NUL 2>&1 +if "%ERRORLEVEL%" == "0" goto init + +echo. +echo ERROR: JAVA_HOME is not set and no 'java' command could be found in your PATH. +echo. +echo Please set the JAVA_HOME variable in your environment to match the +echo location of your Java installation. + +goto fail + +:findJavaFromJavaHome +set JAVA_HOME=%JAVA_HOME:"=% +set JAVA_EXE=%JAVA_HOME%/bin/java.exe + +if exist "%JAVA_EXE%" goto init + +echo. +echo ERROR: JAVA_HOME is set to an invalid directory: %JAVA_HOME% +echo. +echo Please set the JAVA_HOME variable in your environment to match the +echo location of your Java installation. + +goto fail + +:init +@rem Get command-line arguments, handling Windows variants + +if not "%OS%" == "Windows_NT" goto win9xME_args + +:win9xME_args +@rem Slurp the command line arguments. +set CMD_LINE_ARGS= +set _SKIP=2 + +:win9xME_args_slurp +if "x%~1" == "x" goto execute + +set CMD_LINE_ARGS=%* + +:execute +@rem Setup the command line + +set CLASSPATH=%APP_HOME%\gradle\wrapper\gradle-wrapper.jar + +@rem Execute Gradle +"%JAVA_EXE%" %DEFAULT_JVM_OPTS% %JAVA_OPTS% %GRADLE_OPTS% "-Dorg.gradle.appname=%APP_BASE_NAME%" -classpath "%CLASSPATH%" org.gradle.wrapper.GradleWrapperMain %CMD_LINE_ARGS% + +:end +@rem End local scope for the variables with windows NT shell +if "%ERRORLEVEL%"=="0" goto mainEnd + +:fail +rem Set variable GRADLE_EXIT_CONSOLE if you need the _script_ return code instead of +rem the _cmd.exe /c_ return code! +if not "" == "%GRADLE_EXIT_CONSOLE%" exit 1 +exit /b 1 + +:mainEnd +if "%OS%"=="Windows_NT" endlocal + +:omega diff --git a/kmath-commons/build.gradle.kts b/kmath-commons/build.gradle.kts new file mode 100644 index 000000000..305761c01 --- /dev/null +++ b/kmath-commons/build.gradle.kts @@ -0,0 +1,13 @@ +plugins { + kotlin("jvm") +} + +description = "Commons math binding for kmath" + +dependencies { + api(project(":kmath-core")) + api(project(":kmath-sequential")) + api("org.apache.commons:commons-math3:3.6.1") + testImplementation("org.jetbrains.kotlin:kotlin-test") + testImplementation("org.jetbrains.kotlin:kotlin-test-junit") +} diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/expressions/DiffExpression.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/expressions/DiffExpression.kt new file mode 100644 index 000000000..6c1759d54 --- /dev/null +++ b/kmath-commons/src/main/kotlin/scientifik/kmath/expressions/DiffExpression.kt @@ -0,0 +1,129 @@ +package scientifik.kmath.expressions + +import org.apache.commons.math3.analysis.differentiation.DerivativeStructure +import scientifik.kmath.operations.ExtendedField +import scientifik.kmath.operations.ExtendedFieldOperations +import scientifik.kmath.operations.Field +import kotlin.properties.ReadOnlyProperty +import kotlin.reflect.KProperty + +/** + * A field wrapping commons-math derivative structures + */ +class DerivativeStructureField(val order: Int, val parameters: Map) : + ExtendedField { + + override val zero: DerivativeStructure by lazy { DerivativeStructure(order, parameters.size) } + + override val one: DerivativeStructure by lazy { DerivativeStructure(order, parameters.size, 1.0) } + + private val variables: Map = parameters.mapValues { (key, value) -> + DerivativeStructure(parameters.size, order, parameters.keys.indexOf(key), value) + } + + val variable = object : ReadOnlyProperty { + override fun getValue(thisRef: Any?, property: KProperty<*>): DerivativeStructure { + return variables[property.name] ?: error("A variable with name ${property.name} does not exist") + } + } + + fun variable(name: String, default: DerivativeStructure? = null): DerivativeStructure = + variables[name] ?: default ?: error("A variable with name $name does not exist") + + + fun Number.const() = DerivativeStructure(order, parameters.size, toDouble()) + + fun DerivativeStructure.deriv(parName: String, order: Int = 1): Double { + return deriv(mapOf(parName to order)) + } + + fun DerivativeStructure.deriv(orders: Map): Double { + return getPartialDerivative(*parameters.keys.map { orders[it] ?: 0 }.toIntArray()) + } + + fun DerivativeStructure.deriv(vararg orders: Pair): Double = deriv(mapOf(*orders)) + + override fun add(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.add(b) + + override fun multiply(a: DerivativeStructure, k: Number): DerivativeStructure = when (k) { + is Double -> a.multiply(k) + is Int -> a.multiply(k) + else -> a.multiply(k.toDouble()) + } + + override fun multiply(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.multiply(b) + + override fun divide(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.divide(b) + + override fun sin(arg: DerivativeStructure): DerivativeStructure = arg.sin() + + override fun cos(arg: DerivativeStructure): DerivativeStructure = arg.cos() + + override fun power(arg: DerivativeStructure, pow: Number): DerivativeStructure = when (pow) { + is Double -> arg.pow(pow) + is Int -> arg.pow(pow) + else -> arg.pow(pow.toDouble()) + } + + fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure = arg.pow(pow) + + override fun exp(arg: DerivativeStructure): DerivativeStructure = arg.exp() + + override fun ln(arg: DerivativeStructure): DerivativeStructure = arg.log() + + operator fun DerivativeStructure.plus(n: Number): DerivativeStructure = add(n.toDouble()) + operator fun DerivativeStructure.minus(n: Number): DerivativeStructure = subtract(n.toDouble()) + operator fun Number.plus(s: DerivativeStructure) = s + this + operator fun Number.minus(s: DerivativeStructure) = s - this +} + +/** + * A constructs that creates a derivative structure with required order on-demand + */ +class DiffExpression(val function: DerivativeStructureField.() -> DerivativeStructure) : Expression { + override fun invoke(arguments: Map): Double = DerivativeStructureField(0, arguments) + .run(function).value + + /** + * Get the derivative expression with given orders + * TODO make result [DiffExpression] + */ + fun derivative(orders: Map): Expression { + return object : Expression { + override fun invoke(arguments: Map): Double = + DerivativeStructureField(orders.values.max() ?: 0, arguments) + .run { + function().deriv(orders) + } + } + } + + //TODO add gradient and maybe other vector operators +} + +fun DiffExpression.derivative(vararg orders: Pair) = derivative(mapOf(*orders)) +fun DiffExpression.derivative(name: String) = derivative(name to 1) + +/** + * A context for [DiffExpression] (not to be confused with [DerivativeStructure]) + */ +object DiffExpressionContext : ExpressionContext, Field { + override fun variable(name: String, default: Double?) = DiffExpression { variable(name, default?.const()) } + + override fun const(value: Double): DiffExpression = DiffExpression { value.const() } + + override fun add(a: DiffExpression, b: DiffExpression) = DiffExpression { a.function(this) + b.function(this) } + + override val zero = DiffExpression { 0.0.const() } + + override fun multiply(a: DiffExpression, k: Number) = DiffExpression { a.function(this) * k } + + override val one = DiffExpression { 1.0.const() } + + override fun multiply(a: DiffExpression, b: DiffExpression) = DiffExpression { a.function(this) * b.function(this) } + + override fun divide(a: DiffExpression, b: DiffExpression) = DiffExpression { a.function(this) / b.function(this) } +} + + + diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt new file mode 100644 index 000000000..808b2768b --- /dev/null +++ b/kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt @@ -0,0 +1,93 @@ +package scientifik.kmath.linear + +import org.apache.commons.math3.linear.* +import org.apache.commons.math3.linear.RealMatrix +import org.apache.commons.math3.linear.RealVector + +class CMMatrix(val origin: RealMatrix, features: Set? = null) : Matrix { + override val rowNum: Int get() = origin.rowDimension + override val colNum: Int get() = origin.columnDimension + + override val features: Set = features ?: sequence { + if(origin is DiagonalMatrix) yield(DiagonalFeature) + }.toSet() + + override fun suggestFeature(vararg features: MatrixFeature) = + CMMatrix(origin, this.features + features) + + override fun get(i: Int, j: Int): Double = origin.getEntry(i, j) +} + +fun Matrix.toCM(): CMMatrix = if (this is CMMatrix) { + this +} else { + //TODO add feature analysis + val array = Array(rowNum) { i -> DoubleArray(colNum) { j -> get(i, j) } } + CMMatrix(Array2DRowRealMatrix(array)) +} + +fun RealMatrix.toMatrix() = CMMatrix(this) + +class CMVector(val origin: RealVector) : Point { + override val size: Int get() = origin.dimension + + override fun get(index: Int): Double = origin.getEntry(index) + + override fun iterator(): Iterator = origin.toArray().iterator() +} + +fun Point.toCM(): CMVector = if (this is CMVector) { + this +} else { + val array = DoubleArray(size) { this[it] } + CMVector(ArrayRealVector(array)) +} + +fun RealVector.toPoint() = CMVector(this) + +object CMMatrixContext : MatrixContext, LinearSolver { + + override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix { + val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } } + return CMMatrix(Array2DRowRealMatrix(array)) + } + + override fun solve(a: Matrix, b: Matrix): CMMatrix { + val decomposition = LUDecomposition(a.toCM().origin) + return decomposition.solver.solve(b.toCM().origin).toMatrix() + } + + override fun solve(a: Matrix, b: Point): CMVector { + val decomposition = LUDecomposition(a.toCM().origin) + return decomposition.solver.solve(b.toCM().origin).toPoint() + } + + override fun inverse(a: Matrix): CMMatrix { + val decomposition = LUDecomposition(a.toCM().origin) + return decomposition.solver.inverse.toMatrix() + } + + override fun Matrix.dot(other: Matrix) = + CMMatrix(this.toCM().origin.multiply(other.toCM().origin)) + + override fun Matrix.dot(vector: Point): CMVector = + CMVector(this.toCM().origin.preMultiply(vector.toCM().origin)) + + + override fun Matrix.unaryMinus(): CMMatrix = + produce(rowNum, colNum) { i, j -> -get(i, j) } + + override fun Matrix.plus(b: Matrix) = + CMMatrix(this.toCM().origin.multiply(b.toCM().origin)) + + override fun Matrix.minus(b: Matrix) = + CMMatrix(this.toCM().origin.subtract(b.toCM().origin)) + + override fun Matrix.times(value: Double) = + CMMatrix(this.toCM().origin.scalarMultiply(value.toDouble())) +} + +operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.add(other.origin)) +operator fun CMMatrix.minus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.subtract(other.origin)) + +infix fun CMMatrix.dot(other: CMMatrix): CMMatrix = CMMatrix(this.origin.multiply(other.origin)) \ No newline at end of file diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/transform/Transformations.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/transform/Transformations.kt new file mode 100644 index 000000000..17907adbe --- /dev/null +++ b/kmath-commons/src/main/kotlin/scientifik/kmath/transform/Transformations.kt @@ -0,0 +1,86 @@ +package scientifik.kmath.transform + +import org.apache.commons.math3.transform.* +import scientifik.kmath.operations.Complex +import scientifik.kmath.sequential.Processor +import scientifik.kmath.sequential.Producer +import scientifik.kmath.sequential.map +import scientifik.kmath.structures.* + + +/** + * + */ +object Transformations { + + private fun Buffer.toArray(): Array = + Array(size) { org.apache.commons.math3.complex.Complex(get(it).re, get(it).im) } + + private fun Buffer.asArray() = if (this is DoubleBuffer) { + array + } else { + DoubleArray(size) { i -> get(i) } + } + + /** + * Create a virtual buffer on top of array + */ + private fun Array.asBuffer() = VirtualBuffer(size) { + val value = get(it) + Complex(value.real, value.imaginary) + } + + fun fourier( + normalization: DftNormalization = DftNormalization.STANDARD, + direction: TransformType = TransformType.FORWARD + ): BufferTransform = { + FastFourierTransformer(normalization).transform(it.toArray(), direction).asBuffer() + } + + fun realFourier( + normalization: DftNormalization = DftNormalization.STANDARD, + direction: TransformType = TransformType.FORWARD + ): BufferTransform = { + FastFourierTransformer(normalization).transform(it.asArray(), direction).asBuffer() + } + + fun sine( + normalization: DstNormalization = DstNormalization.STANDARD_DST_I, + direction: TransformType = TransformType.FORWARD + ): BufferTransform = { + FastSineTransformer(normalization).transform(it.asArray(), direction).asBuffer() + } + + fun cosine( + normalization: DctNormalization = DctNormalization.STANDARD_DCT_I, + direction: TransformType = TransformType.FORWARD + ): BufferTransform = { + FastCosineTransformer(normalization).transform(it.asArray(), direction).asBuffer() + } + + fun hadamard( + direction: TransformType = TransformType.FORWARD + ): BufferTransform = { + FastHadamardTransformer().transform(it.asArray(), direction).asBuffer() + } +} + +/** + * Process given [Producer] with commons-math fft transformation + */ +fun Producer>.FFT( + normalization: DftNormalization = DftNormalization.STANDARD, + direction: TransformType = TransformType.FORWARD +): Processor, Buffer> { + val transform = Transformations.fourier(normalization, direction) + return map { transform(it) } +} + +@JvmName("realFFT") +fun Producer>.FFT( + normalization: DftNormalization = DftNormalization.STANDARD, + direction: TransformType = TransformType.FORWARD +): Processor, Buffer> { + val transform = Transformations.realFourier(normalization, direction) + return map { transform(it) } +} \ No newline at end of file diff --git a/kmath-commons/src/test/kotlin/scientifik/kmath/expressions/AutoDiffTest.kt b/kmath-commons/src/test/kotlin/scientifik/kmath/expressions/AutoDiffTest.kt new file mode 100644 index 000000000..2695bb918 --- /dev/null +++ b/kmath-commons/src/test/kotlin/scientifik/kmath/expressions/AutoDiffTest.kt @@ -0,0 +1,31 @@ +package scientifik.kmath.expressions + +import org.junit.Test +import kotlin.test.assertEquals + +inline fun diff(order: Int, vararg parameters: Pair, block: DerivativeStructureField.() -> R) = + DerivativeStructureField(order, mapOf(*parameters)).run(block) + +class AutoDiffTest { + @Test + fun derivativeStructureFieldTest() { + val res = diff(3, "x" to 1.0, "y" to 1.0) { + val x by variable + val y = variable("y") + val z = x * (-sin(x * y) + y) + z.deriv("x") + } + } + + @Test + fun autoDifTest() { + val f = DiffExpression { + val x by variable + val y by variable + x.pow(2) + 2 * x * y + y.pow(2) + 1 + } + + assertEquals(10.0, f("x" to 1.0, "y" to 2.0)) + assertEquals(6.0, f.derivative("x")("x" to 1.0, "y" to 2.0)) + } +} \ No newline at end of file diff --git a/kmath-core/build.gradle b/kmath-core/build.gradle deleted file mode 100644 index d51ccdb65..000000000 --- a/kmath-core/build.gradle +++ /dev/null @@ -1,52 +0,0 @@ -plugins { - id "org.jetbrains.kotlin.multiplatform" -} - -kotlin { - targets { - fromPreset(presets.jvm, 'jvm') - fromPreset(presets.js, 'js') - // For ARM, preset should be changed to presets.iosArm32 or presets.iosArm64 - // For Linux, preset should be changed to e.g. presets.linuxX64 - // For MacOS, preset should be changed to e.g. presets.macosX64 - //fromPreset(presets.mingwX64, 'mingw') - } - sourceSets { - commonMain { - dependencies { - api 'org.jetbrains.kotlin:kotlin-stdlib-common' - } - } - commonTest { - dependencies { - implementation 'org.jetbrains.kotlin:kotlin-test-common' - implementation 'org.jetbrains.kotlin:kotlin-test-annotations-common' - } - } - jvmMain { - dependencies { - api 'org.jetbrains.kotlin:kotlin-stdlib-jdk8' - } - } - jvmTest { - dependencies { - implementation 'org.jetbrains.kotlin:kotlin-test' - implementation 'org.jetbrains.kotlin:kotlin-test-junit' - } - } - jsMain { - dependencies { - api 'org.jetbrains.kotlin:kotlin-stdlib-js' - } - } - jsTest { - dependencies { - implementation 'org.jetbrains.kotlin:kotlin-test-js' - } - } -// mingwMain { -// } -// mingwTest { -// } - } -} diff --git a/kmath-core/build.gradle.kts b/kmath-core/build.gradle.kts new file mode 100644 index 000000000..bd125b373 --- /dev/null +++ b/kmath-core/build.gradle.kts @@ -0,0 +1,51 @@ +plugins { + kotlin("multiplatform") +} + +val ioVersion: String by rootProject.extra + + +kotlin { + jvm() + js() + + sourceSets { + val commonMain by getting { + dependencies { + api(project(":kmath-memory")) + api(kotlin("stdlib")) + } + } + val commonTest by getting { + dependencies { + implementation(kotlin("test-common")) + implementation(kotlin("test-annotations-common")) + } + } + val jvmMain by getting { + dependencies { + api(kotlin("stdlib-jdk8")) + } + } + val jvmTest by getting { + dependencies { + implementation(kotlin("test")) + implementation(kotlin("test-junit")) + } + } + val jsMain by getting { + dependencies { + api(kotlin("stdlib-js")) + } + } + val jsTest by getting { + dependencies { + implementation(kotlin("test-js")) + } + } +// mingwMain { +// } +// mingwTest { +// } + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt index 0a34b536c..ad345bc5a 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt @@ -1,48 +1,64 @@ package scientifik.kmath.expressions import scientifik.kmath.operations.Field +import scientifik.kmath.operations.Ring import scientifik.kmath.operations.Space - +/** + * An elementary function that could be invoked on a map of arguments + */ interface Expression { operator fun invoke(arguments: Map): T } operator fun Expression.invoke(vararg pairs: Pair): T = invoke(mapOf(*pairs)) +/** + * A context for expression construction + */ interface ExpressionContext { + /** + * Introduce a variable into expression context + */ fun variable(name: String, default: T? = null): Expression + /** + * A constant expression which does not depend on arguments + */ fun const(value: T): Expression } internal class VariableExpression(val name: String, val default: T? = null) : Expression { - override fun invoke(arguments: Map): T { - return arguments[name] ?: default ?: error("The parameter not found: $name") - } + override fun invoke(arguments: Map): T = + arguments[name] ?: default ?: error("Parameter not found: $name") } internal class ConstantExpression(val value: T) : Expression { override fun invoke(arguments: Map): T = value } -internal class SumExpression(val context: Space, val first: Expression, val second: Expression) : Expression { +internal class SumExpression(val context: Space, val first: Expression, val second: Expression) : + Expression { override fun invoke(arguments: Map): T = context.add(first.invoke(arguments), second.invoke(arguments)) } -internal class ProductExpression(val context: Field, val first: Expression, val second: Expression) : Expression { - override fun invoke(arguments: Map): T = context.multiply(first.invoke(arguments), second.invoke(arguments)) +internal class ProductExpression(val context: Ring, val first: Expression, val second: Expression) : + Expression { + override fun invoke(arguments: Map): T = + context.multiply(first.invoke(arguments), second.invoke(arguments)) } -internal class ConstProductExpession(val context: Field, val expr: Expression, val const: Double) : Expression { +internal class ConstProductExpession(val context: Space, val expr: Expression, val const: Number) : + Expression { override fun invoke(arguments: Map): T = context.multiply(expr.invoke(arguments), const) } -internal class DivExpession(val context: Field, val expr: Expression, val second: Expression) : Expression { +internal class DivExpession(val context: Field, val expr: Expression, val second: Expression) : + Expression { override fun invoke(arguments: Map): T = context.divide(expr.invoke(arguments), second.invoke(arguments)) } -class FieldExpressionContext(val field: Field) : Field>, ExpressionContext { +class ExpressionField(val field: Field) : Field>, ExpressionContext { override val zero: Expression = ConstantExpression(field.zero) @@ -54,9 +70,20 @@ class FieldExpressionContext(val field: Field) : Field>, Exp override fun add(a: Expression, b: Expression): Expression = SumExpression(field, a, b) - override fun multiply(a: Expression, k: Double): Expression = ConstProductExpession(field, a, k) + override fun multiply(a: Expression, k: Number): Expression = ConstProductExpession(field, a, k) override fun multiply(a: Expression, b: Expression): Expression = ProductExpression(field, a, b) override fun divide(a: Expression, b: Expression): Expression = DivExpession(field, a, b) + + + operator fun Expression.plus(arg: T) = this + const(arg) + operator fun Expression.minus(arg: T) = this - const(arg) + operator fun Expression.times(arg: T) = this * const(arg) + operator fun Expression.div(arg: T) = this / const(arg) + + operator fun T.plus(arg: Expression) = arg + this + operator fun T.minus(arg: Expression) = arg - this + operator fun T.times(arg: Expression) = arg * this + operator fun T.div(arg: Expression) = arg / this } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Counters.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Counters.kt deleted file mode 100644 index f6cc1f822..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Counters.kt +++ /dev/null @@ -1,20 +0,0 @@ -package scientifik.kmath.histogram - -/* - * Common representation for atomic counters - */ - - -expect class LongCounter(){ - fun decrement() - fun increment() - fun reset() - fun sum(): Long - fun add(l:Long) -} - -expect class DoubleCounter(){ - fun reset() - fun sum(): Double - fun add(d: Double) -} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt deleted file mode 100644 index 08214142e..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt +++ /dev/null @@ -1,63 +0,0 @@ -package scientifik.kmath.histogram - -import scientifik.kmath.structures.ArrayBuffer -import scientifik.kmath.structures.Buffer -import scientifik.kmath.structures.DoubleBuffer - -typealias Point = Buffer - -typealias RealPoint = Buffer - -/** - * A simple geometric domain - * TODO move to geometry module - */ -interface Domain { - operator fun contains(vector: Point): Boolean - val dimension: Int -} - -/** - * The bin in the histogram. The histogram is by definition always done in the real space - */ -interface Bin : Domain { - /** - * The value of this bin - */ - val value: Number - val center: Point -} - -interface Histogram> : Iterable { - - /** - * Find existing bin, corresponding to given coordinates - */ - operator fun get(point: Point): B? - - /** - * Dimension of the histogram - */ - val dimension: Int - -} - -interface MutableHistogram>: Histogram{ - - /** - * Increment appropriate bin - */ - fun put(point: Point, weight: Double = 1.0) -} - -fun MutableHistogram.put(vararg point: T) = put(ArrayBuffer(point)) - -fun MutableHistogram.put(vararg point: Number) = put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray())) -fun MutableHistogram.put(vararg point: Double) = put(DoubleBuffer(point)) - -fun MutableHistogram.fill(sequence: Iterable>) = sequence.forEach { put(it) } - -/** - * Pass a sequence builder into histogram - */ -fun MutableHistogram.fill(buider: suspend SequenceScope>.() -> Unit) = fill(sequence(buider).asIterable()) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt deleted file mode 100644 index ffffb0d7d..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt +++ /dev/null @@ -1,63 +0,0 @@ -package scientifik.kmath.histogram - -import scientifik.kmath.linear.Vector -import scientifik.kmath.operations.Space -import scientifik.kmath.structures.NDStructure -import scientifik.kmath.structures.asSequence - -data class BinTemplate>(val center: Vector, val sizes: Point) { - fun contains(vector: Point): Boolean { - if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}") - val upper = center.context.run { center + sizes / 2.0} - val lower = center.context.run {center - sizes / 2.0} - return vector.asSequence().mapIndexed { i, value -> - value in lower[i]..upper[i] - }.all { it } - } -} - -/** - * A space to perform arithmetic operations on histograms - */ -interface HistogramSpace, H : Histogram> : Space { - /** - * Rules for performing operations on bins - */ - val binSpace: Space> -} - -class PhantomBin>(val template: BinTemplate, override val value: Number) : Bin { - - override fun contains(vector: Point): Boolean = template.contains(vector) - - override val dimension: Int - get() = template.center.size - - override val center: Point - get() = template.center - -} - -/** - * Immutable histogram with explicit structure for content and additional external bin description. - * Bin search is slow, but full histogram algebra is supported. - * @param bins map a template into structure index - */ -class PhantomHistogram>( - val bins: Map, IntArray>, - val data: NDStructure -) : Histogram> { - - override val dimension: Int - get() = data.dimension - - override fun iterator(): Iterator> { - return bins.asSequence().map { entry -> PhantomBin(entry.key, data[entry.value]) }.iterator() - } - - override fun get(point: Point): PhantomBin? { - val template = bins.keys.find { it.contains(point) } - return template?.let { PhantomBin(it, data[bins[it]!!]) } - } - -} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt new file mode 100644 index 000000000..3752f6db5 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt @@ -0,0 +1,100 @@ +package scientifik.kmath.linear + +import scientifik.kmath.operations.Ring +import scientifik.kmath.structures.* + +/** + * Basic implementation of Matrix space based on [NDStructure] + */ +class BufferMatrixContext>( + override val elementContext: R, + private val bufferFactory: BufferFactory +) : GenericMatrixContext { + + override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix { + val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) } + return BufferMatrix(rows, columns, buffer) + } + + override fun point(size: Int, initializer: (Int) -> T): Point = bufferFactory(size, initializer) +} + +class BufferMatrix( + override val rowNum: Int, + override val colNum: Int, + val buffer: Buffer, + override val features: Set = emptySet() +) : Matrix { + + init { + if (buffer.size != rowNum * colNum) { + error("Dimension mismatch for matrix structure") + } + } + + + override val shape: IntArray get() = intArrayOf(rowNum, colNum) + + override fun suggestFeature(vararg features: MatrixFeature) = + BufferMatrix(rowNum, colNum, buffer, this.features + features) + + override fun get(index: IntArray): T = get(index[0], index[1]) + + override fun get(i: Int, j: Int): T = buffer[i * colNum + j] + + override fun elements(): Sequence> = sequence { + for (i in 0 until rowNum) { + for (j in 0 until colNum) { + yield(intArrayOf(i, j) to get(i, j)) + } + } + } + + override fun equals(other: Any?): Boolean { + if (this === other) return true + return when (other) { + is NDStructure<*> -> return NDStructure.equals(this, other) + else -> false + } + } + + override fun hashCode(): Int { + var result = buffer.hashCode() + result = 31 * result + features.hashCode() + return result + } + + override fun toString(): String { + return if (rowNum <= 5 && colNum <= 5) { + "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)\n" + + rows.asSequence().joinToString(prefix = "(", postfix = ")", separator = "\n ") { + it.asSequence().joinToString(separator = "\t") { it.toString() } + } + } else { + "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)" + } + } +} + +/** + * Optimized dot product for real matrices + */ +infix fun BufferMatrix.dot(other: BufferMatrix): BufferMatrix { + if (this.colNum != other.rowNum) error("Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})") + + val array = DoubleArray(this.rowNum * other.colNum) + + val a = this.buffer.array + val b = other.buffer.array + + for (i in (0 until rowNum)) { + for (j in (0 until other.colNum)) { + for (k in (0 until colNum)) { + array[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j] + } + } + } + + val buffer = DoubleBuffer(array) + return BufferMatrix(rowNum, other.colNum, buffer) +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt deleted file mode 100644 index 223e63419..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt +++ /dev/null @@ -1,251 +0,0 @@ -package scientifik.kmath.linear - -import scientifik.kmath.operations.DoubleField -import scientifik.kmath.operations.Field -import scientifik.kmath.structures.MutableNDStructure -import scientifik.kmath.structures.NDStructure -import scientifik.kmath.structures.genericNdStructure -import scientifik.kmath.structures.get -import kotlin.math.absoluteValue - -/** - * Implementation based on Apache common-maths LU-decomposition - */ -abstract class LUDecomposition, F : Field>(val matrix: Matrix) { - - private val field get() = matrix.context.field - /** Entries of LU decomposition. */ - internal val lu: NDStructure - /** Pivot permutation associated with LU decomposition. */ - internal val pivot: IntArray - /** Parity of the permutation associated with the LU decomposition. */ - private var even: Boolean = false - - init { - val pair = calculateLU() - lu = pair.first - pivot = pair.second - } - - /** - * Returns the matrix L of the decomposition. - * - * L is a lower-triangular matrix - * @return the L matrix (or null if decomposed matrix is singular) - */ - val l: Matrix by lazy { - matrix.context.produce { i, j -> - when { - j < i -> lu[i, j] - j == i -> matrix.context.field.one - else -> matrix.context.field.zero - } - } - } - - - /** - * Returns the matrix U of the decomposition. - * - * U is an upper-triangular matrix - * @return the U matrix (or null if decomposed matrix is singular) - */ - val u: Matrix by lazy { - matrix.context.produce { i, j -> - if (j >= i) lu[i, j] else field.zero - } - } - - /** - * Returns the P rows permutation matrix. - * - * P is a sparse matrix with exactly one element set to 1.0 in - * each row and each column, all other elements being set to 0.0. - * - * The positions of the 1 elements are given by the [ pivot permutation vector][.getPivot]. - * @return the P rows permutation matrix (or null if decomposed matrix is singular) - * @see .getPivot - */ - val p: Matrix by lazy { - matrix.context.produce { i, j -> - //TODO ineffective. Need sparse matrix for that - if (j == pivot[i]) field.one else field.zero - } - } - - /** - * Return the determinant of the matrix - * @return determinant of the matrix - */ - val determinant: T - get() { - with(matrix.context.field) { - var determinant = if (even) one else -one - for (i in 0 until matrix.rows) { - determinant *= lu[i, i] - } - return determinant - } - } - - /** - * In-place transformation for [MutableNDArray], using given transformation for each element - */ - operator fun MutableNDStructure.set(i: Int, j: Int, value: T) { - this[intArrayOf(i, j)] = value - } - - abstract fun isSingular(value: T): Boolean - - private fun abs(value: T) = if (value > matrix.context.field.zero) value else with(matrix.context.field) { -value } - - private fun calculateLU(): Pair, IntArray> { - if (matrix.rows != matrix.columns) { - error("LU decomposition supports only square matrices") - } - - val m = matrix.columns - val pivot = IntArray(matrix.rows) - //TODO fix performance - val lu: MutableNDStructure = genericNdStructure(intArrayOf(matrix.rows, matrix.columns)) { index -> matrix[index[0], index[1]] } - - - with(matrix.context.field) { - // Initialize permutation array and parity - for (row in 0 until m) { - pivot[row] = row - } - even = true - - // Loop over columns - for (col in 0 until m) { - - // upper - for (row in 0 until col) { - var sum = lu[row, col] - for (i in 0 until row) { - sum -= lu[row, i] * lu[i, col] - } - lu[row, col] = sum - } - - // lower - val max = (col until m).maxBy { row -> - var sum = lu[row, col] - for (i in 0 until col) { - sum -= lu[row, i] * lu[i, col] - } - //luRow[col] = sum - lu[row, col] = sum - - abs(sum) - } ?: col - - // Singularity check - if (isSingular(lu[max, col])) { - error("Singular matrix") - } - - // Pivot if necessary - if (max != col) { - //var tmp = zero - //val luMax = lu[max] - //val luCol = lu[col] - for (i in 0 until m) { - lu[max, i] = lu[col, i] - lu[col, i] = lu[max, i] - } - val temp = pivot[max] - pivot[max] = pivot[col] - pivot[col] = temp - even = !even - } - - // Divide the lower elements by the "winning" diagonal elt. - val luDiag = lu[col, col] - for (row in col + 1 until m) { - lu[row, col] = lu[row, col] / luDiag -// lu[row, col] /= luDiag - } - } - } - return Pair(lu, pivot) - } - - /** - * Returns the pivot permutation vector. - * @return the pivot permutation vector - * @see .getP - */ - fun getPivot(): IntArray { - return pivot.copyOf() - } - -} - -class RealLUDecomposition(matrix: RealMatrix, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition(matrix) { - override fun isSingular(value: Double): Boolean { - return value.absoluteValue < singularityThreshold - } - - companion object { - /** Default bound to determine effective singularity in LU decomposition. */ - private const val DEFAULT_TOO_SMALL = 1e-11 - } -} - - -/** Specialized solver. */ -object RealLUSolver : LinearSolver { - - - fun decompose(mat: Matrix, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold) - - override fun solve(a: RealMatrix, b: RealMatrix): RealMatrix { - val decomposition = decompose(a, a.context.field.zero) - - if (b.rows != a.rows) { - error("Matrix dimension mismatch expected ${a.rows}, but got ${b.rows}") - } - - // Apply permutations to b - val bp = Array(a.rows) { DoubleArray(b.columns) } - for (row in 0 until a.rows) { - val bpRow = bp[row] - val pRow = decomposition.pivot[row] - for (col in 0 until b.columns) { - bpRow[col] = b[pRow, col] - } - } - - // Solve LY = b - for (col in 0 until a.rows) { - val bpCol = bp[col] - for (i in col + 1 until a.rows) { - val bpI = bp[i] - val luICol = decomposition.lu[i, col] - for (j in 0 until b.columns) { - bpI[j] -= bpCol[j] * luICol - } - } - } - - // Solve UX = Y - for (col in a.rows - 1 downTo 0) { - val bpCol = bp[col] - val luDiag = decomposition.lu[col, col] - for (j in 0 until b.columns) { - bpCol[j] /= luDiag - } - for (i in 0 until col) { - val bpI = bp[i] - val luICol = decomposition.lu[i, col] - for (j in 0 until b.columns) { - bpI[j] -= bpCol[j] * luICol - } - } - } - - return a.context.produce { i, j -> bp[i][j] } - } -} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt new file mode 100644 index 000000000..85ddb8786 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt @@ -0,0 +1,205 @@ +package scientifik.kmath.linear + +import scientifik.kmath.operations.Field +import scientifik.kmath.operations.Ring +import scientifik.kmath.structures.MutableBuffer +import scientifik.kmath.structures.MutableBuffer.Companion.boxing +import scientifik.kmath.structures.MutableBufferFactory +import scientifik.kmath.structures.NDStructure +import scientifik.kmath.structures.get + +class LUPDecomposition>( + private val elementContext: Ring, + internal val lu: NDStructure, + val pivot: IntArray, + private val even: Boolean +) : LUPDecompositionFeature, DeterminantFeature { + + /** + * Returns the matrix L of the decomposition. + * + * L is a lower-triangular matrix with [Ring.one] in diagonal + */ + override val l: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1], setOf(LFeature)) { i, j -> + when { + j < i -> lu[i, j] + j == i -> elementContext.one + else -> elementContext.zero + } + } + + + /** + * Returns the matrix U of the decomposition. + * + * U is an upper-triangular matrix including the diagonal + */ + override val u: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1], setOf(UFeature)) { i, j -> + if (j >= i) lu[i, j] else elementContext.zero + } + + + /** + * Returns the P rows permutation matrix. + * + * P is a sparse matrix with exactly one element set to [Ring.one] in + * each row and each column, all other elements being set to [Ring.zero]. + */ + override val p: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> + if (j == pivot[i]) elementContext.one else elementContext.zero + } + + + /** + * Return the determinant of the matrix + * @return determinant of the matrix + */ + override val determinant: T by lazy { + with(elementContext) { + (0 until lu.shape[0]).fold(if (even) one else -one) { value, i -> value * lu[i, i] } + } + } + +} + + +class LUSolver, F : Field>( + val context: GenericMatrixContext, + val bufferFactory: MutableBufferFactory = ::boxing, + val singularityCheck: (T) -> Boolean +) : LinearSolver { + + + private fun abs(value: T) = + if (value > context.elementContext.zero) value else with(context.elementContext) { -value } + + fun buildDecomposition(matrix: Matrix): LUPDecomposition { + if (matrix.rowNum != matrix.colNum) { + error("LU decomposition supports only square matrices") + } + + val m = matrix.colNum + val pivot = IntArray(matrix.rowNum) + + val lu = Mutable2DStructure.create(matrix.rowNum, matrix.colNum, bufferFactory) { i, j -> + matrix[i, j] + } + + + with(context.elementContext) { + // Initialize permutation array and parity + for (row in 0 until m) { + pivot[row] = row + } + var even = true + + // Loop over columns + for (col in 0 until m) { + + // upper + for (row in 0 until col) { + var sum = lu[row, col] + for (i in 0 until row) { + sum -= lu[row, i] * lu[i, col] + } + lu[row, col] = sum + } + + // lower + val max = (col until m).maxBy { row -> + var sum = lu[row, col] + for (i in 0 until col) { + sum -= lu[row, i] * lu[i, col] + } + lu[row, col] = sum + + abs(sum) + } ?: col + + // Singularity check + if (singularityCheck(lu[max, col])) { + error("Singular matrix") + } + + // Pivot if necessary + if (max != col) { + for (i in 0 until m) { + lu[max, i] = lu[col, i] + lu[col, i] = lu[max, i] + } + val temp = pivot[max] + pivot[max] = pivot[col] + pivot[col] = temp + even = !even + } + + // Divide the lower elements by the "winning" diagonal elt. + val luDiag = lu[col, col] + for (row in col + 1 until m) { + lu[row, col] = lu[row, col] / luDiag + } + } + return LUPDecomposition(context.elementContext, lu, pivot, even) + } + } + + /** + * Produce a matrix with added decomposition feature + */ + fun decompose(matrix: Matrix): Matrix { + if (matrix.hasFeature>()) { + return matrix + } else { + val decomposition = buildDecomposition(matrix) + return VirtualMatrix.wrap(matrix, decomposition) + } + } + + + override fun solve(a: Matrix, b: Matrix): Matrix { + if (b.rowNum != a.colNum) { + error("Matrix dimension mismatch expected ${a.rowNum}, but got ${b.colNum}") + } + + // Use existing decomposition if it is provided by matrix + val decomposition = a.getFeature() ?: buildDecomposition(a) + + with(decomposition) { + with(context.elementContext) { + // Apply permutations to b + val bp = Mutable2DStructure.create(a.rowNum, a.colNum, bufferFactory) { i, j -> + b[pivot[i], j] + } + + // Solve LY = b + for (col in 0 until a.rowNum) { + for (i in col + 1 until a.rowNum) { + for (j in 0 until b.colNum) { + bp[i, j] -= bp[col, j] * lu[i, col] + } + } + } + + // Solve UX = Y + for (col in a.rowNum - 1 downTo 0) { + for (j in 0 until b.colNum) { + bp[col, j] /= lu[col, col] + } + for (i in 0 until col) { + for (j in 0 until b.colNum) { + bp[i, j] -= bp[col, j] * lu[i, col] + } + } + } + + return context.produce(a.rowNum, a.colNum) { i, j -> bp[i, j] } + } + } + } + + override fun inverse(a: Matrix): Matrix = solve(a, context.one(a.rowNum, a.colNum)) + + companion object { + val real = LUSolver(MatrixContext.real, MutableBuffer.Companion::auto) { it < 1e-11 } + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt index acd27dbd9..51f61e030 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt @@ -1,62 +1,45 @@ package scientifik.kmath.linear -import scientifik.kmath.operations.DoubleField import scientifik.kmath.operations.Field import scientifik.kmath.operations.Norm +import scientifik.kmath.operations.RealField +import scientifik.kmath.structures.VirtualBuffer +import scientifik.kmath.structures.asSequence /** * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors */ -interface LinearSolver> { - fun solve(a: Matrix, b: Matrix): Matrix - fun solve(a: Matrix, b: Vector): Vector = solve(a, b.toMatrix()).toVector() - fun inverse(a: Matrix): Matrix = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field)) +interface LinearSolver { + fun solve(a: Matrix, b: Matrix): Matrix + fun solve(a: Matrix, b: Point): Point = solve(a, b.toMatrix()).toVector() + fun inverse(a: Matrix): Matrix } /** * Convert vector to array (copying content of array) */ -fun Array.toVector(field: Field) = Vector.of(size, field) { this[it] } +fun Array.toVector(field: Field) = Vector.generic(size, field) { this[it] } + +fun DoubleArray.toVector() = Vector.real(this.size) { this[it] } +fun List.toVector() = Vector.real(this.size) { this[it] } + +object VectorL2Norm : Norm, Double> { + override fun norm(arg: Point): Double = + kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() }) +} + +typealias RealVector = Vector +typealias RealMatrix = Matrix + -fun DoubleArray.toVector() = Vector.ofReal(this.size) { this[it] } -fun List.toVector() = Vector.ofReal(this.size) { this[it] } /** * Convert matrix to vector if it is possible */ -fun > Matrix.toVector(): Vector { - return when { - this.columns == 1 -> { -// if (this is ArrayMatrix) { -// //Reuse existing underlying array -// ArrayVector(ArrayVectorSpace(rows, context.field, context.ndFactory), array) -// } else { -// //Generic vector -// vector(rows, context.field) { get(it, 0) } -// } - Vector.of(rows, context.field) { get(it, 0) } - } - else -> error("Can't convert matrix with more than one column to vector") - } -} +fun Matrix.toVector(): Point = + if (this.colNum == 1) { + VirtualBuffer(rowNum){ get(it, 0) } + } else error("Can't convert matrix with more than one column to vector") -fun > Vector.toMatrix(): Matrix { -// return if (this is ArrayVector) { -// //Reuse existing underlying array -// ArrayMatrix(ArrayMatrixSpace(size, 1, context.field, context.ndFactory), array) -// } else { -// //Generic vector -// matrix(size, 1, context.field) { i, j -> get(i) } -// } - return Matrix.of(size, 1, context.space) { i, _ -> get(i) } -} - -object VectorL2Norm : Norm, Double> { - override fun norm(arg: Vector): Double { - return kotlin.math.sqrt(arg.sumByDouble { it.toDouble() }) - } -} - -typealias RealVector = Vector -typealias RealMatrix = Matrix +fun Point.toMatrix(): Matrix = VirtualMatrix(size, 1) { i, _ -> get(i) } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt index 243eb2d6c..92acf2bbb 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt @@ -1,187 +1,214 @@ package scientifik.kmath.linear -import scientifik.kmath.operations.* +import scientifik.kmath.operations.RealField +import scientifik.kmath.operations.Ring +import scientifik.kmath.operations.sum import scientifik.kmath.structures.* +import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory +import scientifik.kmath.structures.Buffer.Companion.boxing +import kotlin.math.sqrt -/** - * The space for linear elements. Supports scalar product alongside with standard linear operations. - * @param T type of individual element of the vector or matrix - * @param V the type of vector space element - */ -abstract class MatrixSpace>(val rows: Int, val columns: Int, val field: F) : Space> { +interface MatrixContext { /** - * Produce the element of this space + * Produce a matrix with this context and given dimensions */ - abstract fun produce(initializer: (Int, Int) -> T): Matrix + fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix - /** - * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements - */ - abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace + infix fun Matrix.dot(other: Matrix): Matrix - override val zero: Matrix by lazy { - produce { _, _ -> field.zero } - } + infix fun Matrix.dot(vector: Point): Point -// val one: Matrix by lazy { -// produce { i, j -> if (i == j) field.one else field.zero } -// } + operator fun Matrix.unaryMinus(): Matrix - override fun add(a: Matrix, b: Matrix): Matrix { - return produce { i, j -> with(field) { a[i, j] + b[i, j] } } - } + operator fun Matrix.plus(b: Matrix): Matrix - override fun multiply(a: Matrix, k: Double): Matrix { - //TODO it is possible to implement scalable linear elements which normed values and adjustable scale to save memory and processing poser - return produce { i, j -> with(field) { a[i, j] * k } } - } + operator fun Matrix.minus(b: Matrix): Matrix - /** - * Dot product. Throws exception on dimension mismatch - */ - fun multiply(a: Matrix, b: Matrix): Matrix { - if (a.rows != b.columns) { - //TODO replace by specific exception - error("Dimension mismatch in linear structure dot product: [${a.rows},${a.columns}]*[${b.rows},${b.columns}]") - } - return produceSpace(a.rows, b.columns).produce { i, j -> - (0 until a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) } - } - } + operator fun Matrix.times(value: T): Matrix - override fun equals(other: Any?): Boolean { - if (this === other) return true - if (other !is MatrixSpace<*,*>) return false + operator fun T.times(m: Matrix): Matrix = m * this - if (rows != other.rows) return false - if (columns != other.columns) return false - if (field != other.field) return false + companion object { + /** + * Non-boxing double matrix + */ + val real = BufferMatrixContext(RealField, DoubleBufferFactory) - return true - } + /** + * A structured matrix with custom buffer + */ + fun > buffered( + ring: R, + bufferFactory: BufferFactory = ::boxing + ): GenericMatrixContext = + BufferMatrixContext(ring, bufferFactory) - override fun hashCode(): Int { - var result = rows - result = 31 * result + columns - result = 31 * result + field.hashCode() - return result + /** + * Automatic buffered matrix, unboxed if it is possible + */ + inline fun > auto(ring: R): GenericMatrixContext = + buffered(ring, Buffer.Companion::auto) } } -infix fun > Matrix.dot(b: Matrix): Matrix = this.context.multiply(this, b) +interface GenericMatrixContext> : MatrixContext { + /** + * The ring context for matrix elements + */ + val elementContext: R + + /** + * Produce a point compatible with matrix space + */ + fun point(size: Int, initializer: (Int) -> T): Point + + override infix fun Matrix.dot(other: Matrix): Matrix { + //TODO add typed error + if (this.colNum != other.rowNum) error("Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})") + return produce(rowNum, other.colNum) { i, j -> + val row = rows[i] + val column = other.columns[j] + with(elementContext) { + sum(row.asSequence().zip(column.asSequence(), ::multiply)) + } + } + } + + override infix fun Matrix.dot(vector: Point): Point { + //TODO add typed error + if (this.colNum != vector.size) error("Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})") + return point(rowNum) { i -> + val row = rows[i] + with(elementContext) { + sum(row.asSequence().zip(vector.asSequence(), ::multiply)) + } + } + } + + override operator fun Matrix.unaryMinus() = + produce(rowNum, colNum) { i, j -> elementContext.run { -get(i, j) } } + + override operator fun Matrix.plus(b: Matrix): Matrix { + if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] + [${b.rowNum},${b.colNum}]") + return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } } + } + + override operator fun Matrix.minus(b: Matrix): Matrix { + if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]") + return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } } + } + + operator fun Matrix.times(number: Number): Matrix = + produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * number } } + + operator fun Number.times(matrix: Matrix): Matrix = matrix * this + + override fun Matrix.times(value: T): Matrix = + produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } } +} /** - * A matrix-like structure + * Specialized 2-d structure */ -interface Matrix> : SpaceElement, MatrixSpace> { - /** - * Number of rows - */ - val rows: Int - /** - * Number of columns - */ - val columns: Int +interface Matrix : Structure2D { + val rowNum: Int + val colNum: Int + + val features: Set /** - * Get element in row [i] and column [j]. Throws error in case of call ounside structure dimensions + * Suggest new feature for this matrix. The result is the new matrix that may or may not reuse existing data structure. + * + * The implementation does not guarantee to check that matrix actually have the feature, so one should be careful to + * add only those features that are valid. */ - operator fun get(i: Int, j: Int): T + fun suggestFeature(vararg features: MatrixFeature): Matrix - override val self: Matrix - get() = this + override fun get(index: IntArray): T = get(index[0], index[1]) - fun transpose(): Matrix { - return object : Matrix { - override val context: MatrixSpace = this@Matrix.context - override val rows: Int = this@Matrix.columns - override val columns: Int = this@Matrix.rows - override fun get(i: Int, j: Int): T = this@Matrix[j, i] + override val shape: IntArray get() = intArrayOf(rowNum, colNum) + + val rows: Point> + get() = VirtualBuffer(rowNum) { i -> + VirtualBuffer(colNum) { j -> get(i, j) } + } + + val columns: Point> + get() = VirtualBuffer(colNum) { j -> + VirtualBuffer(rowNum) { i -> get(i, j) } + } + + override fun elements(): Sequence> = sequence { + for (i in (0 until rowNum)) { + for (j in (0 until colNum)) { + yield(intArrayOf(i, j) to get(i, j)) + } } } companion object { + fun real(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = + MatrixContext.real.produce(rows, columns, initializer) /** - * Create [ArrayMatrix] with custom field + * Build a square matrix from given elements. */ - fun > of(rows: Int, columns: Int, field: F, initializer: (Int, Int) -> T) = - ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer) - - /** - * Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing. - */ - fun ofReal(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = - ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer) - - /** - * Create a diagonal value matrix. By default value equals [Field.one]. - */ - fun > diagonal(rows: Int, columns: Int, field: F, values: (Int) -> T = { field.one }): Matrix { - return of(rows, columns, field) { i, j -> if (i == j) values(i) else field.zero } + fun square(vararg elements: T): Matrix { + val size: Int = sqrt(elements.size.toDouble()).toInt() + if (size * size != elements.size) error("The number of elements ${elements.size} is not a full square") + val buffer = elements.asBuffer() + return BufferMatrix(size, size, buffer) } - /** - * Equality check on two generic matrices - */ - fun equals(mat1: Matrix<*, *>, mat2: Matrix<*, *>): Boolean { - if (mat1 === mat2) return true - if (mat1.context != mat2.context) return false - for (i in 0 until mat1.rows) { - for (j in 0 until mat2.columns) { - if (mat1[i, j] != mat2[i, j]) return false - } - } - return true - } + fun build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns) } } - - - -typealias NDFieldFactory = (IntArray) -> NDField - -internal fun > genericNDFieldFactory(field: F): NDFieldFactory = { index -> GenericNDField(index, field) } -internal val realNDFieldFactory: NDFieldFactory = { index -> ExtendedNDField(index, DoubleField) } - - -/** - * NDArray-based implementation of vector space. By default uses slow [GenericNDField], but could be overridden with custom [NDField] factory. - */ -class ArrayMatrixSpace>( - rows: Int, - columns: Int, - field: F, - val ndFactory: NDFieldFactory = genericNDFieldFactory(field) -) : MatrixSpace(rows, columns, field) { - - val ndField by lazy { - ndFactory(intArrayOf(rows, columns)) - } - - override fun produce(initializer: (Int, Int) -> T): Matrix = ArrayMatrix(this, initializer) - - override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace { - return ArrayMatrixSpace(rows, columns, field, ndFactory) +class MatrixBuilder(val rows: Int, val columns: Int) { + operator fun invoke(vararg elements: T): Matrix { + 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) } } /** - * Member of [ArrayMatrixSpace] which wraps 2-D array + * Check if matrix has the given feature class */ -class ArrayMatrix> internal constructor(override val context: ArrayMatrixSpace, val element: NDElement) : Matrix { +inline fun Matrix<*>.hasFeature(): Boolean = features.find { it is T } != null - constructor(context: ArrayMatrixSpace, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) }) +/** + * Get the first feature matching given class. Does not guarantee that matrix has only one feature matching the criteria + */ +inline fun Matrix<*>.getFeature(): T? = features.filterIsInstance().firstOrNull() - override val rows: Int get() = context.rows - - override val columns: Int get() = context.columns - - override fun get(i: Int, j: Int): T { - return element[i, j] +/** + * Diagonal matrix of ones. The matrix is virtual no actual matrix is created + */ +fun > GenericMatrixContext.one(rows: Int, columns: Int): Matrix = + VirtualMatrix(rows, columns) { i, j -> + if (i == j) elementContext.one else elementContext.zero } - override val self: ArrayMatrix get() = this + +/** + * A virtual matrix of zeroes + */ +fun > GenericMatrixContext.zero(rows: Int, columns: Int): Matrix = + VirtualMatrix(rows, columns) { _, _ -> elementContext.zero } + +class TransposedFeature(val original: Matrix) : MatrixFeature + +/** + * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A` + */ +fun > Matrix.transpose(): Matrix { + return this.getFeature>()?.original ?: VirtualMatrix( + this.colNum, + this.rowNum, + setOf(TransposedFeature(this)) + ) { i, j -> get(j, i) } } + +infix fun Matrix.dot(other: Matrix): Matrix = with(MatrixContext.real) { dot(other) } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt new file mode 100644 index 000000000..6b45a14b1 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixFeatures.kt @@ -0,0 +1,62 @@ +package scientifik.kmath.linear + +/** + * A marker interface representing some matrix feature like diagonal, sparce, zero, etc. Features used to optimize matrix + * operations performance in some cases. + */ +interface MatrixFeature + +/** + * The matrix with this feature is considered to have only diagonal non-null elements + */ +object DiagonalFeature : MatrixFeature + +/** + * Matrix with this feature has all zero elements + */ +object ZeroFeature : MatrixFeature + +/** + * Matrix with this feature have unit elements on diagonal and zero elements in all other places + */ +object UnitFeature : MatrixFeature + +/** + * Inverted matrix feature + */ +interface InverseMatrixFeature : MatrixFeature { + val inverse: Matrix +} + +/** + * A determinant container + */ +interface DeterminantFeature : MatrixFeature { + val determinant: T +} + +@Suppress("FunctionName") +fun DeterminantFeature(determinant: T) = object: DeterminantFeature{ + override val determinant: T = determinant +} + +/** + * Lower triangular matrix + */ +object LFeature: MatrixFeature + +/** + * Upper triangular feature + */ +object UFeature: MatrixFeature + +/** + * TODO add documentation + */ +interface LUPDecompositionFeature : MatrixFeature { + val l: Matrix + val u: Matrix + val p: Matrix +} + +//TODO add sparse matrix feature \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Mutable2DStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Mutable2DStructure.kt new file mode 100644 index 000000000..0d35c57e6 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Mutable2DStructure.kt @@ -0,0 +1,40 @@ +package scientifik.kmath.linear + +import scientifik.kmath.structures.MutableBuffer +import scientifik.kmath.structures.MutableBufferFactory +import scientifik.kmath.structures.MutableNDStructure + +class Mutable2DStructure(val rowNum: Int, val colNum: Int, val buffer: MutableBuffer) : MutableNDStructure { + override val shape: IntArray + get() = intArrayOf(rowNum, colNum) + + operator fun get(i: Int, j: Int): T = buffer[i * colNum + j] + + override fun get(index: IntArray): T = get(index[0], index[1]) + + override fun elements(): Sequence> = sequence { + for (i in 0 until rowNum) { + for (j in 0 until colNum) { + yield(intArrayOf(i, j) to get(i, j)) + } + } + } + + operator fun set(i: Int, j: Int, value: T) { + buffer[i * colNum + j] = value + } + + override fun set(index: IntArray, value: T) = set(index[0], index[1], value) + + companion object { + fun create( + rowNum: Int, + colNum: Int, + bufferFactory: MutableBufferFactory, + init: (i: Int, j: Int) -> T + ): Mutable2DStructure { + val buffer = bufferFactory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) } + return Mutable2DStructure(rowNum, colNum, buffer) + } + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt index 5ab9907de..a7373f647 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt @@ -1,103 +1,125 @@ package scientifik.kmath.linear -import scientifik.kmath.histogram.Point -import scientifik.kmath.operations.DoubleField -import scientifik.kmath.operations.Field +import scientifik.kmath.operations.RealField import scientifik.kmath.operations.Space import scientifik.kmath.operations.SpaceElement -import scientifik.kmath.structures.NDElement -import scientifik.kmath.structures.get +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.BufferFactory +import scientifik.kmath.structures.asSequence + +typealias Point = Buffer /** * A linear space for vectors. * Could be used on any point-like structure */ -abstract class VectorSpace>(val size: Int, val space: S) : Space> { +interface VectorSpace> : Space> { - abstract fun produce(initializer: (Int) -> T): Vector + val size: Int - override val zero: Vector by lazy { produce { space.zero } } + val space: S - override fun add(a: Point, b: Point): Vector = produce { with(space) { a[it] + b[it] } } + fun produce(initializer: (Int) -> T): Point - override fun multiply(a: Point, k: Double): Vector = produce { with(space) { a[it] * k } } + /** + * Produce a space-element of this vector space for expressions + */ + fun produceElement(initializer: (Int) -> T): Vector + + override val zero: Point get() = produce { space.zero } + + override fun add(a: Point, b: Point): Point = produce { with(space) { a[it] + b[it] } } + + override fun multiply(a: Point, k: Number): Point = produce { with(space) { a[it] * k } } + + //TODO add basis + + companion object { + + private val realSpaceCache = HashMap>() + + /** + * Non-boxing double vector space + */ + fun real(size: Int): BufferVectorSpace { + return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, RealField, Buffer.DoubleBufferFactory) } + } + + /** + * A structured vector space with custom buffer + */ + fun > buffered( + size: Int, + space: S, + bufferFactory: BufferFactory = Buffer.Companion::boxing + ): VectorSpace = BufferVectorSpace(size, space, bufferFactory) + + /** + * Automatic buffered vector, unboxed if it is possible + */ + inline fun > smart(size: Int, space: S): VectorSpace = + buffered(size, space, Buffer.Companion::auto) + } } /** * A point coupled to the linear space */ -interface Vector> : SpaceElement, VectorSpace>, Point, Iterable { +interface Vector> : SpaceElement, Vector, VectorSpace>, Point { override val size: Int get() = context.size - override operator fun plus(b: Point): Vector = context.add(self, b) - override operator fun minus(b: Point): Vector = context.add(self, context.multiply(b, -1.0)) - override operator fun times(k: Number): Vector = context.multiply(self, k.toDouble()) - override operator fun div(k: Number): Vector = context.multiply(self, 1.0 / k.toDouble()) + override operator fun plus(b: Point): Vector = context.add(this, b).wrap() + override operator fun minus(b: Point): Vector = context.add(this, context.multiply(b, -1.0)).wrap() + override operator fun times(k: Number): Vector = context.multiply(this, k.toDouble()).wrap() + override operator fun div(k: Number): Vector = context.multiply(this, 1.0 / k.toDouble()).wrap() companion object { /** * Create vector with custom field */ - fun > of(size: Int, field: F, initializer: (Int) -> T) = - ArrayVector(ArrayVectorSpace(size, field), initializer) + fun > generic(size: Int, field: S, initializer: (Int) -> T): Vector = + VectorSpace.buffered(size, field).produceElement(initializer) - private val realSpaceCache = HashMap>() + fun real(size: Int, initializer: (Int) -> Double): Vector = + VectorSpace.real(size).produceElement(initializer) - private fun getRealSpace(size: Int): ArrayVectorSpace { - return realSpaceCache.getOrPut(size){ArrayVectorSpace(size, DoubleField, realNDFieldFactory)} - } + fun ofReal(vararg elements: Double): Vector = + VectorSpace.real(elements.size).produceElement { elements[it] } - /** - * Create vector of [Double] - */ - fun ofReal(size: Int, initializer: (Int) -> Double) = - ArrayVector(getRealSpace(size), initializer) - - fun ofReal(vararg point: Double) = point.toVector() - - fun equals(v1: Vector<*, *>, v2: Vector<*, *>): Boolean { - if (v1 === v2) return true - if (v1.context != v2.context) return false - for (i in 0 until v2.size) { - if (v1[i] != v2[i]) return false - } - return true - } } } -class ArrayVectorSpace>( - size: Int, - field: F, - val ndFactory: NDFieldFactory = genericNDFieldFactory(field) -) : VectorSpace(size, field) { - val ndField by lazy { - ndFactory(intArrayOf(size)) - } - - override fun produce(initializer: (Int) -> T): Vector = ArrayVector(this, initializer) +data class BufferVectorSpace>( + override val size: Int, + override val space: S, + val bufferFactory: BufferFactory +) : VectorSpace { + override fun produce(initializer: (Int) -> T) = bufferFactory(size, initializer) + override fun produceElement(initializer: (Int) -> T): Vector = BufferVector(this, produce(initializer)) } -class ArrayVector> internal constructor(override val context: VectorSpace, val element: NDElement) : Vector { - - constructor(context: ArrayVectorSpace, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) }) +data class BufferVector>(override val context: VectorSpace, val buffer: Buffer) : + Vector { init { - if (context.size != element.shape[0]) { + if (context.size != buffer.size) { error("Array dimension mismatch") } } override fun get(index: Int): T { - return element[index] + return buffer[index] } - override val self: ArrayVector get() = this + override fun unwrap(): Point = this - override fun iterator(): Iterator = (0 until size).map { element[it] }.iterator() + override fun Point.wrap(): Vector = BufferVector(context, this) - override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() } + override fun iterator(): Iterator = (0 until size).map { buffer[it] }.iterator() + + override fun toString(): String = + this.asSequence().joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() } } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt new file mode 100644 index 000000000..1bab52902 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VirtualMatrix.kt @@ -0,0 +1,47 @@ +package scientifik.kmath.linear + +class VirtualMatrix( + override val rowNum: Int, + override val colNum: Int, + override val features: Set = emptySet(), + val generator: (i: Int, j: Int) -> T +) : Matrix { + override fun get(i: Int, j: Int): T = generator(i, j) + + override fun suggestFeature(vararg features: MatrixFeature) = + VirtualMatrix(rowNum, colNum, this.features + features, generator) + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other !is Matrix<*>) return false + + if (rowNum != other.rowNum) return false + if (colNum != other.colNum) return false + + return elements().all { (index, value) -> value == other[index] } + } + + override fun hashCode(): Int { + var result = rowNum + result = 31 * result + colNum + result = 31 * result + features.hashCode() + result = 31 * result + generator.hashCode() + return result + } + + + companion object { + /** + * Wrap a matrix adding additional features to it + */ + fun wrap(matrix: Matrix, vararg features: MatrixFeature): Matrix { + return if (matrix is VirtualMatrix) { + VirtualMatrix(matrix.rowNum, matrix.colNum, matrix.features + features, matrix.generator) + } else { + VirtualMatrix(matrix.rowNum, matrix.colNum, matrix.features + features) { i, j -> + matrix[i, j] + } + } + } + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Cumulative.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Cumulative.kt deleted file mode 100644 index 4d4f8ced6..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Cumulative.kt +++ /dev/null @@ -1,59 +0,0 @@ -package scientifik.kmath.misc - -import kotlin.jvm.JvmName - - -/** - * Generic cumulative operation on iterator - * @param T type of initial iterable - * @param R type of resulting iterable - * @param initial lazy evaluated - */ -fun Iterator.cumulative(initial: R, operation: (T, R) -> R): Iterator = object : Iterator { - var state: R = initial - override fun hasNext(): Boolean = this@cumulative.hasNext() - - override fun next(): R { - state = operation.invoke(this@cumulative.next(), state) - return state - } -} - -fun Iterable.cumulative(initial: R, operation: (T, R) -> R): Iterable = object : Iterable { - override fun iterator(): Iterator = this@cumulative.iterator().cumulative(initial, operation) -} - -fun Sequence.cumulative(initial: R, operation: (T, R) -> R): Sequence = object : Sequence { - override fun iterator(): Iterator = this@cumulative.iterator().cumulative(initial, operation) -} - -fun List.cumulative(initial: R, operation: (T, R) -> R): List = this.iterator().cumulative(initial, operation).asSequence().toList() - -//Cumulative sum - -@JvmName("cumulativeSumOfDouble") -fun Iterable.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfInt") -fun Iterable.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfLong") -fun Iterable.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfDouble") -fun Sequence.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfInt") -fun Sequence.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfLong") -fun Sequence.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfDouble") -fun List.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfInt") -fun List.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element} - -@JvmName("cumulativeSumOfLong") -fun List.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt index b48afbdcc..90ce5da68 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt @@ -8,30 +8,29 @@ package scientifik.kmath.misc * * If step is negative, the same goes from upper boundary downwards */ -fun ClosedFloatingPointRange.toSequence(step: Double): Sequence { - return when { - step == 0.0 -> error("Zero step in double progression") - step > 0 -> sequence { - var current = start - while (current <= endInclusive) { - yield(current) - current += step +fun ClosedFloatingPointRange.toSequence(step: Double): Sequence = + when { + step == 0.0 -> error("Zero step in double progression") + step > 0 -> sequence { + var current = start + while (current <= endInclusive) { + yield(current) + current += step + } + } + else -> sequence { + var current = endInclusive + while (current >= start) { + yield(current) + current += step + } } } - else -> sequence { - var current = endInclusive - while (current >= start) { - yield(current) - current += step - } - } - } -} /** * Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints] */ fun ClosedFloatingPointRange.toGrid(numPoints: Int): DoubleArray { - if (numPoints < 2) error("Can't create grid with less than two points") + if (numPoints < 2) error("Can't create generic grid with less than two points") return DoubleArray(numPoints) { i -> start + (endInclusive - start) / (numPoints - 1) * i } } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt index 9a56a4e91..ee9833623 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt @@ -1,37 +1,7 @@ package scientifik.kmath.operations -/** - * The generic mathematics elements which is able to store its context - * @param T the self type of the element - * @param S the type of mathematical context for this element - */ -interface MathElement { - /** - * Self value. Needed for static type checking. - */ - val self: T - - /** - * The context this element belongs to - */ - val context: S -} - -/** - * A general interface representing linear context of some kind. - * The context defines sum operation for its elements and multiplication by real value. - * One must note that in some cases context is a singleton class, but in some cases it - * works as a context for operations inside it. - * - * TODO do we need commutative context? - */ -interface Space { - /** - * Neutral element for sum operation - */ - val zero: T - +interface SpaceOperations { /** * Addition operation for two context elements */ @@ -40,7 +10,7 @@ interface Space { /** * Multiplication operation for context element and real number */ - fun multiply(a: T, k: Double): T + fun multiply(a: T, k: Number): T //Operation to be performed in this context operator fun T.unaryMinus(): T = multiply(this, -1.0) @@ -50,69 +20,61 @@ interface Space { operator fun T.times(k: Number) = multiply(this, k.toDouble()) operator fun T.div(k: Number) = multiply(this, 1.0 / k.toDouble()) operator fun Number.times(b: T) = b * this - } -/** - * The element of linear context - * @param T self type of the element. Needed for static type checking - * @param S the type of space - */ -interface SpaceElement> : MathElement { - operator fun plus(b: T): T = context.add(self, b) - operator fun minus(b: T): T = context.add(self, context.multiply(b, -1.0)) - operator fun times(k: Number): T = context.multiply(self, k.toDouble()) - operator fun div(k: Number): T = context.multiply(self, 1.0 / k.toDouble()) -} /** - * The same as {@link Space} but with additional multiplication operation + * A general interface representing linear context of some kind. + * The context defines sum operation for its elements and multiplication by real value. + * One must note that in some cases context is a singleton class, but in some cases it + * works as a context for operations inside it. + * + * TODO do we need non-commutative context? */ -interface Ring : Space { +interface Space : SpaceOperations { /** - * neutral operation for multiplication + * Neutral element for sum operation */ - val one: T + val zero: T +} +interface RingOperations : SpaceOperations { /** * Multiplication for two field elements */ fun multiply(a: T, b: T): T operator fun T.times(b: T): T = multiply(this, b) - } /** - * Ring element + * The same as {@link Space} but with additional multiplication operation */ -interface RingElement> : SpaceElement { - override val context: S +interface Ring : Space, RingOperations { + /** + * neutral operation for multiplication + */ + val one: T - operator fun times(b: T): T = context.multiply(self, b) +// operator fun T.plus(b: Number) = this.plus(b * one) +// operator fun Number.plus(b: T) = b + this +// +// operator fun T.minus(b: Number) = this.minus(b * one) +// operator fun Number.minus(b: T) = -b + this +} + +/** + * All ring operations but without neutral elements + */ +interface FieldOperations : RingOperations { + fun divide(a: T, b: T): T + + operator fun T.div(b: T): T = divide(this, b) } /** * Four operations algebra */ -interface Field : Ring { - fun divide(a: T, b: T): T - - operator fun T.div(b: T): T = divide(this, b) +interface Field : Ring, FieldOperations { operator fun Number.div(b: T) = this * divide(one, b) - - operator fun T.plus(b: Number) = this.plus(b * one) - operator fun Number.plus(b: T) = b + this - - operator fun T.minus(b: Number) = this.minus(b * one) - operator fun Number.minus(b: T) = -b + this } - -/** - * Field element - */ -interface FieldElement> : RingElement { - override val context: S - - operator fun div(b: T): T = context.divide(self, b) -} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraElements.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraElements.kt new file mode 100644 index 000000000..093021ae3 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraElements.kt @@ -0,0 +1,48 @@ +package scientifik.kmath.operations + +/** + * The generic mathematics elements which is able to store its context + * @param T the type of space operation results + * @param I self type of the element. Needed for static type checking + * @param C the type of mathematical context for this element + */ +interface MathElement { + /** + * The context this element belongs to + */ + val context: C +} + +interface MathWrapper { + fun unwrap(): T + fun T.wrap(): I +} + +/** + * The element of linear context + * @param T the type of space operation results + * @param I self type of the element. Needed for static type checking + * @param S the type of space + */ +interface SpaceElement, S : Space> : MathElement, MathWrapper { + + operator fun plus(b: T) = context.add(unwrap(), b).wrap() + operator fun minus(b: T) = context.add(unwrap(), context.multiply(b, -1.0)).wrap() + operator fun times(k: Number) = context.multiply(unwrap(), k.toDouble()).wrap() + operator fun div(k: Number) = context.multiply(unwrap(), 1.0 / k.toDouble()).wrap() +} + +/** + * Ring element + */ +interface RingElement, R : Ring> : SpaceElement { + operator fun times(b: T) = context.multiply(unwrap(), b).wrap() +} + +/** + * Field element + */ +interface FieldElement, F : Field> : RingElement { + override val context: F + operator fun div(b: T) = context.divide(unwrap(), b).wrap() +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt new file mode 100644 index 000000000..6dec8bd79 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt @@ -0,0 +1,7 @@ +package scientifik.kmath.operations + +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.asSequence + +fun Space.sum(data : Iterable): T = data.fold(zero) { left, right -> add(left,right) } +fun Space.sum(data : Sequence): T = data.fold(zero) { left, right -> add(left, right) } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt index 8d3c01d28..4b4002cbb 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt @@ -1,45 +1,97 @@ package scientifik.kmath.operations +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.MutableBuffer +import scientifik.kmath.structures.ObjectBuffer +import scientifik.memory.MemoryReader +import scientifik.memory.MemorySpec +import scientifik.memory.MemoryWriter +import kotlin.math.* + /** * A field for complex numbers */ -object ComplexField : Field { +object ComplexField : ExtendedFieldOperations, Field { override val zero: Complex = Complex(0.0, 0.0) - override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im) - - override fun multiply(a: Complex, k: Double): Complex = Complex(a.re * k, a.im * k) - override val one: Complex = Complex(1.0, 0.0) - override fun multiply(a: Complex, b: Complex): Complex = Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re) + val i = Complex(0.0, 1.0) - override fun divide(a: Complex, b: Complex): Complex = Complex(a.re * b.re + a.im * b.im, a.re * b.im - a.im * b.re) / b.square + override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im) + override fun multiply(a: Complex, k: Number): Complex = Complex(a.re * k.toDouble(), a.im * k.toDouble()) + + override fun multiply(a: Complex, b: Complex): Complex = + Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re) + + override fun divide(a: Complex, b: Complex): Complex { + val norm = b.square + return Complex((a.re * b.re + a.im * b.im) / norm, (a.re * b.im - a.im * b.re) / norm) + } + + override fun sin(arg: Complex): Complex = i / 2 * (exp(-i * arg) - exp(i * arg)) + + override fun cos(arg: Complex): Complex = (exp(-i * arg) + exp(i * arg)) / 2 + + override fun power(arg: Complex, pow: Number): Complex = + arg.abs.pow(pow.toDouble()) * (cos(pow.toDouble() * arg.theta) + i * sin(pow.toDouble() * arg.theta)) + + override fun exp(arg: Complex): Complex = exp(arg.re) * (cos(arg.im) + i * sin(arg.im)) + + override fun ln(arg: Complex): Complex = ln(arg.abs) + i * atan2(arg.im, arg.re) + + operator fun Double.plus(c: Complex) = add(this.toComplex(), c) + + operator fun Double.minus(c: Complex) = add(this.toComplex(), -c) + + operator fun Complex.plus(d: Double) = d + this + + operator fun Complex.minus(d: Double) = add(this, -d.toComplex()) + + operator fun Double.times(c: Complex) = Complex(c.re * this, c.im * this) } /** * Complex number class */ -data class Complex(val re: Double, val im: Double) : FieldElement { - override val self: Complex get() = this - override val context: ComplexField - get() = ComplexField +data class Complex(val re: Double, val im: Double) : FieldElement { + override fun unwrap(): Complex = this + + override fun Complex.wrap(): Complex = this + + override val context: ComplexField get() = ComplexField /** * A complex conjugate */ - val conjugate: Complex - get() = Complex(re, -im) + val conjugate: Complex get() = Complex(re, -im) - val square: Double - get() = re * re + im * im + val square: Double get() = re * re + im * im - val abs: Double - get() = kotlin.math.sqrt(square) + val abs: Double get() = sqrt(square) - companion object { + val theta: Double get() = atan(im / re) + companion object : MemorySpec { + override val objectSize: Int = 16 + + override fun MemoryReader.read(offset: Int): Complex = + Complex(readDouble(offset), readDouble(offset + 8)) + + override fun MemoryWriter.write(offset: Int, value: Complex) { + writeDouble(offset, value.re) + writeDouble(offset + 8, value.im) + } } +} -} \ No newline at end of file +fun Double.toComplex() = Complex(this, 0.0) + +fun Buffer.Companion.complex(size: Int, init: (Int) -> Complex): Buffer { + return ObjectBuffer.create(Complex, size, init) +} + +fun MutableBuffer.Companion.complex(size: Int, init: (Int) -> Complex): Buffer { + return ObjectBuffer.create(Complex, size, init) +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt deleted file mode 100644 index d993d37b4..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt +++ /dev/null @@ -1,89 +0,0 @@ -package scientifik.kmath.operations - -import kotlin.math.pow - -/** - * Advanced Number-like field that implements basic operations - */ -interface ExtendedField : - Field, - TrigonometricOperations, - PowerOperations, - ExponentialOperations - - -/** - * Field for real values - */ -object RealField : ExtendedField, Norm { - override val zero: Real = Real(0.0) - override fun add(a: Real, b: Real): Real = Real(a.value + b.value) - override val one: Real = Real(1.0) - override fun multiply(a: Real, b: Real): Real = Real(a.value * b.value) - override fun multiply(a: Real, k: Double): Real = Real(a.value * k) - override fun divide(a: Real, b: Real): Real = Real(a.value / b.value) - - override fun sin(arg: Real): Real = Real(kotlin.math.sin(arg.value)) - override fun cos(arg: Real): Real = Real(kotlin.math.cos(arg.value)) - - override fun power(arg: Real, pow: Double): Real = Real(arg.value.pow(pow)) - - override fun exp(arg: Real): Real = Real(kotlin.math.exp(arg.value)) - - override fun ln(arg: Real): Real = Real(kotlin.math.ln(arg.value)) - - override fun norm(arg: Real): Real = Real(kotlin.math.abs(arg.value)) -} - -/** - * Real field element wrapping double. - * - * TODO inline does not work due to compiler bug. Waiting for fix for KT-27586 - */ -inline class Real(val value: Double) : FieldElement { - - //values are dynamically calculated to save memory - override val self - get() = this - - override val context - get() = RealField - - companion object { - - } -} - -/** - * A field for double without boxing. Does not produce appropriate field element - */ -object DoubleField : ExtendedField, Norm { - override val zero: Double = 0.0 - override fun add(a: Double, b: Double): Double = a + b - override fun multiply(a: Double, @Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE") b: Double): Double = a * b - override val one: Double = 1.0 - override fun divide(a: Double, b: Double): Double = a / b - - override fun sin(arg: Double): Double = kotlin.math.sin(arg) - override fun cos(arg: Double): Double = kotlin.math.cos(arg) - - override fun power(arg: Double, pow: Double): Double = arg.pow(pow) - - override fun exp(arg: Double): Double = kotlin.math.exp(arg) - - override fun ln(arg: Double): Double = kotlin.math.ln(arg) - - override fun norm(arg: Double): Double = kotlin.math.abs(arg) -} - -/** - * A field for double without boxing. Does not produce appropriate field element - */ -object IntField : Field{ - override val zero: Int = 0 - override fun add(a: Int, b: Int): Int = a + b - override fun multiply(a: Int, b: Int): Int = a * b - override fun multiply(a: Int, k: Double): Int = (k*a).toInt() - override val one: Int = 1 - override fun divide(a: Int, b: Int): Int = a / b -} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt new file mode 100644 index 000000000..b80212375 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/NumberAlgebra.kt @@ -0,0 +1,109 @@ +package scientifik.kmath.operations + +import kotlin.math.pow + +/** + * Advanced Number-like field that implements basic operations + */ +interface ExtendedFieldOperations : + FieldOperations, + TrigonometricOperations, + PowerOperations, + ExponentialOperations + +interface ExtendedField : ExtendedFieldOperations, Field + +/** + * Real field element wrapping double. + * + * TODO inline does not work due to compiler bug. Waiting for fix for KT-27586 + */ +inline class Real(val value: Double) : FieldElement { + override fun unwrap(): Double = value + + override fun Double.wrap(): Real = Real(value) + + override val context get() = RealField + + companion object +} + +/** + * A field for double without boxing. Does not produce appropriate field element + */ +@Suppress("EXTENSION_SHADOWED_BY_MEMBER") +object RealField : Field, ExtendedFieldOperations, Norm { + override val zero: Double = 0.0 + override fun add(a: Double, b: Double): Double = a + b + override fun multiply(a: Double, b: Double): Double = a * b + override fun multiply(a: Double, k: Number): Double = a * k.toDouble() + + override val one: Double = 1.0 + override fun divide(a: Double, b: Double): Double = a / b + + override fun sin(arg: Double): Double = kotlin.math.sin(arg) + override fun cos(arg: Double): Double = kotlin.math.cos(arg) + + override fun power(arg: Double, pow: Number): Double = arg.pow(pow.toDouble()) + + override fun exp(arg: Double): Double = kotlin.math.exp(arg) + override fun ln(arg: Double): Double = kotlin.math.ln(arg) + + override fun norm(arg: Double): Double = kotlin.math.abs(arg) + + override fun Double.unaryMinus(): Double = -this + + override fun Double.minus(b: Double): Double = this - b +} + +/** + * A field for [Int] without boxing. Does not produce corresponding field element + */ +object IntRing : Ring, Norm { + override val zero: Int = 0 + override fun add(a: Int, b: Int): Int = a + b + override fun multiply(a: Int, b: Int): Int = a * b + override fun multiply(a: Int, k: Number): Int = (k * a) + override val one: Int = 1 + + override fun norm(arg: Int): Int = arg +} + +/** + * A field for [Short] without boxing. Does not produce appropriate field element + */ +object ShortRing : Ring, Norm { + override val zero: Short = 0 + override fun add(a: Short, b: Short): Short = (a + b).toShort() + override fun multiply(a: Short, b: Short): Short = (a * b).toShort() + override fun multiply(a: Short, k: Number): Short = (a * k) + override val one: Short = 1 + + override fun norm(arg: Short): Short = arg +} + +/** + * A field for [Byte] values + */ +object ByteRing : Ring, Norm { + override val zero: Byte = 0 + override fun add(a: Byte, b: Byte): Byte = (a + b).toByte() + override fun multiply(a: Byte, b: Byte): Byte = (a * b).toByte() + override fun multiply(a: Byte, k: Number): Byte = (a * k) + override val one: Byte = 1 + + override fun norm(arg: Byte): Byte = arg +} + +/** + * A field for [Long] values + */ +object LongRing : Ring, Norm { + override val zero: Long = 0 + override fun add(a: Long, b: Long): Long = (a + b) + override fun multiply(a: Long, b: Long): Long = (a * b) + override fun multiply(a: Long, k: Number): Long = (a * k) + override val one: Long = 1 + + override fun norm(arg: Long): Long = arg +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt index 9e5c8a801..ffbbf69dd 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt @@ -10,7 +10,7 @@ package scientifik.kmath.operations * It also allows to override behavior for optional operations * */ -interface TrigonometricOperations : Field { +interface TrigonometricOperations : FieldOperations { fun sin(arg: T): T fun cos(arg: T): T @@ -19,10 +19,10 @@ interface TrigonometricOperations : Field { fun ctg(arg: T): T = cos(arg) / sin(arg) } -fun >> sin(arg: T): T = arg.context.sin(arg) -fun >> cos(arg: T): T = arg.context.cos(arg) -fun >> tg(arg: T): T = arg.context.tg(arg) -fun >> ctg(arg: T): T = arg.context.ctg(arg) +fun >> sin(arg: T): T = arg.context.sin(arg) +fun >> cos(arg: T): T = arg.context.cos(arg) +fun >> tg(arg: T): T = arg.context.tg(arg) +fun >> ctg(arg: T): T = arg.context.ctg(arg) /* Power and roots */ @@ -30,12 +30,13 @@ fun >> ctg(arg: T): T = arg.c * A context extension to include power operations like square roots, etc */ interface PowerOperations { - fun power(arg: T, pow: Double): T + fun power(arg: T, pow: Number): T + fun sqrt(arg: T) = power(arg, 0.5) } -infix fun >> T.pow(power: Double): T = context.power(this, power) -fun >> sqrt(arg: T): T = arg pow 0.5 -fun >> sqr(arg: T): T = arg pow 2.0 +infix fun >> T.pow(power: Double): T = context.power(this, power) +fun >> sqrt(arg: T): T = arg pow 0.5 +fun >> sqr(arg: T): T = arg pow 2.0 /* Exponential */ @@ -44,11 +45,11 @@ interface ExponentialOperations { fun ln(arg: T): T } -fun >> exp(arg: T): T = arg.context.exp(arg) -fun >> ln(arg: T): T = arg.context.ln(arg) +fun >> exp(arg: T): T = arg.context.exp(arg) +fun >> ln(arg: T): T = arg.context.ln(arg) -interface Norm { +interface Norm { fun norm(arg: T): R } -fun >, R> norm(arg: T): R = arg.context.norm(arg) \ No newline at end of file +fun >, R> norm(arg: T): R = arg.context.norm(arg) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BoxingNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BoxingNDField.kt new file mode 100644 index 000000000..a73c58939 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BoxingNDField.kt @@ -0,0 +1,73 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.Field +import scientifik.kmath.operations.FieldElement + + +class BoxingNDField>( + override val shape: IntArray, + override val elementContext: F, + val bufferFactory: BufferFactory +) : BufferedNDField { + + override val strides: Strides = DefaultStrides(shape) + + fun buildBuffer(size: Int, initializer: (Int) -> T): Buffer = + bufferFactory(size, initializer) + + override fun check(vararg elements: NDBuffer) { + if (!elements.all { it.strides == this.strides }) error("Element strides are not the same as context strides") + } + + override val zero by lazy { produce { zero } } + override val one by lazy { produce { one } } + + override fun produce(initializer: F.(IntArray) -> T) = + BufferedNDFieldElement( + this, + buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }) + + override fun map(arg: NDBuffer, transform: F.(T) -> T): BufferedNDFieldElement { + check(arg) + return BufferedNDFieldElement( + this, + buildBuffer(arg.strides.linearSize) { offset -> elementContext.transform(arg.buffer[offset]) }) + +// val buffer = arg.buffer.transform { _, value -> elementContext.transform(value) } +// return BufferedNDFieldElement(this, buffer) + + } + + override fun mapIndexed( + arg: NDBuffer, + transform: F.(index: IntArray, T) -> T + ): BufferedNDFieldElement { + check(arg) + return BufferedNDFieldElement( + this, + buildBuffer(arg.strides.linearSize) { offset -> + elementContext.transform( + arg.strides.index(offset), + arg.buffer[offset] + ) + }) + +// val buffer = +// arg.buffer.transform { offset, value -> elementContext.transform(arg.strides.index(offset), value) } +// return BufferedNDFieldElement(this, buffer) + } + + override fun combine( + a: NDBuffer, + b: NDBuffer, + transform: F.(T, T) -> T + ): BufferedNDFieldElement { + check(a, b) + return BufferedNDFieldElement( + this, + buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) }) + } + + override fun NDBuffer.toElement(): FieldElement, *, out BufferedNDField> = + BufferedNDFieldElement(this@BoxingNDField, buffer) +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDAlgebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDAlgebra.kt new file mode 100644 index 000000000..9742f3662 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDAlgebra.kt @@ -0,0 +1,43 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.* + +interface BufferedNDAlgebra: NDAlgebra>{ + val strides: Strides + + override fun check(vararg elements: NDBuffer) { + if (!elements.all { it.strides == this.strides }) error("Strides mismatch") + } + + /** + * Convert any [NDStructure] to buffered structure using strides from this context. + * If the structure is already [NDBuffer], conversion is free. If not, it could be expensive because iteration over indexes + * + * If the argument is [NDBuffer] with different strides structure, the new element will be produced. + */ + fun NDStructure.toBuffer(): NDBuffer { + return if (this is NDBuffer && this.strides == this@BufferedNDAlgebra.strides) { + this + } else { + produce { index -> get(index) } + } + } + + /** + * Convert a buffer to element of this algebra + */ + fun NDBuffer.toElement(): MathElement> +} + + +interface BufferedNDSpace> : NDSpace>, BufferedNDAlgebra { + override fun NDBuffer.toElement(): SpaceElement, *, out BufferedNDSpace> +} + +interface BufferedNDRing> : NDRing>, BufferedNDSpace { + override fun NDBuffer.toElement(): RingElement, *, out BufferedNDRing> +} + +interface BufferedNDField> : NDField>, BufferedNDRing { + override fun NDBuffer.toElement(): FieldElement, *, out BufferedNDField> +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt new file mode 100644 index 000000000..04049368a --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferedNDElement.kt @@ -0,0 +1,88 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.* + +/** + * Base interface for an element with context, containing strides + */ +interface BufferedNDElement : NDBuffer, NDElement> { + override val context: BufferedNDAlgebra + + override val strides get() = context.strides + + override val shape: IntArray get() = context.shape +} + +class BufferedNDSpaceElement>( + override val context: BufferedNDSpace, + override val buffer: Buffer +) : BufferedNDElement, SpaceElement, BufferedNDSpaceElement, BufferedNDSpace> { + + override fun unwrap(): NDBuffer = this + + override fun NDBuffer.wrap(): BufferedNDSpaceElement { + context.check(this) + return BufferedNDSpaceElement(context, buffer) + } +} + +class BufferedNDRingElement>( + override val context: BufferedNDRing, + override val buffer: Buffer +) : BufferedNDElement, RingElement, BufferedNDRingElement, BufferedNDRing> { + + override fun unwrap(): NDBuffer = this + + override fun NDBuffer.wrap(): BufferedNDRingElement { + context.check(this) + return BufferedNDRingElement(context, buffer) + } +} + +class BufferedNDFieldElement>( + override val context: BufferedNDField, + override val buffer: Buffer +) : BufferedNDElement, FieldElement, BufferedNDFieldElement, BufferedNDField> { + + override fun unwrap(): NDBuffer = this + + override fun NDBuffer.wrap(): BufferedNDFieldElement { + context.check(this) + return BufferedNDFieldElement(context, buffer) + } +} + + +/** + * Element by element application of any operation on elements to the whole array. Just like in numpy + */ +operator fun > Function1.invoke(ndElement: BufferedNDElement) = + ndElement.context.run { map(ndElement) { invoke(it) }.toElement() } + +/* plus and minus */ + +/** + * Summation operation for [BufferedNDElement] and single element + */ +operator fun > BufferedNDElement.plus(arg: T) = + context.map(this) { it + arg }.wrap() + +/** + * Subtraction operation between [BufferedNDElement] and single element + */ +operator fun > BufferedNDElement.minus(arg: T) = + context.map(this) { it - arg }.wrap() + +/* prod and div */ + +/** + * Product operation for [BufferedNDElement] and single element + */ +operator fun > BufferedNDElement.times(arg: T) = + context.map(this) { it * arg }.wrap() + +/** + * Division operation between [BufferedNDElement] and single element + */ +operator fun > BufferedNDElement.div(arg: T) = + context.map(this) { it / arg }.wrap() \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt index 69e8163c6..a89729bae 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt @@ -1,20 +1,70 @@ package scientifik.kmath.structures +typealias BufferFactory = (Int, (Int) -> T) -> Buffer +typealias MutableBufferFactory = (Int, (Int) -> T) -> MutableBuffer + + /** * A generic random access structure for both primitives and objects */ interface Buffer { + /** + * The size of the buffer + */ val size: Int + /** + * Get element at given index + */ operator fun get(index: Int): T + /** + * Iterate over all elements + */ operator fun iterator(): Iterator + + /** + * Check content eqiality with another buffer + */ + fun contentEquals(other: Buffer<*>): Boolean = + asSequence().mapIndexed { index, value -> value == other[index] }.all { it } + + companion object { + + /** + * Create a boxing buffer of given type + */ + inline fun boxing(size: Int, initializer: (Int) -> T): Buffer = ListBuffer(List(size, initializer)) + + /** + * Create most appropriate immutable buffer for given type avoiding boxing wherever possible + */ + @Suppress("UNCHECKED_CAST") + inline fun auto(size: Int, crossinline initializer: (Int) -> T): Buffer { + return when (T::class) { + Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer + Short::class -> ShortBuffer(ShortArray(size) { initializer(it) as Short }) as Buffer + Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer + Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as Buffer + else -> boxing(size, initializer) + } + } + + val DoubleBufferFactory: BufferFactory = + { size, initializer -> DoubleBuffer(DoubleArray(size, initializer)) } + val ShortBufferFactory: BufferFactory = + { size, initializer -> ShortBuffer(ShortArray(size, initializer)) } + val IntBufferFactory: BufferFactory = { size, initializer -> IntBuffer(IntArray(size, initializer)) } + val LongBufferFactory: BufferFactory = { size, initializer -> LongBuffer(LongArray(size, initializer)) } + } } fun Buffer.asSequence(): Sequence = iterator().asSequence() +fun Buffer.asIterable(): Iterable = iterator().asSequence().asIterable() + interface MutableBuffer : Buffer { operator fun set(index: Int, value: T) @@ -22,10 +72,32 @@ interface MutableBuffer : Buffer { * A shallow copy of the buffer */ fun copy(): MutableBuffer + + companion object { + /** + * Create a boxing mutable buffer of given type + */ + inline fun boxing(size: Int, initializer: (Int) -> T): MutableBuffer = + MutableListBuffer(MutableList(size, initializer)) + + /** + * Create most appropriate mutable buffer for given type avoiding boxing wherever possible + */ + @Suppress("UNCHECKED_CAST") + inline fun auto(size: Int, initializer: (Int) -> T): MutableBuffer { + return when (T::class) { + Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer + Short::class -> ShortBuffer(ShortArray(size) { initializer(it) as Short }) as MutableBuffer + Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer + Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as MutableBuffer + else -> boxing(size, initializer) + } + } + } } -inline class ListBuffer(private val list: List) : Buffer { +inline class ListBuffer(val list: List) : Buffer { override val size: Int get() = list.size @@ -35,7 +107,12 @@ inline class ListBuffer(private val list: List) : Buffer { override fun iterator(): Iterator = list.iterator() } -inline class MutableListBuffer(private val list: MutableList) : MutableBuffer { +fun List.asBuffer() = ListBuffer(this) + +@Suppress("FunctionName") +inline fun ListBuffer(size: Int, init: (Int) -> T) = List(size, init).asBuffer() + +inline class MutableListBuffer(val list: MutableList) : MutableBuffer { override val size: Int get() = list.size @@ -47,12 +124,11 @@ inline class MutableListBuffer(private val list: MutableList) : MutableBuf } override fun iterator(): Iterator = list.iterator() - override fun copy(): MutableBuffer = MutableListBuffer(ArrayList(list)) } class ArrayBuffer(private val array: Array) : MutableBuffer { - //Can't inline because array invariant + //Can't inline because array is invariant override val size: Int get() = array.size @@ -67,9 +143,10 @@ class ArrayBuffer(private val array: Array) : MutableBuffer { override fun copy(): MutableBuffer = ArrayBuffer(array.copyOf()) } -inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer { - override val size: Int - get() = array.size +fun Array.asBuffer() = ArrayBuffer(this) + +inline class DoubleBuffer(val array: DoubleArray) : MutableBuffer { + override val size: Int get() = array.size override fun get(index: Int): Double = array[index] @@ -77,14 +154,45 @@ inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer = array.iterator() + override fun iterator() = array.iterator() override fun copy(): MutableBuffer = DoubleBuffer(array.copyOf()) } -inline class IntBuffer(private val array: IntArray) : MutableBuffer { - override val size: Int - get() = array.size +@Suppress("FunctionName") +inline fun DoubleBuffer(size: Int, init: (Int) -> Double) = DoubleBuffer(DoubleArray(size) { init(it) }) + +/** + * Transform buffer of doubles into array for high performance operations + */ +val Buffer.array: DoubleArray + get() = if (this is DoubleBuffer) { + array + } else { + DoubleArray(size) { get(it) } + } + +fun DoubleArray.asBuffer() = DoubleBuffer(this) + +inline class ShortBuffer(val array: ShortArray) : MutableBuffer { + override val size: Int get() = array.size + + override fun get(index: Int): Short = array[index] + + override fun set(index: Int, value: Short) { + array[index] = value + } + + override fun iterator() = array.iterator() + + override fun copy(): MutableBuffer = ShortBuffer(array.copyOf()) + +} + +fun ShortArray.asBuffer() = ShortBuffer(this) + +inline class IntBuffer(val array: IntArray) : MutableBuffer { + override val size: Int get() = array.size override fun get(index: Int): Int = array[index] @@ -92,17 +200,55 @@ inline class IntBuffer(private val array: IntArray) : MutableBuffer { array[index] = value } - override fun iterator(): Iterator = array.iterator() + override fun iterator() = array.iterator() override fun copy(): MutableBuffer = IntBuffer(array.copyOf()) + } -inline class ReadOnlyBuffer(private val buffer: MutableBuffer) : Buffer { +fun IntArray.asBuffer() = IntBuffer(this) + +inline class LongBuffer(val array: LongArray) : MutableBuffer { + override val size: Int get() = array.size + + override fun get(index: Int): Long = array[index] + + override fun set(index: Int, value: Long) { + array[index] = value + } + + override fun iterator() = array.iterator() + + override fun copy(): MutableBuffer = LongBuffer(array.copyOf()) + +} + +fun LongArray.asBuffer() = LongBuffer(this) + +inline class ReadOnlyBuffer(val buffer: MutableBuffer) : Buffer { override val size: Int get() = buffer.size override fun get(index: Int): T = buffer.get(index) - override fun iterator(): Iterator = buffer.iterator() + override fun iterator() = buffer.iterator() +} + +/** + * A buffer with content calculated on-demand. The calculated contect is not stored, so it is recalculated on each call. + * Useful when one needs single element from the buffer. + */ +class VirtualBuffer(override val size: Int, private val generator: (Int) -> T) : Buffer { + override fun get(index: Int): T = generator(index) + + override fun iterator(): Iterator = (0 until size).asSequence().map(generator).iterator() + + override fun contentEquals(other: Buffer<*>): Boolean { + return if (other is VirtualBuffer) { + this.size == other.size && this.generator == other.generator + } else { + super.contentEquals(other) + } + } } /** @@ -115,25 +261,6 @@ fun Buffer.asReadOnly(): Buffer = if (this is MutableBuffer) { } /** - * Create most appropriate immutable buffer for given type avoiding boxing wherever possible + * Typealias for buffer transformations */ -@Suppress("UNCHECKED_CAST") -inline fun buffer(size: Int, noinline initializer: (Int) -> T): Buffer { - return when (T::class) { - Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer - Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer - else -> ArrayBuffer(Array(size, initializer)) - } -} - -/** - * Create most appropriate mutable buffer for given type avoiding boxing wherever possible - */ -@Suppress("UNCHECKED_CAST") -inline fun mutableBuffer(size: Int, noinline initializer: (Int) -> T): MutableBuffer { - return when (T::class) { - Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer - Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer - else -> ArrayBuffer(Array(size, initializer)) - } -} \ No newline at end of file +typealias BufferTransform = (Buffer) -> Buffer \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ComplexNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ComplexNDField.kt new file mode 100644 index 000000000..6444f667d --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ComplexNDField.kt @@ -0,0 +1,134 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.Complex +import scientifik.kmath.operations.ComplexField +import scientifik.kmath.operations.FieldElement +import scientifik.kmath.operations.complex + +typealias ComplexNDElement = BufferedNDFieldElement + +/** + * An optimized nd-field for complex numbers + */ +class ComplexNDField(override val shape: IntArray) : + BufferedNDField, + ExtendedNDField> { + + override val strides: Strides = DefaultStrides(shape) + + override val elementContext: ComplexField get() = ComplexField + override val zero by lazy { produce { zero } } + override val one by lazy { produce { one } } + + inline fun buildBuffer(size: Int, crossinline initializer: (Int) -> Complex): Buffer = + Buffer.complex(size) { initializer(it) } + + /** + * Inline transform an NDStructure to another structure + */ + override fun map( + arg: NDBuffer, + transform: ComplexField.(Complex) -> Complex + ): ComplexNDElement { + check(arg) + val array = buildBuffer(arg.strides.linearSize) { offset -> ComplexField.transform(arg.buffer[offset]) } + return BufferedNDFieldElement(this, array) + } + + override fun produce(initializer: ComplexField.(IntArray) -> Complex): ComplexNDElement { + val array = buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) } + return BufferedNDFieldElement(this, array) + } + + override fun mapIndexed( + arg: NDBuffer, + transform: ComplexField.(index: IntArray, Complex) -> Complex + ): ComplexNDElement { + check(arg) + return BufferedNDFieldElement( + this, + buildBuffer(arg.strides.linearSize) { offset -> + elementContext.transform( + arg.strides.index(offset), + arg.buffer[offset] + ) + }) + } + + override fun combine( + a: NDBuffer, + b: NDBuffer, + transform: ComplexField.(Complex, Complex) -> Complex + ): ComplexNDElement { + check(a, b) + return BufferedNDFieldElement( + this, + buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) }) + } + + override fun NDBuffer.toElement(): FieldElement, *, out BufferedNDField> = + BufferedNDFieldElement(this@ComplexNDField, buffer) + + override fun power(arg: NDBuffer, pow: Number) = map(arg) { power(it, pow) } + + override fun exp(arg: NDBuffer) = map(arg) { exp(it) } + + override fun ln(arg: NDBuffer) = map(arg) { ln(it) } + + override fun sin(arg: NDBuffer) = map(arg) { sin(it) } + + override fun cos(arg: NDBuffer) = map(arg) { cos(it) } + +} + + +/** + * Fast element production using function inlining + */ +inline fun BufferedNDField.produceInline(crossinline initializer: ComplexField.(Int) -> Complex): ComplexNDElement { + val buffer = Buffer.complex(strides.linearSize) { offset -> ComplexField.initializer(offset) } + return BufferedNDFieldElement(this, buffer) +} + +/** + * Map one [ComplexNDElement] using function with indexes + */ +inline fun ComplexNDElement.mapIndexed(crossinline transform: ComplexField.(index: IntArray, Complex) -> Complex) = + context.produceInline { offset -> transform(strides.index(offset), buffer[offset]) } + +/** + * Map one [ComplexNDElement] using function without indexes + */ +inline fun ComplexNDElement.map(crossinline transform: ComplexField.(Complex) -> Complex): ComplexNDElement { + val buffer = Buffer.complex(strides.linearSize) { offset -> ComplexField.transform(buffer[offset]) } + return BufferedNDFieldElement(context, buffer) +} + +/** + * Element by element application of any operation on elements to the whole array. Just like in numpy + */ +operator fun Function1.invoke(ndElement: ComplexNDElement) = + ndElement.map { this@invoke(it) } + + +/* plus and minus */ + +/** + * Summation operation for [BufferedNDElement] and single element + */ +operator fun ComplexNDElement.plus(arg: Complex) = + map { it + arg } + +/** + * Subtraction operation between [BufferedNDElement] and single element + */ +operator fun ComplexNDElement.minus(arg: Complex) = + map { it - arg } + +operator fun ComplexNDElement.plus(arg: Double) = + map { it + arg } + +operator fun ComplexNDElement.minus(arg: Double) = + map { it - arg } + +fun NDField.Companion.complex(vararg shape: Int) = ComplexNDField(shape) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt index 93ffad5a3..370dc646e 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt @@ -1,42 +1,46 @@ package scientifik.kmath.structures -import scientifik.kmath.operations.ExponentialOperations -import scientifik.kmath.operations.ExtendedField -import scientifik.kmath.operations.PowerOperations -import scientifik.kmath.operations.TrigonometricOperations +import scientifik.kmath.operations.* -/** - * NDField that supports [ExtendedField] operations on its elements - */ -class ExtendedNDField>(shape: IntArray, field: F) : NDField(shape, field), - TrigonometricOperations>, - PowerOperations>, - ExponentialOperations> { - - override fun produceStructure(initializer: F.(IntArray) -> N): NDStructure { - return genericNdStructure(shape) { field.initializer(it) } - } - - override fun power(arg: NDElement, pow: Double): NDElement { - return arg.transform { d -> with(field) { power(d, pow) } } - } - - override fun exp(arg: NDElement): NDElement { - return arg.transform { d -> with(field) { exp(d) } } - } - - override fun ln(arg: NDElement): NDElement { - return arg.transform { d -> with(field) { ln(d) } } - } - - override fun sin(arg: NDElement): NDElement { - return arg.transform { d -> with(field) { sin(d) } } - } - - override fun cos(arg: NDElement): NDElement { - return arg.transform { d -> with(field) { cos(d) } } - } -} +interface ExtendedNDField> : + NDField, + TrigonometricOperations, + PowerOperations, + ExponentialOperations + where F : ExtendedFieldOperations, F : Field + + +///** +// * NDField that supports [ExtendedField] operations on its elements +// */ +//class ExtendedNDFieldWrapper, N : NDStructure>(private val ndField: NDField) : +// ExtendedNDField, NDField by ndField { +// +// override val shape: IntArray get() = ndField.shape +// override val elementContext: F get() = ndField.elementContext +// +// override fun produce(initializer: F.(IntArray) -> T) = ndField.produce(initializer) +// +// override fun power(arg: N, pow: Double): N { +// return produce { with(elementContext) { power(arg[it], pow) } } +// } +// +// override fun exp(arg: N): N { +// return produce { with(elementContext) { exp(arg[it]) } } +// } +// +// override fun ln(arg: N): N { +// return produce { with(elementContext) { ln(arg[it]) } } +// } +// +// override fun sin(arg: N): N { +// return produce { with(elementContext) { sin(arg[it]) } } +// } +// +// override fun cos(arg: N): N { +// return produce { with(elementContext) { cos(arg[it]) } } +// } +//} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt deleted file mode 100644 index 0428535f9..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt +++ /dev/null @@ -1,20 +0,0 @@ -package scientifik.kmath.structures - -// -//class LazyStructureField(val field: Field): Field>{ -// -//} -// -//class LazyStructure : NDStructure { -// -// override val shape: IntArray -// get() = TODO("not implemented") //To change initializer of created properties use File | Settings | File Templates. -// -// override fun get(index: IntArray): T { -// TODO("not implemented") //To change body of created functions use File | Settings | File Templates. -// } -// -// override fun iterator(): Iterator> { -// TODO("not implemented") //To change body of created functions use File | Settings | File Templates. -// } -//} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt new file mode 100644 index 000000000..887a480bb --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt @@ -0,0 +1,144 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.Field +import scientifik.kmath.operations.Ring +import scientifik.kmath.operations.Space +import kotlin.jvm.JvmName + + +/** + * An exception is thrown when the expected ans actual shape of NDArray differs + */ +class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : RuntimeException() + + +/** + * The base interface for all nd-algebra implementations + * @param T the type of nd-structure element + * @param C the type of the element context + * @param N the type of the structure + */ +interface NDAlgebra> { + val shape: IntArray + val elementContext: C + + /** + * Produce a new [N] structure using given initializer function + */ + fun produce(initializer: C.(IntArray) -> T): N + + /** + * Map elements from one structure to another one + */ + fun map(arg: N, transform: C.(T) -> T): N + + /** + * Map indexed elements + */ + fun mapIndexed(arg: N, transform: C.(index: IntArray, T) -> T): N + + /** + * Combine two structures into one + */ + fun combine(a: N, b: N, transform: C.(T, T) -> T): N + + /** + * Check if given elements are consistent with this context + */ + fun check(vararg elements: N) { + elements.forEach { + if (!shape.contentEquals(it.shape)) { + throw ShapeMismatchException(shape, it.shape) + } + } + } + + /** + * element-by-element invoke a function working on [T] on a [NDStructure] + */ + operator fun Function1.invoke(structure: N) = map(structure) { value -> this@invoke(value) } +} + +/** + * An nd-space over element space + */ +interface NDSpace, N : NDStructure> : Space, NDAlgebra { + /** + * Element-by-element addition + */ + override fun add(a: N, b: N): N = combine(a, b) { aValue, bValue -> add(aValue, bValue) } + + /** + * Multiply all elements by constant + */ + override fun multiply(a: N, k: Number): N = map(a) { multiply(it, k) } + + operator fun N.plus(arg: T) = map(this) { value -> add(arg, value) } + operator fun N.minus(arg: T) = map(this) { value -> add(arg, -value) } + + operator fun T.plus(arg: N) = map(arg) { value -> add(this@plus, value) } + operator fun T.minus(arg: N) = map(arg) { value -> add(-this@minus, value) } +} + +/** + * An nd-ring over element ring + */ +interface NDRing, N : NDStructure> : Ring, NDSpace { + + /** + * Element-by-element multiplication + */ + override fun multiply(a: N, b: N): N = combine(a, b) { aValue, bValue -> multiply(aValue, bValue) } + + operator fun N.times(arg: T) = map(this) { value -> multiply(arg, value) } + operator fun T.times(arg: N) = map(arg) { value -> multiply(this@times, value) } +} + +/** + * Field for n-dimensional structures. + * @param shape - the list of dimensions of the array + * @param elementField - operations field defined on individual array element + * @param T - the type of the element contained in ND structure + * @param F - field of structure elements + * @param R - actual nd-element type of this field + */ +interface NDField, N : NDStructure> : Field, NDRing { + + /** + * Element-by-element division + */ + override fun divide(a: N, b: N): N = combine(a, b) { aValue, bValue -> divide(aValue, bValue) } + + operator fun N.div(arg: T) = map(this) { value -> divide(arg, value) } + operator fun T.div(arg: N) = map(arg) { divide(it, this@div) } + + companion object { + + private val realNDFieldCache = HashMap() + + /** + * Create a nd-field for [Double] values or pull it from cache if it was created previously + */ + fun real(vararg shape: Int) = realNDFieldCache.getOrPut(shape) { RealNDField(shape) } + + /** + * Create a nd-field with boxing generic buffer + */ + fun > buffered( + shape: IntArray, + field: F, + bufferFactory: BufferFactory = Buffer.Companion::boxing + ) = + BoxingNDField(shape, field, bufferFactory) + + /** + * Create a most suitable implementation for nd-field using reified class. + */ + @Suppress("UNCHECKED_CAST") + inline fun > auto(field: F, vararg shape: Int): BufferedNDField = + when { + T::class == Double::class -> real(*shape) as BufferedNDField + else -> BoxingNDField(shape, field, Buffer.Companion::auto) + } + } +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt new file mode 100644 index 000000000..59cf98df3 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt @@ -0,0 +1,132 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.Field +import scientifik.kmath.operations.RealField +import scientifik.kmath.operations.Ring +import scientifik.kmath.operations.Space + +/** + * The root for all [NDStructure] based algebra elements. Does not implement algebra element root because of problems with recursive self-types + * @param T the type of the element of the structure + * @param C the type of the context for the element + * @param N the type of the underlying [NDStructure] + */ +interface NDElement> : NDStructure { + + val context: NDAlgebra + + fun unwrap(): N + + fun N.wrap(): NDElement + + companion object { + /** + * Create a optimized NDArray of doubles + */ + fun real(shape: IntArray, initializer: RealField.(IntArray) -> Double = { 0.0 }) = + NDField.real(*shape).produce(initializer) + + + fun real1D(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }) = + real(intArrayOf(dim)) { initializer(it[0]) } + + + fun real2D(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }) = + real(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) } + + fun real3D(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }) = + real(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) } + + + /** + * Simple boxing NDArray + */ + fun > buffered( + shape: IntArray, + field: F, + initializer: F.(IntArray) -> T + ): BufferedNDElement { + val ndField = BoxingNDField(shape, field, Buffer.Companion::boxing) + return ndField.produce(initializer) + } + + inline fun > auto( + shape: IntArray, + field: F, + noinline initializer: F.(IntArray) -> T + ): BufferedNDFieldElement { + val ndField = NDField.auto(field, *shape) + return BufferedNDFieldElement(ndField, ndField.produce(initializer).buffer) + } + } +} + + +fun > NDElement.mapIndexed(transform: C.(index: IntArray, T) -> T) = + context.mapIndexed(unwrap(), transform).wrap() + +fun > NDElement.map(transform: C.(T) -> T) = context.map(unwrap(), transform).wrap() + + +/** + * Element by element application of any operation on elements to the whole [NDElement] + */ +operator fun > Function1.invoke(ndElement: NDElement) = + ndElement.map { value -> this@invoke(value) } + +/* plus and minus */ + +/** + * Summation operation for [NDElement] and single element + */ +operator fun , N : NDStructure> NDElement.plus(arg: T) = + map { value -> arg + value } + +/** + * Subtraction operation between [NDElement] and single element + */ +operator fun , N : NDStructure> NDElement.minus(arg: T) = + map { value -> arg - value } + +/* prod and div */ + +/** + * Product operation for [NDElement] and single element + */ +operator fun , N : NDStructure> NDElement.times(arg: T) = + map { value -> arg * value } + +/** + * Division operation between [NDElement] and single element + */ +operator fun , N : NDStructure> NDElement.div(arg: T) = + map { value -> arg / value } + + +// /** +// * Reverse sum operation +// */ +// operator fun T.plus(arg: NDStructure): NDElement = produce { index -> +// field.run { this@plus + arg[index] } +// } +// +// /** +// * Reverse minus operation +// */ +// operator fun T.minus(arg: NDStructure): NDElement = produce { index -> +// field.run { this@minus - arg[index] } +// } +// +// /** +// * Reverse product operation +// */ +// operator fun T.times(arg: NDStructure): NDElement = produce { index -> +// field.run { this@times * arg[index] } +// } +// +// /** +// * Reverse division operation +// */ +// operator fun T.div(arg: NDStructure): NDElement = produce { index -> +// field.run { this@div / arg[index] } +// } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt deleted file mode 100644 index cddac5801..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt +++ /dev/null @@ -1,209 +0,0 @@ -package scientifik.kmath.structures - -import scientifik.kmath.operations.DoubleField -import scientifik.kmath.operations.Field -import scientifik.kmath.operations.FieldElement - -/** - * An exception is thrown when the expected ans actual shape of NDArray differs - */ -class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : RuntimeException() - -/** - * Field for n-dimensional arrays. - * @param shape - the list of dimensions of the array - * @param field - operations field defined on individual array element - * @param T the type of the element contained in NDArray - */ -abstract class NDField>(val shape: IntArray, val field: F) : Field> { - - abstract fun produceStructure(initializer: F.(IntArray) -> T): NDStructure - - /** - * Create new instance of NDArray using field shape and given initializer - * The producer takes list of indices as argument and returns contained value - */ - fun produce(initializer: F.(IntArray) -> T): NDElement = NDStructureElement(this, produceStructure(initializer)) - - override val zero: NDElement by lazy { - produce { zero } - } - - /** - * Check the shape of given NDArray and throw exception if it does not coincide with shape of the field - */ - private fun checkShape(vararg elements: NDElement) { - elements.forEach { - if (!shape.contentEquals(it.shape)) { - throw ShapeMismatchException(shape, it.shape) - } - } - } - - /** - * Element-by-element addition - */ - override fun add(a: NDElement, b: NDElement): NDElement { - checkShape(a, b) - return produce { with(field) { a[it] + b[it] } } - } - - /** - * Multiply all elements by cinstant - */ - override fun multiply(a: NDElement, k: Double): NDElement { - checkShape(a) - return produce { with(field) { a[it] * k } } - } - - override val one: NDElement - get() = produce { one } - - /** - * Element-by-element multiplication - */ - override fun multiply(a: NDElement, b: NDElement): NDElement { - checkShape(a) - return produce { with(field) { a[it] * b[it] } } - } - - /** - * Element-by-element division - */ - override fun divide(a: NDElement, b: NDElement): NDElement { - checkShape(a) - return produce { with(field) { a[it] / b[it] } } - } - -// /** -// * Reverse sum operation -// */ -// operator fun T.plus(arg: NDElement): NDElement = arg + this -// -// /** -// * Reverse minus operation -// */ -// operator fun T.minus(arg: NDElement): NDElement = arg.transformIndexed { _, value -> -// with(arg.context.field) { -// this@minus - value -// } -// } -// -// /** -// * Reverse product operation -// */ -// operator fun T.times(arg: NDElement): NDElement = arg * this -// -// /** -// * Reverse division operation -// */ -// operator fun T.div(arg: NDElement): NDElement = arg.transformIndexed { _, value -> -// with(arg.context.field) { -// this@div / value -// } -// } -} - - -interface NDElement>: FieldElement, NDField>, NDStructure - -inline fun > NDElement.transformIndexed(crossinline action: F.(IntArray, T) -> T): NDElement = context.produce { action(it, get(*it)) } -inline fun > NDElement.transform(crossinline action: F.(T) -> T): NDElement = context.produce { action(get(*it)) } - - -/** - * Read-only [NDStructure] coupled to the context. - */ -class NDStructureElement>(override val context: NDField, private val structure: NDStructure) : NDElement, NDStructure by structure { - - //TODO ensure structure is immutable - - override val self: NDElement get() = this -} - -/** - * Element by element application of any operation on elements to the whole array. Just like in numpy - */ -operator fun > Function1.invoke(ndElement: NDElement): NDElement = ndElement.transform {value -> this@invoke(value) } - -/* plus and minus */ - -/** - * Summation operation for [NDElement] and single element - */ -operator fun > NDElement.plus(arg: T): NDElement = transform {value -> - with(context.field) { - arg + value - } -} - -/** - * Subtraction operation between [NDElement] and single element - */ -operator fun > NDElement.minus(arg: T): NDElement = transform {value -> - with(context.field) { - arg - value - } -} - -/* prod and div */ - -/** - * Product operation for [NDElement] and single element - */ -operator fun > NDElement.times(arg: T): NDElement = transform { value -> - with(context.field) { - arg * value - } -} - -/** - * Division operation between [NDElement] and single element - */ -operator fun > NDElement.div(arg: T): NDElement = transform { value -> - with(context.field) { - arg / value - } -} - -class GenericNDField>(shape: IntArray, field: F) : NDField(shape, field) { - override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure = genericNdStructure(shape) { field.initializer(it) } -} - -//typealias NDFieldFactory = (IntArray)->NDField - -object NDArrays { - /** - * Create a platform-optimized NDArray of doubles - */ - fun realNDArray(shape: IntArray, initializer: DoubleField.(IntArray) -> Double = { 0.0 }): NDElement { - return ExtendedNDField(shape, DoubleField).produce(initializer) - } - - fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement { - return realNDArray(intArrayOf(dim)) { initializer(it[0]) } - } - - fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement { - return realNDArray(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) } - } - - fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDElement { - return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) } - } - - inline fun produceReal(shape: IntArray, block: ExtendedNDField.() -> NDElement) = - ExtendedNDField(shape, DoubleField).run(block) - -// /** -// * Simple boxing NDField -// */ -// fun fieldFactory(field: Field): NDFieldFactory = { shape -> GenericNDField(shape, field) } - - /** - * Simple boxing NDArray - */ - fun > create(field: F, shape: IntArray, initializer: (IntArray) -> T): NDElement { - return GenericNDField(shape, field).produce { initializer(it) } - } -} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt index f9521cf93..2cf9af6fb 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt @@ -11,6 +11,45 @@ interface NDStructure { operator fun get(index: IntArray): T fun elements(): Sequence> + + companion object { + fun equals(st1: NDStructure<*>, st2: NDStructure<*>): Boolean { + return when { + st1 === st2 -> true + st1 is BufferNDStructure<*> && st2 is BufferNDStructure<*> && st1.strides == st2.strides -> st1.buffer.contentEquals( + st2.buffer + ) + else -> st1.elements().all { (index, value) -> value == st2[index] } + } + } + + /** + * Create a NDStructure with explicit buffer factory + * + * Strides should be reused if possible + */ + fun build( + strides: Strides, + bufferFactory: BufferFactory = Buffer.Companion::boxing, + initializer: (IntArray) -> T + ) = + BufferNDStructure(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) }) + + /** + * Inline create NDStructure with non-boxing buffer implementation if it is possible + */ + inline fun auto(strides: Strides, crossinline initializer: (IntArray) -> T) = + BufferNDStructure(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) }) + + fun build( + shape: IntArray, + bufferFactory: BufferFactory = Buffer.Companion::boxing, + initializer: (IntArray) -> T + ) = build(DefaultStrides(shape), bufferFactory, initializer) + + inline fun auto(shape: IntArray, crossinline initializer: (IntArray) -> T) = + auto(DefaultStrides(shape), initializer) + } } operator fun NDStructure.get(vararg index: Int): T = get(index) @@ -19,7 +58,7 @@ interface MutableNDStructure : NDStructure { operator fun set(index: IntArray, value: T) } -fun MutableNDStructure.transformInPlace(action: (IntArray, T) -> T) { +fun MutableNDStructure.mapInPlace(action: (IntArray, T) -> T) { elements().forEach { (index, oldValue) -> this[index] = action(index, oldValue) } @@ -49,6 +88,9 @@ interface Strides { */ fun index(offset: Int): IntArray + /** + * The size of linear buffer to accommodate all elements of ND-structure corresponding to strides + */ val linearSize: Int /** @@ -60,7 +102,7 @@ interface Strides { } } -class DefaultStrides(override val shape: IntArray) : Strides { +class DefaultStrides private constructor(override val shape: IntArray) : Strides { /** * Strides for memory access */ @@ -77,7 +119,7 @@ class DefaultStrides(override val shape: IntArray) : Strides { override fun offset(index: IntArray): Int { return index.mapIndexed { i, value -> - if (value < 0 || value >= shape[i]) { + if (value < 0 || value >= this.shape[i]) { throw RuntimeException("Index $value out of shape bounds: (0,${this.shape[i]})") } value * strides[i] @@ -98,50 +140,100 @@ class DefaultStrides(override val shape: IntArray) : Strides { override val linearSize: Int get() = strides[shape.size] + + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other !is DefaultStrides) return false + + if (!shape.contentEquals(other.shape)) return false + + return true + } + + override fun hashCode(): Int { + return shape.contentHashCode() + } + + companion object { + private val defaultStridesCache = HashMap() + + /** + * Cached builder for default strides + */ + operator fun invoke(shape: IntArray): Strides = defaultStridesCache.getOrPut(shape) { DefaultStrides(shape) } + } } -abstract class GenericNDStructure> : NDStructure { - protected abstract val buffer: B - protected abstract val strides: Strides +interface NDBuffer : NDStructure { + val buffer: Buffer + val strides: Strides override fun get(index: IntArray): T = buffer[strides.offset(index)] - override val shape: IntArray - get() = strides.shape + override val shape: IntArray get() = strides.shape - override fun elements()= - strides.indices().map { it to this[it] } + override fun elements() = strides.indices().map { it to this[it] } } /** * Boxing generic [NDStructure] */ class BufferNDStructure( - override val strides: Strides, - override val buffer: Buffer -) : GenericNDStructure>() { + override val strides: Strides, + override val buffer: Buffer +) : NDBuffer { init { if (strides.linearSize != buffer.size) { error("Expected buffer side of ${strides.linearSize}, but found ${buffer.size}") } } + + override fun get(index: IntArray): T = buffer[strides.offset(index)] + + override val shape: IntArray get() = strides.shape + + override fun elements() = strides.indices().map { it to this[it] } + + override fun equals(other: Any?): Boolean { + return when { + this === other -> true + other is BufferNDStructure<*> && this.strides == other.strides -> this.buffer.contentEquals(other.buffer) + other is NDStructure<*> -> elements().all { (index, value) -> value == other[index] } + else -> false + } + } + + override fun hashCode(): Int { + var result = strides.hashCode() + result = 31 * result + buffer.hashCode() + return result + } } -inline fun ndStructure(strides: Strides, noinline initializer: (IntArray) -> T) = - BufferNDStructure(strides, buffer(strides.linearSize) { i -> initializer(strides.index(i)) }) - -inline fun ndStructure(shape: IntArray, noinline initializer: (IntArray) -> T) = - ndStructure(DefaultStrides(shape), initializer) - +/** + * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferNDStructure] + */ +inline fun NDStructure.mapToBuffer( + factory: BufferFactory = Buffer.Companion::auto, + crossinline transform: (T) -> R +): BufferNDStructure { + return if (this is BufferNDStructure) { + BufferNDStructure(this.strides, factory.invoke(strides.linearSize) { transform(buffer[it]) }) + } else { + val strides = DefaultStrides(shape) + BufferNDStructure(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) }) + } +} /** - * Mutable ND buffer based on linear [Buffer] + * Mutable ND buffer based on linear [autoBuffer] */ class MutableBufferNDStructure( - override val strides: Strides, - override val buffer: MutableBuffer -) : GenericNDStructure>(), MutableNDStructure { + override val strides: Strides, + override val buffer: MutableBuffer +) : NDBuffer, MutableNDStructure { init { if (strides.linearSize != buffer.size) { @@ -152,25 +244,10 @@ class MutableBufferNDStructure( override fun set(index: IntArray, value: T) = buffer.set(strides.offset(index), value) } -/** - * Create optimized mutable structure for given type - */ -inline fun mutableNdStructure(strides: Strides, noinline initializer: (IntArray) -> T) = - MutableBufferNDStructure(strides, mutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) }) - -inline fun mutableNdStructure(shape: IntArray, noinline initializer: (IntArray) -> T) = - mutableNdStructure(DefaultStrides(shape), initializer) - -/** - * Create universal mutable structure - */ -fun genericNdStructure(shape: IntArray, initializer: (IntArray) -> T): MutableBufferNDStructure { - val strides = DefaultStrides(shape) - val sequence = sequence { - strides.indices().forEach { - yield(initializer(it)) - } - } - val buffer = MutableListBuffer(sequence.toMutableList()) - return MutableBufferNDStructure(strides, buffer) -} +inline fun NDStructure.combine( + struct: NDStructure, + crossinline block: (T, T) -> T +): NDStructure { + if (!this.shape.contentEquals(struct.shape)) error("Shape mismatch in structure combination") + return NDStructure.auto(shape) { block(this[it], struct[it]) } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt new file mode 100644 index 000000000..badf8a483 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt @@ -0,0 +1,48 @@ +package scientifik.kmath.structures + +import scientifik.memory.* + +/** + * A non-boxing buffer based on [ByteBuffer] storage + */ +open class ObjectBuffer(protected val memory: Memory, protected val spec: MemorySpec) : Buffer { + + override val size: Int get() = memory.size / spec.objectSize + + private val reader = memory.reader() + + override fun get(index: Int): T = reader.read(spec, spec.objectSize * index) + + override fun iterator(): Iterator = (0 until size).asSequence().map { get(it) }.iterator() + + + companion object { + fun create(spec: MemorySpec, size: Int) = + ObjectBuffer(Memory.allocate(size * spec.objectSize), spec) + + inline fun create( + spec: MemorySpec, + size: Int, + crossinline initializer: (Int) -> T + ): ObjectBuffer = + MutableObjectBuffer(Memory.allocate(size * spec.objectSize), spec).also { buffer -> + (0 until size).forEach { + buffer[it] = initializer(it) + } + } + } +} + +class MutableObjectBuffer(memory: Memory, spec: MemorySpec) : ObjectBuffer(memory, spec), + MutableBuffer { + + private val writer = memory.writer() + + override fun set(index: Int, value: T) = writer.write(spec, spec.objectSize * index, value) + + override fun copy(): MutableBuffer = MutableObjectBuffer(memory.copy(), spec) + + companion object { + + } +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealBufferField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealBufferField.kt new file mode 100644 index 000000000..b56f42717 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealBufferField.kt @@ -0,0 +1,153 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.ExtendedFieldOperations +import scientifik.kmath.operations.Field +import kotlin.math.* + + +/** + * A simple field over linear buffers of [Double] + */ +object RealBufferFieldOperations : ExtendedFieldOperations> { + override fun add(a: Buffer, b: Buffer): DoubleBuffer { + require(b.size == a.size) { "The size of the first buffer ${a.size} should be the same as for second one: ${b.size} " } + return if (a is DoubleBuffer && b is DoubleBuffer) { + val aArray = a.array + val bArray = b.array + DoubleBuffer(DoubleArray(a.size) { aArray[it] + bArray[it] }) + } else { + DoubleBuffer(DoubleArray(a.size) { a[it] + b[it] }) + } + } + + override fun multiply(a: Buffer, k: Number): DoubleBuffer { + val kValue = k.toDouble() + return if (a is DoubleBuffer) { + val aArray = a.array + DoubleBuffer(DoubleArray(a.size) { aArray[it] * kValue }) + } else { + DoubleBuffer(DoubleArray(a.size) { a[it] * kValue }) + } + } + + override fun multiply(a: Buffer, b: Buffer): DoubleBuffer { + require(b.size == a.size) { "The size of the first buffer ${a.size} should be the same as for second one: ${b.size} " } + return if (a is DoubleBuffer && b is DoubleBuffer) { + val aArray = a.array + val bArray = b.array + DoubleBuffer(DoubleArray(a.size) { aArray[it] * bArray[it] }) + } else { + DoubleBuffer(DoubleArray(a.size) { a[it] * b[it] }) + } + } + + override fun divide(a: Buffer, b: Buffer): DoubleBuffer { + require(b.size == a.size) { "The size of the first buffer ${a.size} should be the same as for second one: ${b.size} " } + return if (a is DoubleBuffer && b is DoubleBuffer) { + val aArray = a.array + val bArray = b.array + DoubleBuffer(DoubleArray(a.size) { aArray[it] / bArray[it] }) + } else { + DoubleBuffer(DoubleArray(a.size) { a[it] / b[it] }) + } + } + + override fun sin(arg: Buffer): DoubleBuffer { + return if (arg is DoubleBuffer) { + val array = arg.array + DoubleBuffer(DoubleArray(arg.size) { sin(array[it]) }) + } else { + DoubleBuffer(DoubleArray(arg.size) { sin(arg[it]) }) + } + } + + override fun cos(arg: Buffer): DoubleBuffer { + return if (arg is DoubleBuffer) { + val array = arg.array + DoubleBuffer(DoubleArray(arg.size) { cos(array[it]) }) + } else { + DoubleBuffer(DoubleArray(arg.size) { cos(arg[it]) }) + } + } + + override fun power(arg: Buffer, pow: Number): DoubleBuffer { + return if (arg is DoubleBuffer) { + val array = arg.array + DoubleBuffer(DoubleArray(arg.size) { array[it].pow(pow.toDouble()) }) + } else { + DoubleBuffer(DoubleArray(arg.size) { arg[it].pow(pow.toDouble()) }) + } + } + + override fun exp(arg: Buffer): DoubleBuffer { + return if (arg is DoubleBuffer) { + val array = arg.array + DoubleBuffer(DoubleArray(arg.size) { exp(array[it]) }) + } else { + DoubleBuffer(DoubleArray(arg.size) { exp(arg[it]) }) + } + } + + override fun ln(arg: Buffer): DoubleBuffer { + return if (arg is DoubleBuffer) { + val array = arg.array + DoubleBuffer(DoubleArray(arg.size) { ln(array[it]) }) + } else { + DoubleBuffer(DoubleArray(arg.size) { ln(arg[it]) }) + } + } +} + +class RealBufferField(val size: Int) : Field>, ExtendedFieldOperations> { + + override val zero: Buffer by lazy { DoubleBuffer(size) { 0.0 } } + + override val one: Buffer by lazy { DoubleBuffer(size) { 1.0 } } + + override fun add(a: Buffer, b: Buffer): DoubleBuffer { + require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } + return RealBufferFieldOperations.add(a, b) + } + + override fun multiply(a: Buffer, k: Number): DoubleBuffer { + require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } + return RealBufferFieldOperations.multiply(a, k) + } + + override fun multiply(a: Buffer, b: Buffer): DoubleBuffer { + require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } + return RealBufferFieldOperations.multiply(a, b) + } + + + override fun divide(a: Buffer, b: Buffer): DoubleBuffer { + require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } + return RealBufferFieldOperations.divide(a, b) + } + + override fun sin(arg: Buffer): DoubleBuffer { + require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } + return RealBufferFieldOperations.sin(arg) + } + + override fun cos(arg: Buffer): DoubleBuffer { + require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } + return RealBufferFieldOperations.cos(arg) + } + + override fun power(arg: Buffer, pow: Number): DoubleBuffer { + require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } + return RealBufferFieldOperations.power(arg, pow) + } + + override fun exp(arg: Buffer): DoubleBuffer { + require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } + return RealBufferFieldOperations.exp(arg) + } + + override fun ln(arg: Buffer): DoubleBuffer { + require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } + return RealBufferFieldOperations.ln(arg) + } + +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt new file mode 100644 index 000000000..82a237817 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt @@ -0,0 +1,121 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.FieldElement +import scientifik.kmath.operations.RealField + +typealias RealNDElement = BufferedNDFieldElement + +class RealNDField(override val shape: IntArray) : + BufferedNDField, + ExtendedNDField> { + + override val strides: Strides = DefaultStrides(shape) + + override val elementContext: RealField get() = RealField + override val zero by lazy { produce { zero } } + override val one by lazy { produce { one } } + + inline fun buildBuffer(size: Int, crossinline initializer: (Int) -> Double): Buffer = + DoubleBuffer(DoubleArray(size) { initializer(it) }) + + /** + * Inline transform an NDStructure to + */ + override fun map( + arg: NDBuffer, + transform: RealField.(Double) -> Double + ): RealNDElement { + check(arg) + val array = buildBuffer(arg.strides.linearSize) { offset -> RealField.transform(arg.buffer[offset]) } + return BufferedNDFieldElement(this, array) + } + + override fun produce(initializer: RealField.(IntArray) -> Double): RealNDElement { + val array = buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) } + return BufferedNDFieldElement(this, array) + } + + override fun mapIndexed( + arg: NDBuffer, + transform: RealField.(index: IntArray, Double) -> Double + ): RealNDElement { + check(arg) + return BufferedNDFieldElement( + this, + buildBuffer(arg.strides.linearSize) { offset -> + elementContext.transform( + arg.strides.index(offset), + arg.buffer[offset] + ) + }) + } + + override fun combine( + a: NDBuffer, + b: NDBuffer, + transform: RealField.(Double, Double) -> Double + ): RealNDElement { + check(a, b) + return BufferedNDFieldElement( + this, + buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) }) + } + + override fun NDBuffer.toElement(): FieldElement, *, out BufferedNDField> = + BufferedNDFieldElement(this@RealNDField, buffer) + + override fun power(arg: NDBuffer, pow: Number) = map(arg) { power(it, pow) } + + override fun exp(arg: NDBuffer) = map(arg) { exp(it) } + + override fun ln(arg: NDBuffer) = map(arg) { ln(it) } + + override fun sin(arg: NDBuffer) = map(arg) { sin(it) } + + override fun cos(arg: NDBuffer) = map(arg) { cos(it) } + +} + + +/** + * Fast element production using function inlining + */ +inline fun BufferedNDField.produceInline(crossinline initializer: RealField.(Int) -> Double): RealNDElement { + val array = DoubleArray(strides.linearSize) { offset -> RealField.initializer(offset) } + return BufferedNDFieldElement(this, DoubleBuffer(array)) +} + +/** + * Map one [RealNDElement] using function with indexes + */ +inline fun RealNDElement.mapIndexed(crossinline transform: RealField.(index: IntArray, Double) -> Double) = + context.produceInline { offset -> transform(strides.index(offset), buffer[offset]) } + +/** + * Map one [RealNDElement] using function without indexes + */ +inline fun RealNDElement.map(crossinline transform: RealField.(Double) -> Double): RealNDElement { + val array = DoubleArray(strides.linearSize) { offset -> RealField.transform(buffer[offset]) } + return BufferedNDFieldElement(context, DoubleBuffer(array)) +} + +/** + * Element by element application of any operation on elements to the whole array. Just like in numpy + */ +operator fun Function1.invoke(ndElement: RealNDElement) = + ndElement.map { this@invoke(it) } + + +/* plus and minus */ + +/** + * Summation operation for [BufferedNDElement] and single element + */ +operator fun RealNDElement.plus(arg: Double) = + map { it + arg } + +/** + * Subtraction operation between [BufferedNDElement] and single element + */ +operator fun RealNDElement.minus(arg: Double) = + map { it - arg } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt new file mode 100644 index 000000000..6b09c91de --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ShortNDRing.kt @@ -0,0 +1,96 @@ +package scientifik.kmath.structures + +import scientifik.kmath.operations.RingElement +import scientifik.kmath.operations.ShortRing + + +typealias ShortNDElement = BufferedNDRingElement + +class ShortNDRing(override val shape: IntArray) : + BufferedNDRing { + + override val strides: Strides = DefaultStrides(shape) + + override val elementContext: ShortRing get() = ShortRing + override val zero by lazy { produce { ShortRing.zero } } + override val one by lazy { produce { ShortRing.one } } + + inline fun buildBuffer(size: Int, crossinline initializer: (Int) -> Short): Buffer = + ShortBuffer(ShortArray(size) { initializer(it) }) + + /** + * Inline transform an NDStructure to + */ + override fun map( + arg: NDBuffer, + transform: ShortRing.(Short) -> Short + ): ShortNDElement { + check(arg) + val array = buildBuffer(arg.strides.linearSize) { offset -> ShortRing.transform(arg.buffer[offset]) } + return BufferedNDRingElement(this, array) + } + + override fun produce(initializer: ShortRing.(IntArray) -> Short): ShortNDElement { + val array = buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) } + return BufferedNDRingElement(this, array) + } + + override fun mapIndexed( + arg: NDBuffer, + transform: ShortRing.(index: IntArray, Short) -> Short + ): ShortNDElement { + check(arg) + return BufferedNDRingElement( + this, + buildBuffer(arg.strides.linearSize) { offset -> + elementContext.transform( + arg.strides.index(offset), + arg.buffer[offset] + ) + }) + } + + override fun combine( + a: NDBuffer, + b: NDBuffer, + transform: ShortRing.(Short, Short) -> Short + ): ShortNDElement { + check(a, b) + return BufferedNDRingElement( + this, + buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) }) + } + + override fun NDBuffer.toElement(): RingElement, *, out BufferedNDRing> = + BufferedNDRingElement(this@ShortNDRing, buffer) +} + + +/** + * Fast element production using function inlining + */ +inline fun BufferedNDRing.produceInline(crossinline initializer: ShortRing.(Int) -> Short): ShortNDElement { + val array = ShortArray(strides.linearSize) { offset -> ShortRing.initializer(offset) } + return BufferedNDRingElement(this, ShortBuffer(array)) +} + +/** + * Element by element application of any operation on elements to the whole array. Just like in numpy + */ +operator fun Function1.invoke(ndElement: ShortNDElement) = + ndElement.context.produceInline { i -> invoke(ndElement.buffer[i]) } + + +/* plus and minus */ + +/** + * Summation operation for [StridedNDFieldElement] and single element + */ +operator fun ShortNDElement.plus(arg: Short) = + context.produceInline { i -> (buffer[i] + arg).toShort() } + +/** + * Subtraction operation between [StridedNDFieldElement] and single element + */ +operator fun ShortNDElement.minus(arg: Short) = + context.produceInline { i -> (buffer[i] - arg).toShort() } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/SpecializedStructures.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/SpecializedStructures.kt new file mode 100644 index 000000000..734aba5ac --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/SpecializedStructures.kt @@ -0,0 +1,94 @@ +package scientifik.kmath.structures + +/** + * A structure that is guaranteed to be one-dimensional + */ +interface Structure1D : NDStructure, Buffer { + override val dimension: Int get() = 1 + + override fun get(index: IntArray): T { + if (index.size != 1) error("Index dimension mismatch. Expected 1 but found ${index.size}") + return get(index[0]) + } + + override fun iterator(): Iterator = (0 until size).asSequence().map { get(it) }.iterator() +} + +/** + * A 1D wrapper for nd-structure + */ +private inline class Structure1DWrapper(val structure: NDStructure) : Structure1D { + + override val shape: IntArray get() = structure.shape + override val size: Int get() = structure.shape[0] + + override fun get(index: Int): T = structure[index] + + override fun elements(): Sequence> = structure.elements() +} + +/** + * Represent a [NDStructure] as [Structure1D]. Throw error in case of dimension mismatch + */ +fun NDStructure.as1D(): Structure1D = if (shape.size == 1) { + Structure1DWrapper(this) +} else { + error("Can't create 1d-structure from ${shape.size}d-structure") +} + +fun NDBuffer.as1D(): Structure1D = if (shape.size == 1) { + Buffer1DWrapper(this.buffer) +} else { + error("Can't create 1d-structure from ${shape.size}d-structure") +} + +/** + * A structure wrapper for buffer + */ +private inline class Buffer1DWrapper(val buffer: Buffer) : Structure1D { + override val shape: IntArray get() = intArrayOf(buffer.size) + + override val size: Int get() = buffer.size + + override fun elements(): Sequence> = + asSequence().mapIndexed { index, value -> intArrayOf(index) to value } + + override fun get(index: Int): T = buffer.get(index) +} + +/** + * Represent this buffer as 1D structure + */ +fun Buffer.asND(): Structure1D = Buffer1DWrapper(this) + +/** + * A structure that is guaranteed to be two-dimensional + */ +interface Structure2D : NDStructure { + operator fun get(i: Int, j: Int): T + + override fun get(index: IntArray): T { + if (index.size != 2) error("Index dimension mismatch. Expected 2 but found ${index.size}") + return get(index[0], index[1]) + } +} + +/** + * A 2D wrapper for nd-structure + */ +private inline class Structure2DWrapper(val structure: NDStructure) : Structure2D { + override fun get(i: Int, j: Int): T = structure[i, j] + + override val shape: IntArray get() = structure.shape + + override fun elements(): Sequence> = structure.elements() +} + +/** + * Represent a [NDStructure] as [Structure1D]. Throw error in case of dimension mismatch + */ +fun NDStructure.as2D(): Structure2D = if (shape.size == 2) { + Structure2DWrapper(this) +} else { + error("Can't create 2d-structure from ${shape.size}d-structure") +} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt similarity index 59% rename from kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt index 8e2845b9e..033b2792f 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt @@ -2,17 +2,17 @@ package scientifik.kmath.expressions import scientifik.kmath.operations.Complex import scientifik.kmath.operations.ComplexField -import scientifik.kmath.operations.DoubleField +import scientifik.kmath.operations.RealField import kotlin.test.Test import kotlin.test.assertEquals -class FieldExpressionContextTest { +class ExpressionFieldTest { @Test fun testExpression() { - val context = FieldExpressionContext(DoubleField) + val context = ExpressionField(RealField) val expression = with(context) { val x = variable("x", 2.0) - x * x + 2 * x + 1.0 + x * x + 2 * x + one } assertEquals(expression("x" to 1.0), 4.0) assertEquals(expression(), 9.0) @@ -20,10 +20,10 @@ class FieldExpressionContextTest { @Test fun testComplex() { - val context = FieldExpressionContext(ComplexField) + val context = ExpressionField(ComplexField) val expression = with(context) { val x = variable("x", Complex(2.0, 0.0)) - x * x + 2 * x + 1.0 + x * x + 2 * x + one } assertEquals(expression("x" to Complex(1.0, 0.0)), Complex(4.0, 0.0)) assertEquals(expression(), Complex(9.0, 0.0)) @@ -31,23 +31,23 @@ class FieldExpressionContextTest { @Test fun separateContext() { - fun FieldExpressionContext.expression(): Expression{ + fun ExpressionField.expression(): Expression { val x = variable("x") - return x * x + 2 * x + 1.0 + return x * x + 2 * x + one } - val expression = FieldExpressionContext(DoubleField).expression() + val expression = ExpressionField(RealField).expression() assertEquals(expression("x" to 1.0), 4.0) } @Test fun valueExpression() { - val expressionBuilder: FieldExpressionContext.()->Expression = { + val expressionBuilder: ExpressionField.() -> Expression = { val x = variable("x") - x * x + 2 * x + 1.0 + x * x + 2 * x + one } - val expression = FieldExpressionContext(DoubleField).expressionBuilder() + val expression = ExpressionField(RealField).expressionBuilder() assertEquals(expression("x" to 1.0), 4.0) } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt deleted file mode 100644 index f767d682b..000000000 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt +++ /dev/null @@ -1,34 +0,0 @@ -package scientifik.kmath.linear - -import kotlin.test.Test -import kotlin.test.assertEquals - -class ArrayMatrixTest { - - @Test - fun testSum() { - val vector1 = Vector.ofReal(5) { it.toDouble() } - val vector2 = Vector.ofReal(5) { 5 - it.toDouble() } - val sum = vector1 + vector2 - assertEquals(5.0, sum[2]) - } - - @Test - fun testVectorToMatrix() { - val vector = Vector.ofReal(5) { it.toDouble() } - val matrix = vector.toMatrix() - assertEquals(4.0, matrix[4, 0]) - } - - - @Test - fun testDot() { - val vector1 = Vector.ofReal(5) { it.toDouble() } - val vector2 = Vector.ofReal(5) { 5 - it.toDouble() } - val product = vector1.toMatrix() dot (vector2.toMatrix().transpose()) - - - assertEquals(5.0, product[1, 0]) - assertEquals(6.0, product[2, 2]) - } -} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt new file mode 100644 index 000000000..61aa506c4 --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt @@ -0,0 +1,54 @@ +package scientifik.kmath.linear + +import kotlin.test.Test +import kotlin.test.assertEquals + +class MatrixTest { + + @Test + fun testSum() { + val vector1 = Vector.real(5) { it.toDouble() } + val vector2 = Vector.real(5) { 5 - it.toDouble() } + val sum = vector1 + vector2 + assertEquals(5.0, sum[2]) + } + + @Test + fun testVectorToMatrix() { + val vector = Vector.real(5) { it.toDouble() } + val matrix = vector.toMatrix() + assertEquals(4.0, matrix[4, 0]) + } + + @Test + fun testTranspose() { + val matrix = MatrixContext.real.one(3, 3) + val transposed = matrix.transpose() + assertEquals(matrix, transposed) + } + + + @Test + fun testDot() { + val vector1 = Vector.real(5) { it.toDouble() } + val vector2 = Vector.real(5) { 5 - it.toDouble() } + + val matrix1 = vector1.toMatrix() + val matrix2 = vector2.toMatrix().transpose() + val product = MatrixContext.real.run { matrix1 dot matrix2 } + + + assertEquals(5.0, product[1, 0]) + assertEquals(6.0, product[2, 2]) + } + + @Test + fun testBuilder() { + val matrix = Matrix.build(2, 3)( + 1.0, 0.0, 0.0, + 0.0, 1.0, 2.0 + ) + + assertEquals(2.0, matrix[1, 2]) + } +} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt index 472d5262f..bfa720369 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt @@ -1,21 +1,41 @@ package scientifik.kmath.linear -import scientifik.kmath.operations.DoubleField import kotlin.test.Test -import kotlin.test.assertTrue +import kotlin.test.assertEquals class RealLUSolverTest { @Test fun testInvertOne() { - val matrix = Matrix.diagonal(2, 2, DoubleField) - val inverted = RealLUSolver.inverse(matrix) - assertTrue { Matrix.equals(matrix,inverted) } + val matrix = MatrixContext.real.one(2, 2) + val inverted = LUSolver.real.inverse(matrix) + assertEquals(matrix, inverted) } -// @Test -// fun testInvert() { -// val matrix = realMatrix(2,2){} -// val inverted = RealLUSolver.inverse(matrix) -// assertTrue { Matrix.equals(matrix,inverted) } -// } + @Test + fun testInvert() { + val matrix = Matrix.square( + 3.0, 1.0, + 1.0, 3.0 + ) + + val decomposed = LUSolver.real.decompose(matrix) + val decomposition = decomposed.getFeature>()!! + + //Check determinant + assertEquals(8.0, decomposition.determinant) + + //Check decomposition + with(MatrixContext.real) { + assertEquals(decomposition.p dot matrix, decomposition.l dot decomposition.u) + } + + val inverted = LUSolver.real.inverse(decomposed) + + val expected = Matrix.square( + 0.375, -0.125, + -0.125, 0.375 + ) + + assertEquals(expected, inverted) + } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt index 0d33204df..aedbf6655 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt @@ -6,11 +6,9 @@ import kotlin.test.assertEquals class RealFieldTest { @Test fun testSqrt() { - - //fails because KT-27586 val sqrt = with(RealField) { - sqrt( 25 * one) + sqrt(25 * one) } - assertEquals(5.0, sqrt.value) + assertEquals(5.0, sqrt) } } \ No newline at end of file diff --git a/kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt similarity index 57% rename from kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt index 350f06848..454683dac 100644 --- a/kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt @@ -1,17 +1,14 @@ package scientifik.kmath.structures -import org.junit.Test import scientifik.kmath.operations.Complex +import scientifik.kmath.operations.complex +import kotlin.test.Test import kotlin.test.assertEquals class ComplexBufferSpecTest { @Test fun testComplexBuffer() { - val buffer = Complex.createBuffer(20) - (0 until 20).forEach { - buffer[it] = Complex(it.toDouble(), -it.toDouble()) - } - + val buffer = Buffer.complex(20) { Complex(it.toDouble(), -it.toDouble()) } assertEquals(Complex(5.0, -5.0), buffer[5]) } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/GenericNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/GenericNDFieldTest.kt deleted file mode 100644 index 338f5f052..000000000 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/GenericNDFieldTest.kt +++ /dev/null @@ -1,16 +0,0 @@ -package scientifik.kmath.structures - -import scientifik.kmath.operations.DoubleField -import scientifik.kmath.structures.NDArrays.create -import kotlin.test.Test -import kotlin.test.assertEquals - - -class GenericNDFieldTest{ - @Test - fun testStrides(){ - val ndArray = create(DoubleField, intArrayOf(10,10)){(it[0]+it[1]).toDouble()} - assertEquals(ndArray[5,5], 10.0) - } - -} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt new file mode 100644 index 000000000..39cce5c67 --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt @@ -0,0 +1,13 @@ +package scientifik.kmath.structures + +import kotlin.test.Test +import kotlin.test.assertEquals + + +class NDFieldTest { + @Test + fun testStrides() { + val ndArray = NDElement.real(intArrayOf(10, 10)) { (it[0] + it[1]).toDouble() } + assertEquals(ndArray[5, 5], 10.0) + } +} \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt index 990adb0c7..60f1f9979 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt @@ -1,16 +1,15 @@ package scientifik.kmath.structures import scientifik.kmath.operations.Norm -import scientifik.kmath.structures.NDArrays.produceReal -import scientifik.kmath.structures.NDArrays.real2DArray +import scientifik.kmath.structures.NDElement.Companion.real2D import kotlin.math.abs import kotlin.math.pow import kotlin.test.Test import kotlin.test.assertEquals class NumberNDFieldTest { - val array1 = real2DArray(3, 3) { i, j -> (i + j).toDouble() } - val array2 = real2DArray(3, 3) { i, j -> (i - j).toDouble() } + val array1 = real2D(3, 3) { i, j -> (i + j).toDouble() } + val array2 = real2D(3, 3) { i, j -> (i - j).toDouble() } @Test fun testSum() { @@ -27,7 +26,7 @@ class NumberNDFieldTest { @Test fun testGeneration() { - val array = real2DArray(3, 3) { i, j -> (i * 10 + j).toDouble() } + val array = real2D(3, 3) { i, j -> (i * 10 + j).toDouble() } for (i in 0..2) { for (j in 0..2) { @@ -51,15 +50,20 @@ class NumberNDFieldTest { assertEquals(2.0, result[0, 2]) } - object L2Norm : Norm, Double> { - override fun norm(arg: NDElement): Double { - return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() }) + @Test + fun combineTest() { + val division = array1.combine(array2, Double::div) + } + + object L2Norm : Norm, Double> { + override fun norm(arg: NDStructure): Double { + return kotlin.math.sqrt(arg.elements().sumByDouble { it.second.toDouble() }) } } @Test fun testInternalContext() { - produceReal(array1.shape) { + NDField.real(*array1.shape).run { with(L2Norm) { 1 + norm(array1) + exp(array2) } diff --git a/kmath-core/src/jsMain/kotlin/scientifik/kmath/histogram/Counters.kt b/kmath-core/src/jsMain/kotlin/scientifik/kmath/histogram/Counters.kt deleted file mode 100644 index 3c2915587..000000000 --- a/kmath-core/src/jsMain/kotlin/scientifik/kmath/histogram/Counters.kt +++ /dev/null @@ -1,16 +0,0 @@ -package scientifik.kmath.histogram - -actual class LongCounter{ - private var sum: Long = 0 - actual fun decrement() {sum--} - actual fun increment() {sum++} - actual fun reset() {sum = 0} - actual fun sum(): Long = sum - actual fun add(l: Long) {sum+=l} -} -actual class DoubleCounter{ - private var sum: Double = 0.0 - actual fun reset() {sum = 0.0} - actual fun sum(): Double = sum - actual fun add(d: Double) {sum+=d} -} \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt deleted file mode 100644 index a5cc8c5de..000000000 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt +++ /dev/null @@ -1,55 +0,0 @@ -package scientifik.kmath.structures - -import java.nio.ByteBuffer - - -/** - * A specification for serialization and deserialization objects to buffer - */ -interface BufferSpec { - fun fromBuffer(buffer: ByteBuffer): T - fun toBuffer(value: T): ByteBuffer -} - -/** - * A [BufferSpec] with fixed unit size. Allows storage of any object without boxing. - */ -interface FixedSizeBufferSpec : BufferSpec { - val unitSize: Int - - /** - * Read an object from buffer in current position - */ - fun ByteBuffer.readObject(): T { - val buffer = ByteArray(unitSize) - get(buffer) - return fromBuffer(ByteBuffer.wrap(buffer)) - } - - /** - * Read an object from buffer in given index (not buffer position - */ - fun ByteBuffer.readObject(index: Int): T { - val dup = duplicate() - dup.position(index*unitSize) - return dup.readObject() - } - - /** - * Write object to [ByteBuffer] in current buffer position - */ - fun ByteBuffer.writeObject(obj: T) { - val buffer = toBuffer(obj).apply { rewind() } - assert(buffer.limit() == unitSize) - put(buffer) - } - - /** - * Put an object in given index - */ - fun ByteBuffer.writeObject(index: Int, obj: T) { - val dup = duplicate() - dup.position(index*unitSize) - dup.writeObject(obj) - } -} \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt deleted file mode 100644 index 68e989e81..000000000 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt +++ /dev/null @@ -1,25 +0,0 @@ -package scientifik.kmath.structures - -import scientifik.kmath.operations.Complex -import java.nio.ByteBuffer - -object ComplexBufferSpec : FixedSizeBufferSpec { - override val unitSize: Int = 16 - - override fun fromBuffer(buffer: ByteBuffer): Complex { - val re = buffer.getDouble(0) - val im = buffer.getDouble(8) - return Complex(re, im) - } - - override fun toBuffer(value: Complex): ByteBuffer = ByteBuffer.allocate(16).apply { - putDouble(value.re) - putDouble(value.im) - } -} - -/** - * Create a mutable buffer which ignores boxing - */ -fun Complex.Companion.createBuffer(size: Int) = ObjectBuffer.create(ComplexBufferSpec, size) - diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt deleted file mode 100644 index b50ed9674..000000000 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt +++ /dev/null @@ -1,28 +0,0 @@ -package scientifik.kmath.structures - -import java.nio.ByteBuffer - -class ObjectBuffer(private val buffer: ByteBuffer, private val spec: FixedSizeBufferSpec) : MutableBuffer { - override val size: Int - get() = buffer.limit() / spec.unitSize - - override fun get(index: Int): T = with(spec) { buffer.readObject(index) } - - override fun iterator(): Iterator = (0 until size).asSequence().map { get(it) }.iterator() - - override fun set(index: Int, value: T) = with(spec) { buffer.writeObject(index, value) } - - override fun copy(): MutableBuffer { - val dup = buffer.duplicate() - val copy = ByteBuffer.allocate(dup.capacity()) - dup.rewind() - copy.put(dup) - copy.flip() - return ObjectBuffer(copy, spec) - } - - companion object { - fun create(spec: FixedSizeBufferSpec, size: Int) = - ObjectBuffer(ByteBuffer.allocate(size * spec.unitSize), spec) - } -} \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt deleted file mode 100644 index 761cfc2db..000000000 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt +++ /dev/null @@ -1,23 +0,0 @@ -package scientifik.kmath.structures - -import scientifik.kmath.operations.Real -import java.nio.ByteBuffer - -object RealBufferSpec : FixedSizeBufferSpec { - override val unitSize: Int = 8 - - override fun fromBuffer(buffer: ByteBuffer): Real = Real(buffer.double) - - override fun toBuffer(value: Real): ByteBuffer = ByteBuffer.allocate(8).apply { putDouble(value.value) } -} - -object DoubleBufferSpec : FixedSizeBufferSpec { - override val unitSize: Int = 8 - - override fun fromBuffer(buffer: ByteBuffer): Double = buffer.double - - override fun toBuffer(value: Double): ByteBuffer = ByteBuffer.allocate(8).apply { putDouble(value) } -} - -fun Double.Companion.createBuffer(size: Int) = ObjectBuffer.create(DoubleBufferSpec, size) -fun Real.Companion.createBuffer(size: Int) = ObjectBuffer.create(RealBufferSpec, size) \ No newline at end of file diff --git a/kmath-coroutines/build.gradle b/kmath-coroutines/build.gradle deleted file mode 100644 index 7149e8bb9..000000000 --- a/kmath-coroutines/build.gradle +++ /dev/null @@ -1,42 +0,0 @@ -plugins { - id "org.jetbrains.kotlin.multiplatform" -} - -kotlin { - targets { - fromPreset(presets.jvm, 'jvm') - // For ARM, preset should be changed to presets.iosArm32 or presets.iosArm64 - // For Linux, preset should be changed to e.g. presets.linuxX64 - // For MacOS, preset should be changed to e.g. presets.macosX64 - //fromPreset(presets.mingwX64, 'mingw') - } - sourceSets { - commonMain { - dependencies { - api project(":kmath-core") - api "org.jetbrains.kotlinx:kotlinx-coroutines-core-common:$coroutinesVersion" - } - } - commonTest { - dependencies { - api 'org.jetbrains.kotlin:kotlin-test-common' - api 'org.jetbrains.kotlin:kotlin-test-annotations-common' - } - } - jvmMain { - dependencies { - api "org.jetbrains.kotlinx:kotlinx-coroutines-core:$coroutinesVersion" - } - } - jvmTest { - dependencies { - implementation 'org.jetbrains.kotlin:kotlin-test' - implementation 'org.jetbrains.kotlin:kotlin-test-junit' - } - } -// mingwMain { -// } -// mingwTest { -// } - } -} diff --git a/kmath-coroutines/build.gradle.kts b/kmath-coroutines/build.gradle.kts new file mode 100644 index 000000000..c73d4e4dc --- /dev/null +++ b/kmath-coroutines/build.gradle.kts @@ -0,0 +1,46 @@ +plugins { + kotlin("multiplatform") +} + +val coroutinesVersion: String by rootProject.extra + +kotlin { + jvm() + js() + + sourceSets { + val commonMain by getting { + dependencies { + api(project(":kmath-core")) + api("org.jetbrains.kotlinx:kotlinx-coroutines-core-common:$coroutinesVersion") + } + } + val commonTest by getting { + dependencies { + implementation(kotlin("test-common")) + implementation(kotlin("test-annotations-common")) + } + } + val jvmMain by getting { + dependencies { + api("org.jetbrains.kotlinx:kotlinx-coroutines-core:$coroutinesVersion") + } + } + val jvmTest by getting { + dependencies { + implementation(kotlin("test")) + implementation(kotlin("test-junit")) + } + } + val jsMain by getting { + dependencies { + api("org.jetbrains.kotlinx:kotlinx-coroutines-core-js:$coroutinesVersion") + } + } + val jsTest by getting { + dependencies { + implementation(kotlin("test-js")) + } + } + } +} diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt index a837cde01..f6f21af09 100644 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt +++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt @@ -6,6 +6,4 @@ import kotlinx.coroutines.Dispatchers import kotlin.coroutines.CoroutineContext import kotlin.coroutines.EmptyCoroutineContext -expect fun runBlocking(context: CoroutineContext = EmptyCoroutineContext, function: suspend CoroutineScope.()->R): R - val Dispatchers.Math: CoroutineDispatcher get() = Dispatchers.Default \ No newline at end of file diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt deleted file mode 100644 index 823245685..000000000 --- a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt +++ /dev/null @@ -1,79 +0,0 @@ -package scientifik.kmath.structures - -import kotlinx.coroutines.* -import scientifik.kmath.operations.Field - -class LazyNDField>(shape: IntArray, field: F, val scope: CoroutineScope = GlobalScope) : NDField(shape, field) { - - override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure = LazyNDStructure(this) { initializer(field, it) } - - - override fun add(a: NDElement, b: NDElement): NDElement { - return LazyNDStructure(this) { index -> - val aDeferred = a.deferred(index) - val bDeferred = b.deferred(index) - aDeferred.await() + bDeferred.await() - } - } - - override fun multiply(a: NDElement, k: Double): NDElement { - return LazyNDStructure(this) { index -> a.await(index) * k } - } - - override fun multiply(a: NDElement, b: NDElement): NDElement { - return LazyNDStructure(this) { index -> - val aDeferred = a.deferred(index) - val bDeferred = b.deferred(index) - aDeferred.await() * bDeferred.await() - } - } - - override fun divide(a: NDElement, b: NDElement): NDElement { - return LazyNDStructure(this) { index -> - val aDeferred = a.deferred(index) - val bDeferred = b.deferred(index) - aDeferred.await() / bDeferred.await() - } - } -} - -class LazyNDStructure>(override val context: LazyNDField, val function: suspend F.(IntArray) -> T) : NDElement, NDStructure { - override val self: NDElement get() = this - override val shape: IntArray get() = context.shape - - private val cache = HashMap>() - - fun deferred(index: IntArray) = cache.getOrPut(index) { context.scope.async(context = Dispatchers.Math) { function.invoke(context.field, index) } } - - suspend fun await(index: IntArray): T = deferred(index).await() - - override fun get(index: IntArray): T = runBlocking { - deferred(index).await() - } - - override fun elements(): Sequence> { - val strides = DefaultStrides(shape) - return strides.indices().map { index -> index to runBlocking { await(index) } } - } -} - -fun NDElement.deferred(index: IntArray) = if (this is LazyNDStructure) this.deferred(index) else CompletableDeferred(get(index)) - -suspend fun NDElement.await(index: IntArray) = if (this is LazyNDStructure) this.await(index) else get(index) - -fun > NDElement.lazy(scope: CoroutineScope = GlobalScope): LazyNDStructure { - return if (this is LazyNDStructure) { - this - } else { - val context = LazyNDField(context.shape, context.field) - LazyNDStructure(context) { get(it) } - } -} - -inline fun > LazyNDStructure.transformIndexed(crossinline action: suspend F.(IntArray, T) -> T) = LazyNDStructure(context) { index -> - action.invoke(this, index, await(index)) -} - -inline fun > LazyNDStructure.transform(crossinline action: suspend F.(T) -> T) = LazyNDStructure(context) { index -> - action.invoke(this, await(index)) -} \ No newline at end of file diff --git a/kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt b/kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt deleted file mode 100644 index 403fe2d31..000000000 --- a/kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt +++ /dev/null @@ -1,20 +0,0 @@ -package scientifik.kmath.structures - -import scientifik.kmath.operations.IntField -import kotlin.test.Test -import kotlin.test.assertEquals - - -class LazyNDFieldTest { - @Test - fun testLazyStructure() { - var counter = 0 - val regularStructure = NDArrays.create(IntField, intArrayOf(2, 2, 2)) { it[0] + it[1] - it[2] } - val result = (regularStructure.lazy() + 2).transform { - counter++ - it * it - } - assertEquals(4, result[0,0,0]) - assertEquals(1, counter) - } -} \ No newline at end of file diff --git a/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt new file mode 100644 index 000000000..b4832827d --- /dev/null +++ b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt @@ -0,0 +1,47 @@ +package scientifik.kmath.structures + +import kotlinx.coroutines.* + +class LazyNDStructure( + val scope: CoroutineScope, + override val shape: IntArray, + val function: suspend (IntArray) -> T +) : NDStructure { + + private val cache = HashMap>() + + fun deferred(index: IntArray) = cache.getOrPut(index) { + scope.async(context = Dispatchers.Math) { + function(index) + } + } + + suspend fun await(index: IntArray): T = deferred(index).await() + + override fun get(index: IntArray): T = runBlocking { + deferred(index).await() + } + + override fun elements(): Sequence> { + val strides = DefaultStrides(shape) + val res = runBlocking { + strides.indices().toList().map { index -> index to await(index) } + } + return res.asSequence() + } +} + +fun NDStructure.deferred(index: IntArray) = + if (this is LazyNDStructure) this.deferred(index) else CompletableDeferred(get(index)) + +suspend fun NDStructure.await(index: IntArray) = + if (this is LazyNDStructure) this.await(index) else get(index) + +/** + * PENDING would benifit from KEEP-176 + */ +fun NDStructure.mapAsyncIndexed(scope: CoroutineScope, function: suspend (T, index: IntArray) -> R) = + LazyNDStructure(scope, shape) { index -> function(get(index), index) } + +fun NDStructure.mapAsync(scope: CoroutineScope, function: suspend (T) -> R) = + LazyNDStructure(scope, shape) { index -> function(get(index)) } \ No newline at end of file diff --git a/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt deleted file mode 100644 index 72a764729..000000000 --- a/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt +++ /dev/null @@ -1,6 +0,0 @@ -package scientifik.kmath.structures - -import kotlinx.coroutines.CoroutineScope -import kotlin.coroutines.CoroutineContext - -actual fun runBlocking(context: CoroutineContext, function: suspend CoroutineScope.() -> R): R = kotlinx.coroutines.runBlocking(context, function) \ No newline at end of file diff --git a/kmath-histograms/build.gradle.kts b/kmath-histograms/build.gradle.kts new file mode 100644 index 000000000..81b3fb83f --- /dev/null +++ b/kmath-histograms/build.gradle.kts @@ -0,0 +1,34 @@ +plugins { + kotlin("multiplatform") +} + +kotlin { + jvm() + js() + + sourceSets { + + val commonMain by getting { + dependencies { + api(project(":kmath-core")) + } + } + val commonTest by getting { + dependencies { + implementation(kotlin("test-common")) + implementation(kotlin("test-annotations-common")) + } + } + val jvmTest by getting { + dependencies { + implementation(kotlin("test")) + implementation(kotlin("test-junit")) + } + } + val jsTest by getting { + dependencies { + implementation(kotlin("test-js")) + } + } + } +} \ No newline at end of file diff --git a/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/Counters.kt b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/Counters.kt new file mode 100644 index 000000000..9c7de3303 --- /dev/null +++ b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/Counters.kt @@ -0,0 +1,20 @@ +package scientifik.kmath.histogram + +/* + * Common representation for atomic counters + * TODO replace with atomics + */ + +expect class LongCounter() { + fun decrement() + fun increment() + fun reset() + fun sum(): Long + fun add(l: Long) +} + +expect class DoubleCounter() { + fun reset() + fun sum(): Double + fun add(d: Double) +} \ No newline at end of file diff --git a/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt new file mode 100644 index 000000000..329af72a1 --- /dev/null +++ b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt @@ -0,0 +1,64 @@ +package scientifik.kmath.histogram + +import scientifik.kmath.linear.Point +import scientifik.kmath.structures.ArrayBuffer +import scientifik.kmath.structures.DoubleBuffer + +/** + * A simple geometric domain + * TODO move to geometry module + */ +interface Domain { + operator fun contains(vector: Point): Boolean + val dimension: Int +} + +/** + * The bin in the histogram. The histogram is by definition always done in the real space + */ +interface Bin : Domain { + /** + * The value of this bin + */ + val value: Number + val center: Point +} + +interface Histogram> : Iterable { + + /** + * Find existing bin, corresponding to given coordinates + */ + operator fun get(point: Point): B? + + /** + * Dimension of the histogram + */ + val dimension: Int + +} + +interface MutableHistogram> : Histogram { + + /** + * Increment appropriate bin + */ + fun putWithWeight(point: Point, weight: Double) + + fun put(point: Point) = putWithWeight(point, 1.0) +} + +fun MutableHistogram.put(vararg point: T) = put(ArrayBuffer(point)) + +fun MutableHistogram.put(vararg point: Number) = + put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray())) + +fun MutableHistogram.put(vararg point: Double) = put(DoubleBuffer(point)) + +fun MutableHistogram.fill(sequence: Iterable>) = sequence.forEach { put(it) } + +/** + * Pass a sequence builder into histogram + */ +fun MutableHistogram.fill(buider: suspend SequenceScope>.() -> Unit) = + fill(sequence(buider).asIterable()) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt similarity index 50% rename from kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt rename to kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt index e48979430..75296ba27 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt +++ b/kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt @@ -1,43 +1,66 @@ package scientifik.kmath.histogram +import scientifik.kmath.linear.Point import scientifik.kmath.linear.toVector +import scientifik.kmath.operations.SpaceOperations import scientifik.kmath.structures.* import kotlin.math.floor -private operator fun RealPoint.minus(other: RealPoint) = ListBuffer((0 until size).map { get(it) - other[it] }) -private inline fun Buffer.mapIndexed(crossinline mapper: (Int, Double) -> T): Sequence = (0 until size).asSequence().map { mapper(it, get(it)) } +data class BinDef>(val space: SpaceOperations>, val center: Point, val sizes: Point) { + fun contains(vector: Point): Boolean { + if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}") + val upper = space.run { center + sizes / 2.0 } + val lower = space.run { center - sizes / 2.0 } + return vector.asSequence().mapIndexed { i, value -> + value in lower[i]..upper[i] + }.all { it } + } +} + + +class MultivariateBin>(val def: BinDef, override val value: Number) : Bin { + + override fun contains(vector: Point): Boolean = def.contains(vector) + + override val dimension: Int + get() = def.center.size + + override val center: Point + get() = def.center + +} /** * Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions. */ -class FastHistogram( - private val lower: RealPoint, - private val upper: RealPoint, - private val binNums: IntArray = IntArray(lower.size) { 20 } -) : MutableHistogram> { +class RealHistogram( + private val lower: Buffer, + private val upper: Buffer, + private val binNums: IntArray = IntArray(lower.size) { 20 } +) : MutableHistogram> { private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 }) - private val values: NDStructure = ndStructure(strides) { LongCounter() } + private val values: NDStructure = NDStructure.auto(strides) { LongCounter() } //private val weight: NDStructure = ndStructure(strides){null} - //TODO optimize binSize performance if needed - private val binSize: RealPoint = ListBuffer((upper - lower).mapIndexed { index, value -> value / binNums[index] }.toList()) + + override val dimension: Int get() = lower.size + + + private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] } init { // argument checks if (lower.size != upper.size) error("Dimension mismatch in histogram lower and upper limits.") if (lower.size != binNums.size) error("Dimension mismatch in bin count.") - if ((upper - lower).asSequence().any { it <= 0 }) error("Range for one of axis is not strictly positive") + if ((0 until dimension).any { upper[it] - lower[it] < 0 }) error("Range for one of axis is not strictly positive") } - override val dimension: Int get() = lower.size - - /** * Get internal [NDStructure] bin index for given axis */ @@ -59,49 +82,41 @@ class FastHistogram( return getValue(getIndex(point)) } - private fun getTemplate(index: IntArray): BinTemplate { + private fun getDef(index: IntArray): BinDef { val center = index.mapIndexed { axis, i -> when (i) { 0 -> Double.NEGATIVE_INFINITY strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY else -> lower[axis] + (i.toDouble() - 0.5) * binSize[axis] } - }.toVector() - return BinTemplate(center, binSize) + }.asBuffer() + return BinDef(RealBufferFieldOperations, center, binSize) } - fun getTemplate(point: Buffer): BinTemplate { - return getTemplate(getIndex(point)) + fun getDef(point: Buffer): BinDef { + return getDef(getIndex(point)) } - override fun get(point: Buffer): PhantomBin? { + override fun get(point: Buffer): MultivariateBin? { val index = getIndex(point) - return PhantomBin(getTemplate(index), getValue(index)) + return MultivariateBin(getDef(index), getValue(index)) } - override fun put(point: Buffer, weight: Double) { + override fun putWithWeight(point: Buffer, weight: Double) { if (weight != 1.0) TODO("Implement weighting") val index = getIndex(point) values[index].increment() } - override fun iterator(): Iterator> = values.elements().map { (index, value) -> - PhantomBin(getTemplate(index), value.sum()) + override fun iterator(): Iterator> = values.elements().map { (index, value) -> + MultivariateBin(getDef(index), value.sum()) }.iterator() /** * Convert this histogram into NDStructure containing bin values but not bin descriptions */ fun asNDStructure(): NDStructure { - return ndStructure(this.values.shape) { values[it].sum() } - } - - /** - * Create a phantom lightweight immutable copy of this histogram - */ - fun asPhantomHistogram(): PhantomHistogram { - val binTemplates = values.elements().associate { (index, _) -> getTemplate(index) to index } - return PhantomHistogram(binTemplates, asNDStructure()) + return NDStructure.auto(this.values.shape) { values[it].sum() } } companion object { @@ -115,8 +130,11 @@ class FastHistogram( *) *``` */ - fun fromRanges(vararg ranges: ClosedFloatingPointRange): FastHistogram { - return FastHistogram(ranges.map { it.start }.toVector(), ranges.map { it.endInclusive }.toVector()) + fun fromRanges(vararg ranges: ClosedFloatingPointRange): RealHistogram { + return RealHistogram( + ranges.map { it.start }.toVector(), + ranges.map { it.endInclusive }.toVector() + ) } /** @@ -128,11 +146,11 @@ class FastHistogram( *) *``` */ - fun fromRanges(vararg ranges: Pair, Int>): FastHistogram { - return FastHistogram( - ListBuffer(ranges.map { it.first.start }), - ListBuffer(ranges.map { it.first.endInclusive }), - ranges.map { it.second }.toIntArray() + fun fromRanges(vararg ranges: Pair, Int>): RealHistogram { + return RealHistogram( + ListBuffer(ranges.map { it.first.start }), + ListBuffer(ranges.map { it.first.endInclusive }), + ranges.map { it.second }.toIntArray() ) } } diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt similarity index 57% rename from kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt rename to kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt index 842a12ae4..678b8eba7 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt +++ b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt @@ -1,5 +1,8 @@ -package scientifik.kmath.histogram +package scietifik.kmath.histogram +import scientifik.kmath.histogram.RealHistogram +import scientifik.kmath.histogram.fill +import scientifik.kmath.histogram.put import scientifik.kmath.linear.Vector import kotlin.random.Random import kotlin.test.Test @@ -10,9 +13,9 @@ import kotlin.test.assertTrue class MultivariateHistogramTest { @Test fun testSinglePutHistogram() { - val histogram = FastHistogram.fromRanges( - (-1.0..1.0), - (-1.0..1.0) + val histogram = RealHistogram.fromRanges( + (-1.0..1.0), + (-1.0..1.0) ) histogram.put(0.55, 0.55) val bin = histogram.find { it.value.toInt() > 0 }!! @@ -22,21 +25,21 @@ class MultivariateHistogramTest { } @Test - fun testSequentialPut(){ - val histogram = FastHistogram.fromRanges( - (-1.0..1.0), - (-1.0..1.0), - (-1.0..1.0) + fun testSequentialPut() { + val histogram = RealHistogram.fromRanges( + (-1.0..1.0), + (-1.0..1.0), + (-1.0..1.0) ) val random = Random(1234) - fun nextDouble() = random.nextDouble(-1.0,1.0) + fun nextDouble() = random.nextDouble(-1.0, 1.0) val n = 10000 histogram.fill { - repeat(n){ - yield(Vector.ofReal(nextDouble(),nextDouble(),nextDouble())) + repeat(n) { + yield(Vector.ofReal(nextDouble(), nextDouble(), nextDouble())) } } assertEquals(n, histogram.sumBy { it.value.toInt() }) diff --git a/kmath-histograms/src/jsMain/kotlin/scientifik/kmath/histogram/Counters.kt b/kmath-histograms/src/jsMain/kotlin/scientifik/kmath/histogram/Counters.kt new file mode 100644 index 000000000..3765220b9 --- /dev/null +++ b/kmath-histograms/src/jsMain/kotlin/scientifik/kmath/histogram/Counters.kt @@ -0,0 +1,33 @@ +package scientifik.kmath.histogram + +actual class LongCounter { + private var sum: Long = 0 + actual fun decrement() { + sum-- + } + + actual fun increment() { + sum++ + } + + actual fun reset() { + sum = 0 + } + + actual fun sum(): Long = sum + actual fun add(l: Long) { + sum += l + } +} + +actual class DoubleCounter { + private var sum: Double = 0.0 + actual fun reset() { + sum = 0.0 + } + + actual fun sum(): Double = sum + actual fun add(d: Double) { + sum += d + } +} \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/Counters.kt b/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/Counters.kt similarity index 100% rename from kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/Counters.kt rename to kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/Counters.kt diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt similarity index 85% rename from kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt rename to kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt index 26ad520e6..53c2da641 100644 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt +++ b/kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt @@ -18,7 +18,7 @@ class UnivariateBin(val position: Double, val size: Double, val counter: LongCou override fun contains(vector: Buffer): Boolean = contains(vector[0]) - internal operator fun inc() = this.also { counter.increment()} + internal operator fun inc() = this.also { counter.increment() } override val dimension: Int get() = 1 } @@ -26,7 +26,8 @@ class UnivariateBin(val position: Double, val size: Double, val counter: LongCou /** * Univariate histogram with log(n) bin search speed */ -class UnivariateHistogram private constructor(private val factory: (Double) -> UnivariateBin) : MutableHistogram { +class UnivariateHistogram private constructor(private val factory: (Double) -> UnivariateBin) : + MutableHistogram { private val bins: TreeMap = TreeMap() @@ -58,7 +59,7 @@ class UnivariateHistogram private constructor(private val factory: (Double) -> U (get(value) ?: createBin(value)).inc() } - override fun put(point: Buffer, weight: Double) { + override fun putWithWeight(point: Buffer, weight: Double) { if (weight != 1.0) TODO("Implement weighting") put(point[0]) } @@ -75,8 +76,14 @@ class UnivariateHistogram private constructor(private val factory: (Double) -> U val sorted = borders.sortedArray() return UnivariateHistogram { value -> when { - value < sorted.first() -> UnivariateBin(Double.NEGATIVE_INFINITY, Double.MAX_VALUE) - value > sorted.last() -> UnivariateBin(Double.POSITIVE_INFINITY, Double.MAX_VALUE) + value < sorted.first() -> UnivariateBin( + Double.NEGATIVE_INFINITY, + Double.MAX_VALUE + ) + value > sorted.last() -> UnivariateBin( + Double.POSITIVE_INFINITY, + Double.MAX_VALUE + ) else -> { val index = (0 until sorted.size).first { value > sorted[it] } val left = sorted[index] diff --git a/kmath-jmh/build.gradle b/kmath-jmh/build.gradle deleted file mode 100644 index 1876d86d7..000000000 --- a/kmath-jmh/build.gradle +++ /dev/null @@ -1,10 +0,0 @@ -plugins { - id "java" - id "kotlin" - id "me.champeau.gradle.jmh" version "0.4.7" -} - -dependencies { - compile project(':kmath-core') - //jmh project(':kmath-core') -} \ No newline at end of file diff --git a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt b/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt deleted file mode 100644 index fe10fbd75..000000000 --- a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt +++ /dev/null @@ -1,49 +0,0 @@ -package scientifik.kmath.structures - -import org.openjdk.jmh.annotations.* -import java.nio.IntBuffer - - -@Warmup(iterations = 1) -@Measurement(iterations = 5) -@State(Scope.Benchmark) -open class ArrayBenchmark { - - lateinit var array: IntArray - lateinit var arrayBuffer: IntBuffer - lateinit var nativeBuffer: IntBuffer - - @Setup - fun setup() { - array = IntArray(10000) { it } - arrayBuffer = IntBuffer.wrap(array) - nativeBuffer = IntBuffer.allocate(10000) - for (i in 0 until 10000) { - nativeBuffer.put(i,i) - } - } - - @Benchmark - fun benchmarkArrayRead() { - var res = 0 - for (i in 1..10000) { - res += array[10000 - i] - } - } - - @Benchmark - fun benchmarkBufferRead() { - var res = 0 - for (i in 1..10000) { - res += arrayBuffer.get(10000 - i) - } - } - - @Benchmark - fun nativeBufferRead() { - var res = 0 - for (i in 1..10000) { - res += nativeBuffer.get(10000 - i) - } - } -} \ No newline at end of file diff --git a/kmath-koma/build.gradle.kts b/kmath-koma/build.gradle.kts new file mode 100644 index 000000000..b95aaf3c8 --- /dev/null +++ b/kmath-koma/build.gradle.kts @@ -0,0 +1,57 @@ +plugins { + kotlin("multiplatform") +} + +repositories { + maven("http://dl.bintray.com/kyonifer/maven") +} + +kotlin { + jvm { + compilations.all { + kotlinOptions { + jvmTarget = "1.8" + freeCompilerArgs += "-progressive" + } + } + } + js() + + sourceSets { + + val commonMain by getting { + dependencies { + api(project(":kmath-core")) + api("com.kyonifer:koma-core-api-common:0.12") + } + } + val commonTest by getting { + dependencies { + implementation(kotlin("test-common")) + implementation(kotlin("test-annotations-common")) + } + } + val jvmMain by getting { + dependencies { + api("com.kyonifer:koma-core-api-jvm:0.12") + } + } + val jvmTest by getting { + dependencies { + implementation(kotlin("test")) + implementation(kotlin("test-junit")) + implementation("com.kyonifer:koma-core-ejml:0.12") + } + } + val jsMain by getting { + dependencies { + api("com.kyonifer:koma-core-api-js:0.12") + } + } + val jsTest by getting { + dependencies { + implementation(kotlin("test-js")) + } + } + } +} \ No newline at end of file diff --git a/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt new file mode 100644 index 000000000..343d5b8b9 --- /dev/null +++ b/kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt @@ -0,0 +1,85 @@ +package scientifik.kmath.linear + +import koma.extensions.fill +import koma.matrix.MatrixFactory + +class KomaMatrixContext(val factory: MatrixFactory>) : MatrixContext, + LinearSolver { + + override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T) = + KomaMatrix(factory.zeros(rows, columns).fill(initializer)) + + fun Matrix.toKoma(): KomaMatrix = if (this is KomaMatrix) { + this + } else { + produce(rowNum, colNum) { i, j -> get(i, j) } + } + + fun Point.toKoma(): KomaVector = if (this is KomaVector) { + this + } else { + KomaVector(factory.zeros(size, 1).fill { i, _ -> get(i) }) + } + + + override fun Matrix.dot(other: Matrix) = + KomaMatrix(this.toKoma().origin * other.toKoma().origin) + + override fun Matrix.dot(vector: Point) = + KomaVector(this.toKoma().origin * vector.toKoma().origin) + + override fun Matrix.unaryMinus() = + KomaMatrix(this.toKoma().origin.unaryMinus()) + + override fun Matrix.plus(b: Matrix) = + KomaMatrix(this.toKoma().origin + b.toKoma().origin) + + override fun Matrix.minus(b: Matrix) = + KomaMatrix(this.toKoma().origin - b.toKoma().origin) + + override fun Matrix.times(value: T) = + KomaMatrix(this.toKoma().origin * value) + + + override fun solve(a: Matrix, b: Matrix) = + KomaMatrix(a.toKoma().origin.solve(b.toKoma().origin)) + + override fun inverse(a: Matrix) = + KomaMatrix(a.toKoma().origin.inv()) +} + +class KomaMatrix(val origin: koma.matrix.Matrix, features: Set? = null) : + Matrix { + override val rowNum: Int get() = origin.numRows() + override val colNum: Int get() = origin.numCols() + + override val features: Set = features ?: setOf( + object : DeterminantFeature { + override val determinant: T get() = origin.det() + }, + object : LUPDecompositionFeature { + private val lup by lazy { origin.LU() } + override val l: Matrix get() = KomaMatrix(lup.second) + override val u: Matrix get() = KomaMatrix(lup.third) + override val p: Matrix get() = KomaMatrix(lup.first) + } + ) + + override fun suggestFeature(vararg features: MatrixFeature): Matrix = + KomaMatrix(this.origin, this.features + features) + + override fun get(i: Int, j: Int): T = origin.getGeneric(i, j) +} + +class KomaVector internal constructor(val origin: koma.matrix.Matrix) : Point { + init { + if (origin.numCols() != 1) error("Only single column matrices are allowed") + } + + override val size: Int get() = origin.numRows() + + override fun get(index: Int): T = origin.getGeneric(index) + + override fun iterator(): Iterator = origin.toIterable().iterator() +} + diff --git a/kmath-memory/build.gradle.kts b/kmath-memory/build.gradle.kts new file mode 100644 index 000000000..f6d10875c --- /dev/null +++ b/kmath-memory/build.gradle.kts @@ -0,0 +1,50 @@ +plugins { + kotlin("multiplatform") +} + +val ioVersion: String by rootProject.extra + + +kotlin { + jvm() + js() + + sourceSets { + val commonMain by getting { + dependencies { + api(kotlin("stdlib")) + } + } + val commonTest by getting { + dependencies { + implementation(kotlin("test-common")) + implementation(kotlin("test-annotations-common")) + } + } + val jvmMain by getting { + dependencies { + api(kotlin("stdlib-jdk8")) + } + } + val jvmTest by getting { + dependencies { + implementation(kotlin("test")) + implementation(kotlin("test-junit")) + } + } + val jsMain by getting { + dependencies { + api(kotlin("stdlib-js")) + } + } + val jsTest by getting { + dependencies { + implementation(kotlin("test-js")) + } + } +// mingwMain { +// } +// mingwTest { +// } + } +} \ No newline at end of file diff --git a/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt b/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt new file mode 100644 index 000000000..f9f938dcc --- /dev/null +++ b/kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt @@ -0,0 +1,71 @@ +package scientifik.memory + +interface Memory { + val size: Int + + /** + * Get a projection of this memory (it reflects the changes in the parent memory block) + */ + fun view(offset: Int, length: Int): Memory + + /** + * Create a copy of this memory, which does not know anything about this memory + */ + fun copy(): Memory + + /** + * Create and possibly register a new reader + */ + fun reader(): MemoryReader + + fun writer(): MemoryWriter + + companion object { + + } +} + +interface MemoryReader { + val memory: Memory + + fun readDouble(offset: Int): Double + fun readFloat(offset: Int): Float + fun readByte(offset: Int): Byte + fun readShort(offset: Int): Short + fun readInt(offset: Int): Int + fun readLong(offset: Int): Long + + fun release() +} + +/** + * Use the memory for read then release the reader + */ +inline fun Memory.read(block: MemoryReader.() -> Unit) { + reader().apply(block).apply { release() } +} + +interface MemoryWriter { + val memory: Memory + + fun writeDouble(offset: Int, value: Double) + fun writeFloat(offset: Int, value: Float) + fun writeByte(offset: Int, value: Byte) + fun writeShort(offset: Int, value: Short) + fun writeInt(offset: Int, value: Int) + fun writeLong(offset: Int, value: Long) + + fun release() +} + +/** + * Use the memory for write then release the writer + */ +inline fun Memory.write(block: MemoryWriter.() -> Unit) { + writer().apply(block).apply { release() } +} + +/** + * Allocate the most effective platform-specific memory + */ +expect fun Memory.Companion.allocate(length: Int): Memory diff --git a/kmath-memory/src/commonMain/kotlin/scientifik/memory/MemorySpec.kt b/kmath-memory/src/commonMain/kotlin/scientifik/memory/MemorySpec.kt new file mode 100644 index 000000000..a3166ecd7 --- /dev/null +++ b/kmath-memory/src/commonMain/kotlin/scientifik/memory/MemorySpec.kt @@ -0,0 +1,34 @@ +package scientifik.memory + +/** + * A specification to read or write custom objects with fixed size in bytes + */ +interface MemorySpec { + /** + * Size of [T] in bytes after serialization + */ + val objectSize: Int + + fun MemoryReader.read(offset: Int): T + fun MemoryWriter.write(offset: Int, value: T) +} + +fun MemoryReader.read(spec: MemorySpec, offset: Int): T = spec.run { read(offset) } +fun MemoryWriter.write(spec: MemorySpec, offset: Int, value: T) = spec.run { write(offset, value) } + +inline fun MemoryReader.readArray(spec: MemorySpec, offset: Int, size: Int) = + Array(size) { i -> + spec.run { + read(offset + i * objectSize) + } + } + +fun MemoryWriter.writeArray(spec: MemorySpec, offset: Int, array: Array) { + spec.run { + for (i in 0 until array.size) { + write(offset + i * objectSize, array[i]) + } + } +} + +//TODO It is possible to add elastic MemorySpec with unknown object size \ No newline at end of file diff --git a/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt b/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt new file mode 100644 index 000000000..843464ab9 --- /dev/null +++ b/kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt @@ -0,0 +1,91 @@ +package scientifik.memory + +import org.khronos.webgl.ArrayBuffer +import org.khronos.webgl.DataView + +/** + * Allocate the most effective platform-specific memory + */ +actual fun Memory.Companion.allocate(length: Int): Memory { + val buffer = ArrayBuffer(length) + return DataViewMemory(DataView(buffer, 0, length)) +} + +class DataViewMemory(val view: DataView) : Memory { + + override val size: Int get() = view.byteLength + + override fun view(offset: Int, length: Int): Memory { + require(offset >= 0) { "offset shouldn't be negative: $offset" } + require(length >= 0) { "length shouldn't be negative: $length" } + if (offset + length > size) { + throw IndexOutOfBoundsException("offset + length > size: $offset + $length > $size") + } + return DataViewMemory(DataView(view.buffer, view.byteOffset + offset, length)) + } + + + override fun copy(): Memory { + TODO("not implemented") //To change body of created functions use File | Settings | File Templates. + } + + private val reader = object : MemoryReader { + override val memory: Memory get() = this@DataViewMemory + + override fun readDouble(offset: Int): Double = view.getFloat64(offset, false) + + override fun readFloat(offset: Int): Float = view.getFloat32(offset, false) + + override fun readByte(offset: Int): Byte = view.getInt8(offset) + + override fun readShort(offset: Int): Short = view.getInt16(offset, false) + + override fun readInt(offset: Int): Int = view.getInt32(offset, false) + + override fun readLong(offset: Int): Long = (view.getInt32(offset, false).toLong() shl 32) or + view.getInt32(offset + 4, false).toLong() + + override fun release() { + // does nothing on JS because of GC + } + } + + override fun reader(): MemoryReader = reader + + private val writer = object : MemoryWriter { + override val memory: Memory get() = this@DataViewMemory + + override fun writeDouble(offset: Int, value: Double) { + view.setFloat64(offset, value, false) + } + + override fun writeFloat(offset: Int, value: Float) { + view.setFloat32(offset, value, false) + } + + override fun writeByte(offset: Int, value: Byte) { + view.setInt8(offset, value) + } + + override fun writeShort(offset: Int, value: Short) { + view.setUint16(offset, value, false) + } + + override fun writeInt(offset: Int, value: Int) { + view.setInt32(offset, value, false) + } + + override fun writeLong(offset: Int, value: Long) { + view.setInt32(offset, (value shr 32).toInt(), littleEndian = false) + view.setInt32(offset + 4, (value and 0xffffffffL).toInt(), littleEndian = false) + } + + override fun release() { + //does nothing on JS + } + + } + + override fun writer(): MemoryWriter = writer + +} \ No newline at end of file diff --git a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt new file mode 100644 index 000000000..30a87228a --- /dev/null +++ b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt @@ -0,0 +1,93 @@ +package scientifik.memory + +import java.nio.ByteBuffer + + +/** + * Allocate the most effective platform-specific memory + */ +actual fun Memory.Companion.allocate(length: Int): Memory { + val buffer = ByteBuffer.allocate(length) + return ByteBufferMemory(buffer) +} + +class ByteBufferMemory( + val buffer: ByteBuffer, + val startOffset: Int = 0, + override val size: Int = buffer.limit() +) : Memory { + + + @Suppress("NOTHING_TO_INLINE") + private inline fun position(o: Int): Int = startOffset + o + + override fun view(offset: Int, length: Int): Memory { + if (offset + length > size) error("Selecting a Memory view outside of memory range") + return ByteBufferMemory(buffer, position(offset), length) + } + + override fun copy(): Memory { + val copy = ByteBuffer.allocate(buffer.capacity()) + buffer.rewind() + copy.put(buffer) + copy.flip() + return ByteBufferMemory(copy) + + } + + private val reader = object : MemoryReader { + override val memory: Memory get() = this@ByteBufferMemory + + override fun readDouble(offset: Int) = buffer.getDouble(position(offset)) + + override fun readFloat(offset: Int) = buffer.getFloat(position(offset)) + + override fun readByte(offset: Int) = buffer.get(position(offset)) + + override fun readShort(offset: Int) = buffer.getShort(position(offset)) + + override fun readInt(offset: Int) = buffer.getInt(position(offset)) + + override fun readLong(offset: Int) = buffer.getLong(position(offset)) + + override fun release() { + //does nothing on JVM + } + } + + override fun reader(): MemoryReader = reader + + private val writer = object : MemoryWriter { + override val memory: Memory get() = this@ByteBufferMemory + + override fun writeDouble(offset: Int, value: Double) { + buffer.putDouble(position(offset), value) + } + + override fun writeFloat(offset: Int, value: Float) { + buffer.putFloat(position(offset), value) + } + + override fun writeByte(offset: Int, value: Byte) { + buffer.put(position(offset), value) + } + + override fun writeShort(offset: Int, value: Short) { + buffer.putShort(position(offset), value) + } + + override fun writeInt(offset: Int, value: Int) { + buffer.putInt(position(offset), value) + } + + override fun writeLong(offset: Int, value: Long) { + buffer.putLong(position(offset), value) + } + + override fun release() { + //does nothing on JVM + } + } + + override fun writer(): MemoryWriter = writer +} \ No newline at end of file diff --git a/kmath-sequential/build.gradle.kts b/kmath-sequential/build.gradle.kts new file mode 100644 index 000000000..391504f23 --- /dev/null +++ b/kmath-sequential/build.gradle.kts @@ -0,0 +1,53 @@ +plugins { + kotlin("multiplatform") + id("kotlinx-atomicfu") version "0.12.1" +} + +val atomicfuVersion: String by rootProject.extra + +kotlin { + jvm () + //js() + + sourceSets { + val commonMain by getting { + dependencies { + api(project(":kmath-core")) + api(project(":kmath-coroutines")) + compileOnly("org.jetbrains.kotlinx:atomicfu-common:$atomicfuVersion") + } + } + val commonTest by getting { + dependencies { + implementation(kotlin("test-common")) + implementation(kotlin("test-annotations-common")) + } + } + val jvmMain by getting { + dependencies { + compileOnly("org.jetbrains.kotlinx:atomicfu:$atomicfuVersion") + } + } + val jvmTest by getting { + dependencies { + implementation(kotlin("test")) + implementation(kotlin("test-junit")) + } + } +// val jsMain by getting { +// dependencies { +// compileOnly("org.jetbrains.kotlinx:atomicfu-js:$atomicfuVersion") +// } +// } +// val jsTest by getting { +// dependencies { +// implementation(kotlin("test-js")) +// } +// } + + } +} + +atomicfu { + variant = "VH" +} \ No newline at end of file diff --git a/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/BufferStreaming.kt b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/BufferStreaming.kt new file mode 100644 index 000000000..b98f48186 --- /dev/null +++ b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/BufferStreaming.kt @@ -0,0 +1,82 @@ +package scientifik.kmath.sequential + +import kotlinx.coroutines.CoroutineScope +import kotlinx.coroutines.ExperimentalCoroutinesApi +import kotlinx.coroutines.channels.Channel +import kotlinx.coroutines.channels.produce +import kotlinx.coroutines.isActive +import kotlinx.coroutines.sync.Mutex +import kotlinx.coroutines.sync.withLock +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.BufferFactory + +/** + * A processor that collects incoming elements into fixed size buffers + */ +@ExperimentalCoroutinesApi +class JoinProcessor( + scope: CoroutineScope, + bufferSize: Int, + bufferFactory: BufferFactory = Buffer.Companion::boxing +) : AbstractProcessor>(scope) { + + private val input = Channel(bufferSize) + + private val output = produce(coroutineContext) { + val list = ArrayList(bufferSize) + while (isActive) { + list.clear() + repeat(bufferSize) { + list.add(input.receive()) + } + val buffer = bufferFactory(bufferSize) { list[it] } + send(buffer) + } + } + + override suspend fun receive(): Buffer = output.receive() + + override suspend fun send(value: T) { + input.send(value) + } +} + +/** + * A processor that splits incoming buffers into individual elements + */ +class SplitProcessor(scope: CoroutineScope) : AbstractProcessor, T>(scope) { + + private val input = Channel>() + + private val mutex = Mutex() + + private var currentBuffer: Buffer? = null + + private var pos = 0 + + + override suspend fun receive(): T { + mutex.withLock { + while (currentBuffer == null || pos == currentBuffer!!.size) { + currentBuffer = input.receive() + pos = 0 + } + return currentBuffer!![pos].also { pos++ } + } + } + + override suspend fun send(value: Buffer) { + input.send(value) + } +} + +@ExperimentalCoroutinesApi +fun Producer.chunked(chunkSize: Int, bufferFactory: BufferFactory) = + JoinProcessor(this, chunkSize, bufferFactory).also { connect(it) } + +@ExperimentalCoroutinesApi +inline fun Producer.chunked(chunkSize: Int) = + JoinProcessor(this, chunkSize, Buffer.Companion::auto).also { connect(it) } + + + diff --git a/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Chain.kt b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Chain.kt new file mode 100644 index 000000000..7633b2223 --- /dev/null +++ b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Chain.kt @@ -0,0 +1,151 @@ +/* + * Copyright 2018 Alexander Nozik. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package scientifik.kmath.sequential + +import kotlinx.atomicfu.atomic +import kotlinx.coroutines.CoroutineScope +import kotlinx.coroutines.ExperimentalCoroutinesApi +import kotlinx.coroutines.channels.ReceiveChannel +import kotlinx.coroutines.channels.produce +import kotlinx.coroutines.isActive + + +/** + * A not-necessary-Markov chain of some type + * @param R - the chain element type + */ +interface Chain { + /** + * Last value of the chain. Returns null if [next] was not called + */ + val value: R? + + /** + * Generate next value, changing state if needed + */ + suspend fun next(): R + + /** + * Create a copy of current chain state. Consuming resulting chain does not affect initial chain + */ + fun fork(): Chain + +} + +/** + * Chain as a coroutine receive channel + */ +@ExperimentalCoroutinesApi +fun Chain.asChannel(scope: CoroutineScope): ReceiveChannel = scope.produce { while (isActive) send(next()) } + +fun Iterator.asChain(): Chain = SimpleChain { next() } +fun Sequence.asChain(): Chain = iterator().asChain() + + +/** + * Map the chain result using suspended transformation. Initial chain result can no longer be safely consumed + * since mapped chain consumes tokens. Accepts regular transformation function + */ +fun Chain.map(func: (T) -> R): Chain { + val parent = this; + return object : Chain { + override val value: R? get() = parent.value?.let(func) + + override suspend fun next(): R { + return func(parent.next()) + } + + override fun fork(): Chain { + return parent.fork().map(func) + } + } +} + +/** + * A simple chain of independent tokens + */ +class SimpleChain(private val gen: suspend () -> R) : Chain { + private val atomicValue = atomic(null) + override val value: R? get() = atomicValue.value + + override suspend fun next(): R = gen().also { atomicValue.lazySet(it) } + + override fun fork(): Chain = this +} + +//TODO force forks on mapping operations? + +/** + * A stateless Markov chain + */ +class MarkovChain(private val seed: () -> R, private val gen: suspend (R) -> R) : + Chain { + + constructor(seed: R, gen: suspend (R) -> R) : this({ seed }, gen) + + private val atomicValue = atomic(null) + override val value: R get() = atomicValue.value ?: seed() + + override suspend fun next(): R { + val newValue = gen(value) + atomicValue.lazySet(newValue) + return value + } + + override fun fork(): Chain { + return MarkovChain(value, gen) + } +} + +/** + * A chain with possibly mutable state. The state must not be changed outside the chain. Two chins should never share the state + * @param S - the state of the chain + */ +class StatefulChain( + private val state: S, + private val seed: S.() -> R, + private val gen: suspend S.(R) -> R +) : Chain { + + constructor(state: S, seed: R, gen: suspend S.(R) -> R) : this(state, { seed }, gen) + + private val atomicValue = atomic(null) + override val value: R get() = atomicValue.value ?: seed(state) + + override suspend fun next(): R { + val newValue = gen(state, value) + atomicValue.lazySet(newValue) + return value + } + + override fun fork(): Chain { + throw RuntimeException("Fork not supported for stateful chain") + } +} + +/** + * A chain that repeats the same value + */ +class ConstantChain(override val value: T) : Chain { + override suspend fun next(): T { + return value + } + + override fun fork(): Chain { + return this + } +} \ No newline at end of file diff --git a/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Cumulative.kt b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Cumulative.kt new file mode 100644 index 000000000..b0e1e9ac5 --- /dev/null +++ b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Cumulative.kt @@ -0,0 +1,73 @@ +package scientifik.kmath.sequential + +import scientifik.kmath.operations.Space +import kotlin.jvm.JvmName + + +/** + * Generic cumulative operation on iterator + * @param T type of initial iterable + * @param R type of resulting iterable + * @param initial lazy evaluated + */ +fun Iterator.cumulative(initial: R, operation: (T, R) -> R): Iterator = object : Iterator { + var state: R = initial + override fun hasNext(): Boolean = this@cumulative.hasNext() + + override fun next(): R { + state = operation.invoke(this@cumulative.next(), state) + return state + } +} + +fun Iterable.cumulative(initial: R, operation: (T, R) -> R): Iterable = object : Iterable { + override fun iterator(): Iterator = this@cumulative.iterator().cumulative(initial, operation) +} + +fun Sequence.cumulative(initial: R, operation: (T, R) -> R): Sequence = object : Sequence { + override fun iterator(): Iterator = this@cumulative.iterator().cumulative(initial, operation) +} + +fun List.cumulative(initial: R, operation: (T, R) -> R): List = + this.iterator().cumulative(initial, operation).asSequence().toList() + +//Cumulative sum + +fun Iterable.cumulativeSum(space: Space) = with(space) { + cumulative(zero) { element: T, sum: T -> sum + element } +} + +@JvmName("cumulativeSumOfDouble") +fun Iterable.cumulativeSum() = this.cumulative(0.0) { element, sum -> sum + element } + +@JvmName("cumulativeSumOfInt") +fun Iterable.cumulativeSum() = this.cumulative(0) { element, sum -> sum + element } + +@JvmName("cumulativeSumOfLong") +fun Iterable.cumulativeSum() = this.cumulative(0L) { element, sum -> sum + element } + +fun Sequence.cumulativeSum(space: Space) = with(space) { + cumulative(zero) { element: T, sum: T -> sum + element } +} + +@JvmName("cumulativeSumOfDouble") +fun Sequence.cumulativeSum() = this.cumulative(0.0) { element, sum -> sum + element } + +@JvmName("cumulativeSumOfInt") +fun Sequence.cumulativeSum() = this.cumulative(0) { element, sum -> sum + element } + +@JvmName("cumulativeSumOfLong") +fun Sequence.cumulativeSum() = this.cumulative(0L) { element, sum -> sum + element } + +fun List.cumulativeSum(space: Space) = with(space) { + cumulative(zero) { element: T, sum: T -> sum + element } +} + +@JvmName("cumulativeSumOfDouble") +fun List.cumulativeSum() = this.cumulative(0.0) { element, sum -> sum + element } + +@JvmName("cumulativeSumOfInt") +fun List.cumulativeSum() = this.cumulative(0) { element, sum -> sum + element } + +@JvmName("cumulativeSumOfLong") +fun List.cumulativeSum() = this.cumulative(0L) { element, sum -> sum + element } \ No newline at end of file diff --git a/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/RingBuffer.kt b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/RingBuffer.kt new file mode 100644 index 000000000..108b7f828 --- /dev/null +++ b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/RingBuffer.kt @@ -0,0 +1,89 @@ +package scientifik.kmath.sequential + +import kotlinx.coroutines.sync.Mutex +import kotlinx.coroutines.sync.withLock +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.MutableBuffer +import scientifik.kmath.structures.VirtualBuffer + +/** + * Thread-safe ring buffer + */ +internal class RingBuffer( + private val buffer: MutableBuffer, + private var startIndex: Int = 0, + size: Int = 0 +) : Buffer { + + private val mutex = Mutex() + + override var size: Int = size + private set + + override fun get(index: Int): T { + require(index >= 0) { "Index must be positive" } + require(index < size) { "Index $index is out of circular buffer size $size" } + return buffer[startIndex.forward(index)] + } + + fun isFull() = size == buffer.size + + /** + * Iterator could provide wrong results if buffer is changed in initialization (iteration is safe) + */ + override fun iterator(): Iterator = object : AbstractIterator() { + private var count = size + private var index = startIndex + val copy = buffer.copy() + + override fun computeNext() { + if (count == 0) { + done() + } else { + setNext(copy[index]) + index = index.forward(1) + count-- + } + } + } + + /** + * A safe snapshot operation + */ + suspend fun snapshot(): Buffer { + mutex.withLock { + val copy = buffer.copy() + return VirtualBuffer(size) { i -> copy[startIndex.forward(i)] } + } + } + + suspend fun push(element: T) { + mutex.withLock { + buffer[startIndex.forward(size)] = element + if (isFull()) { + startIndex++ + } else { + size++ + } + } + } + + + @Suppress("NOTHING_TO_INLINE") + private inline fun Int.forward(n: Int): Int = (this + n) % (buffer.size) + + companion object { + inline fun build(size: Int, empty: T): RingBuffer { + val buffer = MutableBuffer.auto(size) { empty } + return RingBuffer(buffer) + } + + /** + * Slow yet universal buffer + */ + fun boxing(size: Int): RingBuffer { + val buffer: MutableBuffer = MutableBuffer.boxing(size) { null } + return RingBuffer(buffer) + } + } +} \ No newline at end of file diff --git a/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Streaming.kt b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Streaming.kt new file mode 100644 index 000000000..c0332d639 --- /dev/null +++ b/kmath-sequential/src/commonMain/kotlin/scientifik/kmath/sequential/Streaming.kt @@ -0,0 +1,275 @@ +package scientifik.kmath.sequential + +import kotlinx.coroutines.CoroutineScope +import kotlinx.coroutines.GlobalScope +import kotlinx.coroutines.channels.* +import kotlinx.coroutines.isActive +import kotlinx.coroutines.launch +import kotlinx.coroutines.sync.Mutex +import kotlinx.coroutines.sync.withLock +import scientifik.kmath.structures.Buffer +import kotlin.coroutines.CoroutineContext + +/** + * Initial chain block. Could produce an element sequence and be connected to single [Consumer] + * + * The general rule is that channel is created on first call. Also each element is responsible for its connection so + * while the connections are symmetric, the scope, used for making the connection is responsible for cancelation. + * + * Also connections are not reversible. Once connected block stays faithful until it finishes processing. + * Manually putting elements to connected block could lead to undetermined behavior and must be avoided. + */ +interface Producer : CoroutineScope { + fun connect(consumer: Consumer) + + suspend fun receive(): T + + val consumer: Consumer? + + val outputIsConnected: Boolean get() = consumer != null + + //fun close() +} + +/** + * Terminal chain block. Could consume an element sequence and be connected to signle [Producer] + */ +interface Consumer : CoroutineScope { + fun connect(producer: Producer) + + suspend fun send(value: T) + + val producer: Producer? + + val inputIsConnected: Boolean get() = producer != null + + //fun close() +} + +interface Processor : Consumer, Producer + +abstract class AbstractProducer(scope: CoroutineScope) : Producer { + override val coroutineContext: CoroutineContext = scope.coroutineContext + + override var consumer: Consumer? = null + protected set + + override fun connect(consumer: Consumer) { + //Ignore if already connected to specific consumer + if (consumer != this.consumer) { + if (outputIsConnected) error("The output slot of producer is occupied") + if (consumer.inputIsConnected) error("The input slot of consumer is occupied") + this.consumer = consumer + if (consumer.producer != null) { + //No need to save the job, it will be canceled on scope cancel + connectOutput(consumer) + // connect back, consumer is already set so no circular reference + consumer.connect(this) + } else error("Unreachable statement") + } + } + + protected open fun connectOutput(consumer: Consumer) { + launch { + while (this.isActive) { + consumer.send(receive()) + } + } + } +} + +abstract class AbstractConsumer(scope: CoroutineScope) : Consumer { + override val coroutineContext: CoroutineContext = scope.coroutineContext + + override var producer: Producer? = null + protected set + + override fun connect(producer: Producer) { + //Ignore if already connected to specific consumer + if (producer != this.producer) { + if (inputIsConnected) error("The input slot of consumer is occupied") + if (producer.outputIsConnected) error("The input slot of producer is occupied") + this.producer = producer + //No need to save the job, it will be canceled on scope cancel + if (producer.consumer != null) { + connectInput(producer) + // connect back + producer.connect(this) + } else error("Unreachable statement") + } + } + + protected open fun connectInput(producer: Producer) { + launch { + while (isActive) { + send(producer.receive()) + } + } + } +} + +abstract class AbstractProcessor(scope: CoroutineScope) : Processor, AbstractProducer(scope) { + + override var producer: Producer? = null + protected set + + override fun connect(producer: Producer) { + //Ignore if already connected to specific consumer + if (producer != this.producer) { + if (inputIsConnected) error("The input slot of consumer is occupied") + if (producer.outputIsConnected) error("The input slot of producer is occupied") + this.producer = producer + //No need to save the job, it will be canceled on scope cancel + if (producer.consumer != null) { + connectInput(producer) + // connect back + producer.connect(this) + } else error("Unreachable statement") + } + } + + protected open fun connectInput(producer: Producer) { + launch { + while (isActive) { + send(producer.receive()) + } + } + } +} + +/** + * A simple [produce]-based producer + */ +class GenericProducer( + scope: CoroutineScope, + capacity: Int = Channel.UNLIMITED, + block: suspend ProducerScope.() -> Unit +) : AbstractProducer(scope) { + + private val channel: ReceiveChannel by lazy { produce(capacity = capacity, block = block) } + + override suspend fun receive(): T = channel.receive() +} + +/** + * A simple pipeline [Processor] block + */ +class PipeProcessor( + scope: CoroutineScope, + capacity: Int = Channel.RENDEZVOUS, + process: suspend (T) -> R +) : AbstractProcessor(scope) { + + private val input = Channel(capacity) + private val output: ReceiveChannel = input.map(coroutineContext, process) + + override suspend fun receive(): R = output.receive() + + override suspend fun send(value: T) { + input.send(value) + } +} + + +/** + * A moving window [Processor] with circular buffer + */ +class WindowedProcessor( + scope: CoroutineScope, + window: Int, + val process: suspend (Buffer) -> R +) : AbstractProcessor(scope) { + + private val ringBuffer = RingBuffer.boxing(window) + + private val channel = Channel(Channel.RENDEZVOUS) + + override suspend fun receive(): R { + return channel.receive() + } + + override suspend fun send(value: T) { + ringBuffer.push(value) + channel.send(process(ringBuffer.snapshot())) + } +} + +/** + * Thread-safe aggregator of values from input. The aggregator does not store all incoming values, it uses fold procedure + * to incorporate them into state on-arrival. + * The current aggregated state could be accessed by [state]. The input channel is inactive unless requested + * @param T - the type of the input element + * @param S - the type of the aggregator + */ +class Reducer( + scope: CoroutineScope, + initialState: S, + val fold: suspend (S, T) -> S +) : AbstractConsumer(scope) { + + var state: S = initialState + private set + + private val mutex = Mutex() + + override suspend fun send(value: T) = mutex.withLock { + state = fold(state, value) + } +} + +/** + * Collector that accumulates all values in a list. List could be accessed from non-suspending environment via [list] value. + */ +class Collector(scope: CoroutineScope) : AbstractConsumer(scope) { + + private val _list = ArrayList() + private val mutex = Mutex() + val list: List get() = _list + + override suspend fun send(value: T) { + mutex.withLock { + _list.add(value) + } + } +} + +/** + * Convert a sequence to [Producer] + */ +fun Sequence.produce(scope: CoroutineScope = GlobalScope) = + GenericProducer(scope) { forEach { send(it) } } + +/** + * Convert a [ReceiveChannel] to [Producer] + */ +fun ReceiveChannel.produce(scope: CoroutineScope = GlobalScope) = + GenericProducer(scope) { for (e in this@produce) send(e) } + + +fun > Producer.consumer(consumerFactory: () -> C): C = + consumerFactory().also { connect(it) } + +fun Producer.map(capacity: Int = Channel.RENDEZVOUS, process: suspend (T) -> R) = + PipeProcessor(this, capacity, process).also { connect(it) } + +/** + * Create a reducer and connect this producer to reducer + */ +fun Producer.reduce(initialState: S, fold: suspend (S, T) -> S) = + Reducer(this, initialState, fold).also { connect(it) } + +/** + * Create a [Collector] and attach it to this [Producer] + */ +fun Producer.collect() = + Collector(this).also { connect(it) } + +fun > Producer.process(processorBuilder: () -> P): P = + processorBuilder().also { connect(it) } + +fun Producer.process(capacity: Int = Channel.RENDEZVOUS, process: suspend (T) -> R) = + PipeProcessor(this, capacity, process).also { connect(it) } + + +fun Producer.windowed(window: Int, process: suspend (Buffer) -> R) = + WindowedProcessor(this, window, process).also { connect(it) } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt b/kmath-sequential/src/commonTest/kotlin/scientifik/kmath/sequential/CumulativeKtTest.kt similarity index 86% rename from kmath-core/src/commonTest/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt rename to kmath-sequential/src/commonTest/kotlin/scientifik/kmath/sequential/CumulativeKtTest.kt index e7c99e7d0..cafa0526f 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt +++ b/kmath-sequential/src/commonTest/kotlin/scientifik/kmath/sequential/CumulativeKtTest.kt @@ -1,5 +1,6 @@ package scientifik.kmath.misc +import scientifik.kmath.sequential.cumulativeSum import kotlin.test.Test import kotlin.test.assertEquals diff --git a/kmath-sequential/src/jvmMain/kotlin/scientifik/kmath/sequential/ChainExt.kt b/kmath-sequential/src/jvmMain/kotlin/scientifik/kmath/sequential/ChainExt.kt new file mode 100644 index 000000000..74cb6bc6d --- /dev/null +++ b/kmath-sequential/src/jvmMain/kotlin/scientifik/kmath/sequential/ChainExt.kt @@ -0,0 +1,40 @@ +package scientifik.kmath.sequential + +import kotlinx.coroutines.runBlocking +import kotlin.sequences.Sequence + +/** + * Represent a chain as regular iterator (uses blocking calls) + */ +operator fun Chain.iterator() = object : Iterator { + override fun hasNext(): Boolean = true + + override fun next(): R = runBlocking { next() } +} + +/** + * Represent a chain as a sequence + */ +fun Chain.asSequence(): Sequence = object : Sequence { + override fun iterator(): Iterator = this@asSequence.iterator() +} + + +/** + * Map the chain result using suspended transformation. Initial chain result can no longer be safely consumed + * since mapped chain consumes tokens. Accepts suspending transformation function. + */ +fun Chain.map(func: suspend (T) -> R): Chain { + val parent = this; + return object : Chain { + override val value: R? get() = runBlocking { parent.value?.let { func(it) } } + + override suspend fun next(): R { + return func(parent.next()) + } + + override fun fork(): Chain { + return parent.fork().map(func) + } + } +} \ No newline at end of file diff --git a/kmath-sequential/src/jvmTest/kotlin/scientifik.kmath.sequential/RingBufferTest.kt b/kmath-sequential/src/jvmTest/kotlin/scientifik.kmath.sequential/RingBufferTest.kt new file mode 100644 index 000000000..c8f84e7d8 --- /dev/null +++ b/kmath-sequential/src/jvmTest/kotlin/scientifik.kmath.sequential/RingBufferTest.kt @@ -0,0 +1,19 @@ +package scientifik.kmath.sequential + +import kotlinx.coroutines.runBlocking +import scientifik.kmath.structures.asSequence +import kotlin.test.Test +import kotlin.test.assertEquals + +class RingBufferTest { + @Test + fun testPush() { + val buffer = RingBuffer.build(20, Double.NaN) + runBlocking { + for (i in 1..30) { + buffer.push(i.toDouble()) + } + assertEquals(410.0, buffer.asSequence().sum()) + } + } +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index 0c91b8446..8684de3ee 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -1,16 +1,31 @@ pluginManagement { repositories { mavenCentral() - maven("https://plugins.gradle.org/m2/") + gradlePluginPortal() + //maven ("https://dl.bintray.com/kotlin/kotlin-eap") + } + resolutionStrategy { + eachPlugin { + when (requested.id.id) { + "kotlinx-atomicfu" -> { + useModule("org.jetbrains.kotlinx:atomicfu-gradle-plugin:${requested.version}") + } + } + } } } -//enableFeaturePreview("GRADLE_METADATA") +enableFeaturePreview("GRADLE_METADATA") rootProject.name = "kmath" include( - ":kmath-core", - ":kmath-io", - ":kmath-coroutines", - ":kmath-jmh" + ":kmath-memory", + ":kmath-core", +// ":kmath-io", + ":kmath-coroutines", + ":kmath-histograms", + ":kmath-commons", + ":kmath-koma", + ":kmath-sequential", + ":benchmarks" )