pre-0.0.3 #46

Merged
altavir merged 75 commits from dev into master 2019-02-20 13:05:39 +03:00
118 changed files with 5777 additions and 1869 deletions

View File

@ -1,53 +1,54 @@
# KMath # KMath
Kotlin MATHematics library is intended as a kotlin based analog of numpy python library. Contrary to `numpy` 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.
and `scipy` it is modular and has a lightweight core.
## Features ## Features
* **Algebra** * **Algebra**
* Mathematical operation entities like rings, spaces and fields with (**TODO** add example to wiki) * Algebraic structures like rings, spaces and field (**TODO** add example to wiki)
* Basic linear algebra operations (summs products, etc) backed by `Space` API. * Basic linear algebra operations (sums, products, etc.), backed by the `Space` API.
* Complex numbers backed by `Field` API (meaning that they will be useable in any structures like vectors and NDArrays). * 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 inversions. * [In progress] advanced linear algebra operations like matrix inversion and LU decomposition.
* **Array-like structures** Full support of numpy-like ndarray including mixed arithmetic operations and function operations * **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).
on arrays and numbers just like it works in python (with 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 * **Expressions** Expressions are one of the ultimate goals of KMath. By writing a single mathematical expression
expression once an then apply it to different types of objects by providing different context. Expressions could be used once, users will be able to apply different types of objects to the expression by providing a context. Exceptions
for a wide variety of purposes from high performance calculations to code generation. can be used for a wide variety of purposes from high performance calculations to code generation.
## Planned features ## Planned features
* **Common mathematics** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/) * **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 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 future request if you want something to be done first. 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. * **Messaging** A mathematical notation to support multi-language and multi-node communication for mathematical tasks.
## Multi-platform support ## 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 KMath is developed as a multi-platform library, which means that most of interfaces are declared in the [common module](kmath-core/src/commonMain).
platform even if they could be done in common module because of platform performance optimization. Implementation is also done in the common module wherever possible. In some cases, features are delegated to
Currently the main focus of development is the JVM platform, contribution of implementations for Kotlin - Native and platform-specific implementations even if they could be done in the common module for performance reasons.
Kotlin - JS is welcome. Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome.
## Performance ## 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 Calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve
increasing performance for specific cases. We expect the worst KMath performance still be better than natural python, both performance and flexibility. We expect to focus on creating convenient universal API first and then work on
but worse than optimized native/scipy (mostly due to boxing operations on primitive numbers). The best performance increasing performance for specific cases. We expect the worst KMath benchmarks will perform better than native Python,
of optimized parts should be better than scipy. 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 ## 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 ```groovy
repositories { repositories {
maven { url = "http://npm.mipt.ru:8081/artifactory/gradle-dev" } maven { url = "http://npm.mipt.ru:8081/artifactory/gradle-dev" }
mavenCentral() mavenCentral()
} }
``` ```
or for kotlin gradle dsl:
or for the Gradle Kotlin DSL:
```kotlin ```kotlin
repositories { repositories {
@ -56,16 +57,20 @@ repositories {
} }
``` ```
Then use regular dependency like Then use a regular dependency like so:
```groovy ```groovy
compile(group: 'scientifik', name: 'kmath-core', version: '0.0.1-SNAPSHOT') compile(group: 'scientifik', name: 'kmath-core', version: '0.0.1-SNAPSHOT')
``` ```
or in kotlin
or in the Gradle Kotlin DSL:
```kotlin ```kotlin
compile(group = "scientifik", name = "kmath-core", version = "0.0.1-SNAPSHOT") 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 ## 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.

39
benchmarks/build.gradle Normal file
View File

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

View File

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

View File

@ -1,10 +1,10 @@
package scientifik.kmath.structures 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 import scientifik.kmath.operations.Complex
@Warmup(iterations = 1)
@Measurement(iterations = 5)
@State(Scope.Benchmark) @State(Scope.Benchmark)
open class BufferBenchmark { open class BufferBenchmark {
@ -22,17 +22,17 @@ open class BufferBenchmark {
@Benchmark @Benchmark
fun complexBufferReadWrite() { fun complexBufferReadWrite() {
val buffer = Complex.createBuffer(size/2) val buffer = MutableBuffer.complex(size / 2)
(0 until size/2).forEach { (0 until size / 2).forEach {
buffer[it] = Complex(it.toDouble(), -it.toDouble()) buffer[it] = Complex(it.toDouble(), -it.toDouble())
} }
(0 until size/2).forEach { (0 until size / 2).forEach {
buffer[it] buffer[it]
} }
} }
companion object { companion object {
const val size = 1000 const val size = 100
} }
} }

View File

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

View File

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

View File

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

View File

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

View File

@ -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<String>) {
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<Double> = 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<Double> = one
repeat(n) {
res += 1.0
}
}
}
println("Generic addition completed in $genericTime millis")
}

View File

@ -0,0 +1,36 @@
package scientifik.kmath.structures
import kotlin.system.measureTimeMillis
fun main(args: Array<String>) {
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")
}

View File

@ -0,0 +1,39 @@
package scientifik.kmath.structures
import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory
import kotlin.system.measureTimeMillis
fun main(args: Array<String>) {
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")
}

View File

@ -1,13 +1,13 @@
buildscript { import org.jetbrains.kotlin.gradle.dsl.KotlinMultiplatformExtension
extra["kotlinVersion"] = "1.3.11"
extra["ioVersion"] = "0.1.2-dev-2"
extra["coroutinesVersion"] = "1.0.1"
val kotlinVersion: String by extra buildscript {
val ioVersion: String by extra val kotlinVersion: String by rootProject.extra("1.3.21")
val coroutinesVersion: String by extra 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 { repositories {
//maven("https://dl.bintray.com/kotlin/kotlin-eap")
jcenter() jcenter()
} }
@ -18,18 +18,39 @@ buildscript {
} }
plugins { plugins {
id("com.jfrog.artifactory") version "4.8.1" apply false id("com.jfrog.artifactory") version "4.9.1" apply false
// id("org.jetbrains.kotlin.multiplatform") apply false
} }
allprojects { allprojects {
if (project.name.startsWith("kmath")) {
apply(plugin = "maven-publish") apply(plugin = "maven-publish")
apply(plugin = "com.jfrog.artifactory") apply(plugin = "com.jfrog.artifactory")
}
group = "scientifik" 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<KotlinMultiplatformExtension>()?.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") apply(from = "artifactory.gradle")
} }

70
doc/algebra.md Normal file
View File

@ -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<T>`. Next one needs to run actual operation
in the context:
```kotlin
val a: T
val b: T
val space: Space<T>
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.

1
doc/buffers.md Normal file
View File

@ -0,0 +1 @@
**TODO**

73
doc/contexts.md Normal file
View File

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

26
doc/expressions.md Normal file
View File

@ -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<String, T>): 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**

0
doc/features.md Normal file
View File

1
doc/histograms.md Normal file
View File

@ -0,0 +1 @@
**TODO**

1
doc/linear.md Normal file
View File

@ -0,0 +1 @@
**TODO**

127
doc/nd-performance.md Normal file
View File

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

40
doc/operations.md Normal file
View File

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

BIN
gradle/wrapper/gradle-wrapper.jar vendored Normal file

Binary file not shown.

View File

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

172
gradlew vendored Executable file
View File

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

84
gradlew.bat vendored Normal file
View File

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

View File

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

View File

@ -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<String, Double>) :
ExtendedField<DerivativeStructure> {
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<String, DerivativeStructure> = parameters.mapValues { (key, value) ->
DerivativeStructure(parameters.size, order, parameters.keys.indexOf(key), value)
}
val variable = object : ReadOnlyProperty<Any?, DerivativeStructure> {
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<String, Int>): Double {
return getPartialDerivative(*parameters.keys.map { orders[it] ?: 0 }.toIntArray())
}
fun DerivativeStructure.deriv(vararg orders: Pair<String, Int>): 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<Double> {
override fun invoke(arguments: Map<String, Double>): Double = DerivativeStructureField(0, arguments)
.run(function).value
/**
* Get the derivative expression with given orders
* TODO make result [DiffExpression]
*/
fun derivative(orders: Map<String, Int>): Expression<Double> {
return object : Expression<Double> {
override fun invoke(arguments: Map<String, Double>): 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<String, Int>) = 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<Double>, Field<DiffExpression> {
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) }
}

View File

@ -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<MatrixFeature>? = null) : Matrix<Double> {
override val rowNum: Int get() = origin.rowDimension
override val colNum: Int get() = origin.columnDimension
override val features: Set<MatrixFeature> = features ?: sequence<MatrixFeature> {
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<Double>.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<Double> {
override val size: Int get() = origin.dimension
override fun get(index: Int): Double = origin.getEntry(index)
override fun iterator(): Iterator<Double> = origin.toArray().iterator()
}
fun Point<Double>.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<Double>, LinearSolver<Double> {
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<Double>, b: Matrix<Double>): CMMatrix {
val decomposition = LUDecomposition(a.toCM().origin)
return decomposition.solver.solve(b.toCM().origin).toMatrix()
}
override fun solve(a: Matrix<Double>, b: Point<Double>): CMVector {
val decomposition = LUDecomposition(a.toCM().origin)
return decomposition.solver.solve(b.toCM().origin).toPoint()
}
override fun inverse(a: Matrix<Double>): CMMatrix {
val decomposition = LUDecomposition(a.toCM().origin)
return decomposition.solver.inverse.toMatrix()
}
override fun Matrix<Double>.dot(other: Matrix<Double>) =
CMMatrix(this.toCM().origin.multiply(other.toCM().origin))
override fun Matrix<Double>.dot(vector: Point<Double>): CMVector =
CMVector(this.toCM().origin.preMultiply(vector.toCM().origin))
override fun Matrix<Double>.unaryMinus(): CMMatrix =
produce(rowNum, colNum) { i, j -> -get(i, j) }
override fun Matrix<Double>.plus(b: Matrix<Double>) =
CMMatrix(this.toCM().origin.multiply(b.toCM().origin))
override fun Matrix<Double>.minus(b: Matrix<Double>) =
CMMatrix(this.toCM().origin.subtract(b.toCM().origin))
override fun Matrix<Double>.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))

View File

@ -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<Complex>.toArray(): Array<org.apache.commons.math3.complex.Complex> =
Array(size) { org.apache.commons.math3.complex.Complex(get(it).re, get(it).im) }
private fun Buffer<Double>.asArray() = if (this is DoubleBuffer) {
array
} else {
DoubleArray(size) { i -> get(i) }
}
/**
* Create a virtual buffer on top of array
*/
private fun Array<org.apache.commons.math3.complex.Complex>.asBuffer() = VirtualBuffer<Complex>(size) {
val value = get(it)
Complex(value.real, value.imaginary)
}
fun fourier(
normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD
): BufferTransform<Complex, Complex> = {
FastFourierTransformer(normalization).transform(it.toArray(), direction).asBuffer()
}
fun realFourier(
normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD
): BufferTransform<Double, Complex> = {
FastFourierTransformer(normalization).transform(it.asArray(), direction).asBuffer()
}
fun sine(
normalization: DstNormalization = DstNormalization.STANDARD_DST_I,
direction: TransformType = TransformType.FORWARD
): BufferTransform<Double, Double> = {
FastSineTransformer(normalization).transform(it.asArray(), direction).asBuffer()
}
fun cosine(
normalization: DctNormalization = DctNormalization.STANDARD_DCT_I,
direction: TransformType = TransformType.FORWARD
): BufferTransform<Double, Double> = {
FastCosineTransformer(normalization).transform(it.asArray(), direction).asBuffer()
}
fun hadamard(
direction: TransformType = TransformType.FORWARD
): BufferTransform<Double, Double> = {
FastHadamardTransformer().transform(it.asArray(), direction).asBuffer()
}
}
/**
* Process given [Producer] with commons-math fft transformation
*/
fun Producer<Buffer<Complex>>.FFT(
normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD
): Processor<Buffer<Complex>, Buffer<Complex>> {
val transform = Transformations.fourier(normalization, direction)
return map { transform(it) }
}
@JvmName("realFFT")
fun Producer<Buffer<Double>>.FFT(
normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD
): Processor<Buffer<Double>, Buffer<Complex>> {
val transform = Transformations.realFourier(normalization, direction)
return map { transform(it) }
}

View File

@ -0,0 +1,31 @@
package scientifik.kmath.expressions
import org.junit.Test
import kotlin.test.assertEquals
inline fun <R> diff(order: Int, vararg parameters: Pair<String, Double>, 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))
}
}

View File

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

View File

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

View File

@ -1,48 +1,64 @@
package scientifik.kmath.expressions package scientifik.kmath.expressions
import scientifik.kmath.operations.Field import scientifik.kmath.operations.Field
import scientifik.kmath.operations.Ring
import scientifik.kmath.operations.Space import scientifik.kmath.operations.Space
/**
* An elementary function that could be invoked on a map of arguments
*/
interface Expression<T> { interface Expression<T> {
operator fun invoke(arguments: Map<String, T>): T operator fun invoke(arguments: Map<String, T>): T
} }
operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs)) operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs))
/**
* A context for expression construction
*/
interface ExpressionContext<T> { interface ExpressionContext<T> {
/**
* Introduce a variable into expression context
*/
fun variable(name: String, default: T? = null): Expression<T> fun variable(name: String, default: T? = null): Expression<T>
/**
* A constant expression which does not depend on arguments
*/
fun const(value: T): Expression<T> fun const(value: T): Expression<T>
} }
internal class VariableExpression<T>(val name: String, val default: T? = null) : Expression<T> { internal class VariableExpression<T>(val name: String, val default: T? = null) : Expression<T> {
override fun invoke(arguments: Map<String, T>): T { override fun invoke(arguments: Map<String, T>): T =
return arguments[name] ?: default ?: error("The parameter not found: $name") arguments[name] ?: default ?: error("Parameter not found: $name")
}
} }
internal class ConstantExpression<T>(val value: T) : Expression<T> { internal class ConstantExpression<T>(val value: T) : Expression<T> {
override fun invoke(arguments: Map<String, T>): T = value override fun invoke(arguments: Map<String, T>): T = value
} }
internal class SumExpression<T>(val context: Space<T>, val first: Expression<T>, val second: Expression<T>) : Expression<T> { internal class SumExpression<T>(val context: Space<T>, val first: Expression<T>, val second: Expression<T>) :
Expression<T> {
override fun invoke(arguments: Map<String, T>): T = context.add(first.invoke(arguments), second.invoke(arguments)) override fun invoke(arguments: Map<String, T>): T = context.add(first.invoke(arguments), second.invoke(arguments))
} }
internal class ProductExpression<T>(val context: Field<T>, val first: Expression<T>, val second: Expression<T>) : Expression<T> { internal class ProductExpression<T>(val context: Ring<T>, val first: Expression<T>, val second: Expression<T>) :
override fun invoke(arguments: Map<String, T>): T = context.multiply(first.invoke(arguments), second.invoke(arguments)) Expression<T> {
override fun invoke(arguments: Map<String, T>): T =
context.multiply(first.invoke(arguments), second.invoke(arguments))
} }
internal class ConstProductExpession<T>(val context: Field<T>, val expr: Expression<T>, val const: Double) : Expression<T> { internal class ConstProductExpession<T>(val context: Space<T>, val expr: Expression<T>, val const: Number) :
Expression<T> {
override fun invoke(arguments: Map<String, T>): T = context.multiply(expr.invoke(arguments), const) override fun invoke(arguments: Map<String, T>): T = context.multiply(expr.invoke(arguments), const)
} }
internal class DivExpession<T>(val context: Field<T>, val expr: Expression<T>, val second: Expression<T>) : Expression<T> { internal class DivExpession<T>(val context: Field<T>, val expr: Expression<T>, val second: Expression<T>) :
Expression<T> {
override fun invoke(arguments: Map<String, T>): T = context.divide(expr.invoke(arguments), second.invoke(arguments)) override fun invoke(arguments: Map<String, T>): T = context.divide(expr.invoke(arguments), second.invoke(arguments))
} }
class FieldExpressionContext<T>(val field: Field<T>) : Field<Expression<T>>, ExpressionContext<T> { class ExpressionField<T>(val field: Field<T>) : Field<Expression<T>>, ExpressionContext<T> {
override val zero: Expression<T> = ConstantExpression(field.zero) override val zero: Expression<T> = ConstantExpression(field.zero)
@ -54,9 +70,20 @@ class FieldExpressionContext<T>(val field: Field<T>) : Field<Expression<T>>, Exp
override fun add(a: Expression<T>, b: Expression<T>): Expression<T> = SumExpression(field, a, b) override fun add(a: Expression<T>, b: Expression<T>): Expression<T> = SumExpression(field, a, b)
override fun multiply(a: Expression<T>, k: Double): Expression<T> = ConstProductExpession(field, a, k) override fun multiply(a: Expression<T>, k: Number): Expression<T> = ConstProductExpession(field, a, k)
override fun multiply(a: Expression<T>, b: Expression<T>): Expression<T> = ProductExpression(field, a, b) override fun multiply(a: Expression<T>, b: Expression<T>): Expression<T> = ProductExpression(field, a, b)
override fun divide(a: Expression<T>, b: Expression<T>): Expression<T> = DivExpession(field, a, b) override fun divide(a: Expression<T>, b: Expression<T>): Expression<T> = DivExpession(field, a, b)
operator fun Expression<T>.plus(arg: T) = this + const(arg)
operator fun Expression<T>.minus(arg: T) = this - const(arg)
operator fun Expression<T>.times(arg: T) = this * const(arg)
operator fun Expression<T>.div(arg: T) = this / const(arg)
operator fun T.plus(arg: Expression<T>) = arg + this
operator fun T.minus(arg: Expression<T>) = arg - this
operator fun T.times(arg: Expression<T>) = arg * this
operator fun T.div(arg: Expression<T>) = arg / this
} }

View File

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

View File

@ -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<T> = Buffer<T>
typealias RealPoint = Buffer<Double>
/**
* A simple geometric domain
* TODO move to geometry module
*/
interface Domain<T: Any> {
operator fun contains(vector: Point<out T>): Boolean
val dimension: Int
}
/**
* The bin in the histogram. The histogram is by definition always done in the real space
*/
interface Bin<T: Any> : Domain<T> {
/**
* The value of this bin
*/
val value: Number
val center: Point<T>
}
interface Histogram<T: Any, out B : Bin<T>> : Iterable<B> {
/**
* Find existing bin, corresponding to given coordinates
*/
operator fun get(point: Point<out T>): B?
/**
* Dimension of the histogram
*/
val dimension: Int
}
interface MutableHistogram<T: Any, out B : Bin<T>>: Histogram<T,B>{
/**
* Increment appropriate bin
*/
fun put(point: Point<out T>, weight: Double = 1.0)
}
fun <T: Any> MutableHistogram<T,*>.put(vararg point: T) = put(ArrayBuffer(point))
fun MutableHistogram<Double,*>.put(vararg point: Number) = put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray()))
fun MutableHistogram<Double,*>.put(vararg point: Double) = put(DoubleBuffer(point))
fun <T: Any> MutableHistogram<T,*>.fill(sequence: Iterable<Point<T>>) = sequence.forEach { put(it) }
/**
* Pass a sequence builder into histogram
*/
fun <T: Any> MutableHistogram<T, *>.fill(buider: suspend SequenceScope<Point<T>>.() -> Unit) = fill(sequence(buider).asIterable())

View File

@ -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<T : Comparable<T>>(val center: Vector<T, *>, val sizes: Point<T>) {
fun contains(vector: Point<out T>): 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<T : Any, B : Bin<T>, H : Histogram<T, B>> : Space<H> {
/**
* Rules for performing operations on bins
*/
val binSpace: Space<Bin<T>>
}
class PhantomBin<T : Comparable<T>>(val template: BinTemplate<T>, override val value: Number) : Bin<T> {
override fun contains(vector: Point<out T>): Boolean = template.contains(vector)
override val dimension: Int
get() = template.center.size
override val center: Point<T>
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<T : Comparable<T>>(
val bins: Map<BinTemplate<T>, IntArray>,
val data: NDStructure<Number>
) : Histogram<T, PhantomBin<T>> {
override val dimension: Int
get() = data.dimension
override fun iterator(): Iterator<PhantomBin<T>> {
return bins.asSequence().map { entry -> PhantomBin(entry.key, data[entry.value]) }.iterator()
}
override fun get(point: Point<out T>): PhantomBin<T>? {
val template = bins.keys.find { it.contains(point) }
return template?.let { PhantomBin(it, data[bins[it]!!]) }
}
}

View File

@ -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<T : Any, R : Ring<T>>(
override val elementContext: R,
private val bufferFactory: BufferFactory<T>
) : GenericMatrixContext<T, R> {
override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix<T> {
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<T> = bufferFactory(size, initializer)
}
class BufferMatrix<T : Any>(
override val rowNum: Int,
override val colNum: Int,
val buffer: Buffer<out T>,
override val features: Set<MatrixFeature> = emptySet()
) : Matrix<T> {
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<Pair<IntArray, T>> = 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<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
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)
}

View File

@ -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<T : Comparable<T>, F : Field<T>>(val matrix: Matrix<T, F>) {
private val field get() = matrix.context.field
/** Entries of LU decomposition. */
internal val lu: NDStructure<T>
/** 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<out T, F> 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<out T, F> 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<out T, F> 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 <T> MutableNDStructure<T>.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<NDStructure<T>, 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<T> = 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<Double, DoubleField>(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<Double, DoubleField> {
fun decompose(mat: Matrix<Double, DoubleField>, 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] }
}
}

View File

@ -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<T : Comparable<T>>(
private val elementContext: Ring<T>,
internal val lu: NDStructure<T>,
val pivot: IntArray,
private val even: Boolean
) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
/**
* Returns the matrix L of the decomposition.
*
* L is a lower-triangular matrix with [Ring.one] in diagonal
*/
override val l: Matrix<T> = 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<T> = 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<T> = 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<T : Comparable<T>, F : Field<T>>(
val context: GenericMatrixContext<T, F>,
val bufferFactory: MutableBufferFactory<T> = ::boxing,
val singularityCheck: (T) -> Boolean
) : LinearSolver<T> {
private fun abs(value: T) =
if (value > context.elementContext.zero) value else with(context.elementContext) { -value }
fun buildDecomposition(matrix: Matrix<T>): LUPDecomposition<T> {
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<T>): Matrix<T> {
if (matrix.hasFeature<LUPDecomposition<*>>()) {
return matrix
} else {
val decomposition = buildDecomposition(matrix)
return VirtualMatrix.wrap(matrix, decomposition)
}
}
override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
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<T>): Matrix<T> = solve(a, context.one(a.rowNum, a.colNum))
companion object {
val real = LUSolver(MatrixContext.real, MutableBuffer.Companion::auto) { it < 1e-11 }
}
}

View File

@ -1,62 +1,45 @@
package scientifik.kmath.linear package scientifik.kmath.linear
import scientifik.kmath.operations.DoubleField
import scientifik.kmath.operations.Field import scientifik.kmath.operations.Field
import scientifik.kmath.operations.Norm 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 * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
*/ */
interface LinearSolver<T : Any, F : Field<T>> { interface LinearSolver<T : Any> {
fun solve(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F> fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
fun solve(a: Matrix<T, F>, b: Vector<T, F>): Vector<T, F> = solve(a, b.toMatrix()).toVector() fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.toMatrix()).toVector()
fun inverse(a: Matrix<T, F>): Matrix<T, F> = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field)) fun inverse(a: Matrix<T>): Matrix<T>
} }
/** /**
* Convert vector to array (copying content of array) * Convert vector to array (copying content of array)
*/ */
fun <T : Any> Array<T>.toVector(field: Field<T>) = Vector.of(size, field) { this[it] } fun <T : Any> Array<T>.toVector(field: Field<T>) = Vector.generic(size, field) { this[it] }
fun DoubleArray.toVector() = Vector.real(this.size) { this[it] }
fun List<Double>.toVector() = Vector.real(this.size) { this[it] }
object VectorL2Norm : Norm<Point<out Number>, Double> {
override fun norm(arg: Point<out Number>): Double =
kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() })
}
typealias RealVector = Vector<Double, RealField>
typealias RealMatrix = Matrix<Double>
fun DoubleArray.toVector() = Vector.ofReal(this.size) { this[it] }
fun List<Double>.toVector() = Vector.ofReal(this.size) { this[it] }
/** /**
* Convert matrix to vector if it is possible * Convert matrix to vector if it is possible
*/ */
fun <T : Any, F : Field<T>> Matrix<T, F>.toVector(): Vector<T, F> { fun <T: Any> Matrix<T>.toVector(): Point<T> =
return when { if (this.colNum == 1) {
this.columns == 1 -> { VirtualBuffer(rowNum){ get(it, 0) }
// if (this is ArrayMatrix) { } else error("Can't convert matrix with more than one column to vector")
// //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 <T : Any, F : Field<T>> Vector<T, F>.toMatrix(): Matrix<T, F> { fun <T: Any> Point<T>.toMatrix(): Matrix<T> = VirtualMatrix(size, 1) { i, _ -> get(i) }
// 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<Vector<out Number, *>, Double> {
override fun norm(arg: Vector<out Number, *>): Double {
return kotlin.math.sqrt(arg.sumByDouble { it.toDouble() })
}
}
typealias RealVector = Vector<Double, DoubleField>
typealias RealMatrix = Matrix<Double, DoubleField>

View File

@ -1,187 +1,214 @@
package scientifik.kmath.linear 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.*
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. interface MatrixContext<T : Any> {
* @param T type of individual element of the vector or matrix /**
* @param V the type of vector space element * Produce a matrix with this context and given dimensions
*/ */
abstract class MatrixSpace<T : Any, F : Ring<T>>(val rows: Int, val columns: Int, val field: F) : Space<Matrix<T, F>> { fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
infix fun Matrix<T>.dot(vector: Point<T>): Point<T>
operator fun Matrix<T>.unaryMinus(): Matrix<T>
operator fun Matrix<T>.plus(b: Matrix<T>): Matrix<T>
operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T>
operator fun Matrix<T>.times(value: T): Matrix<T>
operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
companion object {
/**
* Non-boxing double matrix
*/
val real = BufferMatrixContext(RealField, DoubleBufferFactory)
/** /**
* Produce the element of this space * A structured matrix with custom buffer
*/ */
abstract fun produce(initializer: (Int, Int) -> T): Matrix<T, F> fun <T : Any, R : Ring<T>> buffered(
ring: R,
bufferFactory: BufferFactory<T> = ::boxing
): GenericMatrixContext<T, R> =
BufferMatrixContext(ring, bufferFactory)
/** /**
* Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements * Automatic buffered matrix, unboxed if it is possible
*/ */
abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace<T, F> inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
buffered(ring, Buffer.Companion::auto)
override val zero: Matrix<T, F> by lazy {
produce { _, _ -> field.zero }
}
// val one: Matrix<T> by lazy {
// produce { i, j -> if (i == j) field.one else field.zero }
// }
override fun add(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F> {
return produce { i, j -> with(field) { a[i, j] + b[i, j] } }
}
override fun multiply(a: Matrix<T, F>, k: Double): Matrix<T, F> {
//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 } }
}
/**
* Dot product. Throws exception on dimension mismatch
*/
fun multiply(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F> {
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) }
}
}
override fun equals(other: Any?): Boolean {
if (this === other) return true
if (other !is MatrixSpace<*,*>) return false
if (rows != other.rows) return false
if (columns != other.columns) return false
if (field != other.field) return false
return true
}
override fun hashCode(): Int {
var result = rows
result = 31 * result + columns
result = 31 * result + field.hashCode()
return result
} }
} }
infix fun <T : Any, F : Field<T>> Matrix<T, F>.dot(b: Matrix<T, F>): Matrix<T, F> = this.context.multiply(this, b) interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
/**
* The ring context for matrix elements
*/
val elementContext: R
/**
* Produce a point compatible with matrix space
*/
fun point(size: Int, initializer: (Int) -> T): Point<T>
override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
//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<T>.dot(vector: Point<T>): Point<T> {
//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<T>.unaryMinus() =
produce(rowNum, colNum) { i, j -> elementContext.run { -get(i, j) } }
override operator fun Matrix<T>.plus(b: Matrix<T>): Matrix<T> {
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<T>.minus(b: Matrix<T>): Matrix<T> {
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<T>.times(number: Number): Matrix<T> =
produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * number } }
operator fun Number.times(matrix: Matrix<T>): Matrix<T> = matrix * this
override fun Matrix<T>.times(value: T): Matrix<T> =
produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } }
}
/** /**
* A matrix-like structure * Specialized 2-d structure
*/ */
interface Matrix<T : Any, F: Ring<T>> : SpaceElement<Matrix<T, F>, MatrixSpace<T, F>> { interface Matrix<T : Any> : Structure2D<T> {
/** val rowNum: Int
* Number of rows val colNum: Int
*/
val rows: Int val features: Set<MatrixFeature>
/**
* Number of columns
*/
val columns: Int
/** /**
* 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<T>
override val self: Matrix<T, F> override fun get(index: IntArray): T = get(index[0], index[1])
get() = this
fun transpose(): Matrix<T, F> { override val shape: IntArray get() = intArrayOf(rowNum, colNum)
return object : Matrix<T, F> {
override val context: MatrixSpace<T, F> = this@Matrix.context val rows: Point<Point<T>>
override val rows: Int = this@Matrix.columns get() = VirtualBuffer(rowNum) { i ->
override val columns: Int = this@Matrix.rows VirtualBuffer(colNum) { j -> get(i, j) }
override fun get(i: Int, j: Int): T = this@Matrix[j, i] }
val columns: Point<Point<T>>
get() = VirtualBuffer(colNum) { j ->
VirtualBuffer(rowNum) { i -> get(i, j) }
}
override fun elements(): Sequence<Pair<IntArray, T>> = sequence {
for (i in (0 until rowNum)) {
for (j in (0 until colNum)) {
yield(intArrayOf(i, j) to get(i, j))
}
} }
} }
companion object { 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 <T : Any, F: Field<T>> of(rows: Int, columns: Int, field: F, initializer: (Int, Int) -> T) = fun <T : Any> square(vararg elements: T): Matrix<T> {
ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer) 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()
* Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing. return BufferMatrix(size, size, buffer)
*/
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 <T : Any, F: Field<T>> diagonal(rows: Int, columns: Int, field: F, values: (Int) -> T = { field.one }): Matrix<T, F> {
return of(rows, columns, field) { i, j -> if (i == j) values(i) else field.zero }
} }
/** fun <T : Any> build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)
* 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
}
} }
} }
class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
operator fun invoke(vararg elements: T): Matrix<T> {
if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns")
typealias NDFieldFactory<T, F> = (IntArray) -> NDField<T, F> val buffer = elements.asBuffer()
return BufferMatrix(rows, columns, buffer)
internal fun <T : Any, F : Field<T>> genericNDFieldFactory(field: F): NDFieldFactory<T, F> = { index -> GenericNDField(index, field) }
internal val realNDFieldFactory: NDFieldFactory<Double, DoubleField> = { 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<T : Any, F : Field<T>>(
rows: Int,
columns: Int,
field: F,
val ndFactory: NDFieldFactory<T, F> = genericNDFieldFactory(field)
) : MatrixSpace<T, F>(rows, columns, field) {
val ndField by lazy {
ndFactory(intArrayOf(rows, columns))
}
override fun produce(initializer: (Int, Int) -> T): Matrix<T, F> = ArrayMatrix(this, initializer)
override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace<T, F> {
return ArrayMatrixSpace(rows, columns, field, ndFactory)
} }
} }
/** /**
* Member of [ArrayMatrixSpace] which wraps 2-D array * Check if matrix has the given feature class
*/ */
class ArrayMatrix<T : Any, F : Field<T>> internal constructor(override val context: ArrayMatrixSpace<T, F>, val element: NDElement<T, F>) : Matrix<T, F> { inline fun <reified T : Any> Matrix<*>.hasFeature(): Boolean = features.find { it is T } != null
constructor(context: ArrayMatrixSpace<T, F>, 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 <reified T : Any> Matrix<*>.getFeature(): T? = features.filterIsInstance<T>().firstOrNull()
override val rows: Int get() = context.rows /**
* Diagonal matrix of ones. The matrix is virtual no actual matrix is created
override val columns: Int get() = context.columns */
fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, columns: Int): Matrix<T> =
override fun get(i: Int, j: Int): T { VirtualMatrix<T>(rows, columns) { i, j ->
return element[i, j] if (i == j) elementContext.one else elementContext.zero
} }
override val self: ArrayMatrix<T, F> get() = this
/**
* A virtual matrix of zeroes
*/
fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.zero(rows: Int, columns: Int): Matrix<T> =
VirtualMatrix<T>(rows, columns) { _, _ -> elementContext.zero }
class TransposedFeature<T : Any>(val original: Matrix<T>) : MatrixFeature
/**
* Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A`
*/
fun <T : Any, R : Ring<T>> Matrix<T>.transpose(): Matrix<T> {
return this.getFeature<TransposedFeature<T>>()?.original ?: VirtualMatrix(
this.colNum,
this.rowNum,
setOf(TransposedFeature(this))
) { i, j -> get(j, i) }
} }
infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = with(MatrixContext.real) { dot(other) }

View File

@ -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<T : Any> : MatrixFeature {
val inverse: Matrix<T>
}
/**
* A determinant container
*/
interface DeterminantFeature<T : Any> : MatrixFeature {
val determinant: T
}
@Suppress("FunctionName")
fun <T: Any> DeterminantFeature(determinant: T) = object: DeterminantFeature<T>{
override val determinant: T = determinant
}
/**
* Lower triangular matrix
*/
object LFeature: MatrixFeature
/**
* Upper triangular feature
*/
object UFeature: MatrixFeature
/**
* TODO add documentation
*/
interface LUPDecompositionFeature<T : Any> : MatrixFeature {
val l: Matrix<T>
val u: Matrix<T>
val p: Matrix<T>
}
//TODO add sparse matrix feature

View File

@ -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<T>(val rowNum: Int, val colNum: Int, val buffer: MutableBuffer<T>) : MutableNDStructure<T> {
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<Pair<IntArray, T>> = 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 <T> create(
rowNum: Int,
colNum: Int,
bufferFactory: MutableBufferFactory<T>,
init: (i: Int, j: Int) -> T
): Mutable2DStructure<T> {
val buffer = bufferFactory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
return Mutable2DStructure(rowNum, colNum, buffer)
}
}
}

View File

@ -1,103 +1,125 @@
package scientifik.kmath.linear package scientifik.kmath.linear
import scientifik.kmath.histogram.Point import scientifik.kmath.operations.RealField
import scientifik.kmath.operations.DoubleField
import scientifik.kmath.operations.Field
import scientifik.kmath.operations.Space import scientifik.kmath.operations.Space
import scientifik.kmath.operations.SpaceElement import scientifik.kmath.operations.SpaceElement
import scientifik.kmath.structures.NDElement import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.get import scientifik.kmath.structures.BufferFactory
import scientifik.kmath.structures.asSequence
typealias Point<T> = Buffer<T>
/** /**
* A linear space for vectors. * A linear space for vectors.
* Could be used on any point-like structure * Could be used on any point-like structure
*/ */
abstract class VectorSpace<T : Any, S : Space<T>>(val size: Int, val space: S) : Space<Point<T>> { interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
abstract fun produce(initializer: (Int) -> T): Vector<T, S> val size: Int
override val zero: Vector<T, S> by lazy { produce { space.zero } } val space: S
override fun add(a: Point<T>, b: Point<T>): Vector<T, S> = produce { with(space) { a[it] + b[it] } } fun produce(initializer: (Int) -> T): Point<T>
override fun multiply(a: Point<T>, k: Double): Vector<T, S> = produce { with(space) { a[it] * k } } /**
* Produce a space-element of this vector space for expressions
*/
fun produceElement(initializer: (Int) -> T): Vector<T, S>
override val zero: Point<T> get() = produce { space.zero }
override fun add(a: Point<T>, b: Point<T>): Point<T> = produce { with(space) { a[it] + b[it] } }
override fun multiply(a: Point<T>, k: Number): Point<T> = produce { with(space) { a[it] * k } }
//TODO add basis
companion object {
private val realSpaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
/**
* Non-boxing double vector space
*/
fun real(size: Int): BufferVectorSpace<Double, RealField> {
return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, RealField, Buffer.DoubleBufferFactory) }
}
/**
* A structured vector space with custom buffer
*/
fun <T : Any, S : Space<T>> buffered(
size: Int,
space: S,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
): VectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
/**
* Automatic buffered vector, unboxed if it is possible
*/
inline fun <reified T : Any, S : Space<T>> smart(size: Int, space: S): VectorSpace<T, S> =
buffered(size, space, Buffer.Companion::auto)
}
} }
/** /**
* A point coupled to the linear space * A point coupled to the linear space
*/ */
interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, VectorSpace<T, S>>, Point<T>, Iterable<T> { interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, Vector<T, S>, VectorSpace<T, S>>, Point<T> {
override val size: Int get() = context.size override val size: Int get() = context.size
override operator fun plus(b: Point<T>): Vector<T, S> = context.add(self, b) override operator fun plus(b: Point<T>): Vector<T, S> = context.add(this, b).wrap()
override operator fun minus(b: Point<T>): Vector<T, S> = context.add(self, context.multiply(b, -1.0)) override operator fun minus(b: Point<T>): Vector<T, S> = context.add(this, context.multiply(b, -1.0)).wrap()
override operator fun times(k: Number): Vector<T, S> = context.multiply(self, k.toDouble()) override operator fun times(k: Number): Vector<T, S> = context.multiply(this, k.toDouble()).wrap()
override operator fun div(k: Number): Vector<T, S> = context.multiply(self, 1.0 / k.toDouble()) override operator fun div(k: Number): Vector<T, S> = context.multiply(this, 1.0 / k.toDouble()).wrap()
companion object { companion object {
/** /**
* Create vector with custom field * Create vector with custom field
*/ */
fun <T : Any, F : Field<T>> of(size: Int, field: F, initializer: (Int) -> T) = fun <T : Any, S : Space<T>> generic(size: Int, field: S, initializer: (Int) -> T): Vector<T, S> =
ArrayVector(ArrayVectorSpace(size, field), initializer) VectorSpace.buffered(size, field).produceElement(initializer)
private val realSpaceCache = HashMap<Int, ArrayVectorSpace<Double, DoubleField>>() fun real(size: Int, initializer: (Int) -> Double): Vector<Double, RealField> =
VectorSpace.real(size).produceElement(initializer)
private fun getRealSpace(size: Int): ArrayVectorSpace<Double, DoubleField> { fun ofReal(vararg elements: Double): Vector<Double, RealField> =
return realSpaceCache.getOrPut(size){ArrayVectorSpace(size, DoubleField, realNDFieldFactory)} 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<T : Any, F : Field<T>>( data class BufferVectorSpace<T : Any, S : Space<T>>(
size: Int, override val size: Int,
field: F, override val space: S,
val ndFactory: NDFieldFactory<T, F> = genericNDFieldFactory(field) val bufferFactory: BufferFactory<T>
) : VectorSpace<T, F>(size, field) { ) : VectorSpace<T, S> {
val ndField by lazy { override fun produce(initializer: (Int) -> T) = bufferFactory(size, initializer)
ndFactory(intArrayOf(size)) override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
}
override fun produce(initializer: (Int) -> T): Vector<T, F> = ArrayVector(this, initializer)
} }
class ArrayVector<T : Any, F : Field<T>> internal constructor(override val context: VectorSpace<T, F>, val element: NDElement<T, F>) : Vector<T, F> { data class BufferVector<T : Any, S : Space<T>>(override val context: VectorSpace<T, S>, val buffer: Buffer<T>) :
Vector<T, S> {
constructor(context: ArrayVectorSpace<T, F>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
init { init {
if (context.size != element.shape[0]) { if (context.size != buffer.size) {
error("Array dimension mismatch") error("Array dimension mismatch")
} }
} }
override fun get(index: Int): T { override fun get(index: Int): T {
return element[index] return buffer[index]
} }
override val self: ArrayVector<T, F> get() = this override fun unwrap(): Point<T> = this
override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator() override fun Point<T>.wrap(): Vector<T, S> = BufferVector(context, this)
override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() } override fun iterator(): Iterator<T> = (0 until size).map { buffer[it] }.iterator()
override fun toString(): String =
this.asSequence().joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
} }

View File

@ -0,0 +1,47 @@
package scientifik.kmath.linear
class VirtualMatrix<T : Any>(
override val rowNum: Int,
override val colNum: Int,
override val features: Set<MatrixFeature> = emptySet(),
val generator: (i: Int, j: Int) -> T
) : Matrix<T> {
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 <T : Any> wrap(matrix: Matrix<T>, vararg features: MatrixFeature): Matrix<T> {
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]
}
}
}
}
}

View File

@ -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 <T, R> Iterator<T>.cumulative(initial: R, operation: (T, R) -> R): Iterator<R> = object : Iterator<R> {
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 <T, R> Iterable<T>.cumulative(initial: R, operation: (T, R) -> R): Iterable<R> = object : Iterable<R> {
override fun iterator(): Iterator<R> = this@cumulative.iterator().cumulative(initial, operation)
}
fun <T, R> Sequence<T>.cumulative(initial: R, operation: (T, R) -> R): Sequence<R> = object : Sequence<R> {
override fun iterator(): Iterator<R> = this@cumulative.iterator().cumulative(initial, operation)
}
fun <T, R> List<T>.cumulative(initial: R, operation: (T, R) -> R): List<R> = this.iterator().cumulative(initial, operation).asSequence().toList()
//Cumulative sum
@JvmName("cumulativeSumOfDouble")
fun Iterable<Double>.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element}
@JvmName("cumulativeSumOfInt")
fun Iterable<Int>.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element}
@JvmName("cumulativeSumOfLong")
fun Iterable<Long>.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element}
@JvmName("cumulativeSumOfDouble")
fun Sequence<Double>.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element}
@JvmName("cumulativeSumOfInt")
fun Sequence<Int>.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element}
@JvmName("cumulativeSumOfLong")
fun Sequence<Long>.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element}
@JvmName("cumulativeSumOfDouble")
fun List<Double>.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element}
@JvmName("cumulativeSumOfInt")
fun List<Int>.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element}
@JvmName("cumulativeSumOfLong")
fun List<Long>.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element}

View File

@ -8,8 +8,8 @@ package scientifik.kmath.misc
* *
* If step is negative, the same goes from upper boundary downwards * If step is negative, the same goes from upper boundary downwards
*/ */
fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double> { fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double> =
return when { when {
step == 0.0 -> error("Zero step in double progression") step == 0.0 -> error("Zero step in double progression")
step > 0 -> sequence { step > 0 -> sequence {
var current = start var current = start
@ -26,12 +26,11 @@ fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double>
} }
} }
} }
}
/** /**
* Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints] * Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints]
*/ */
fun ClosedFloatingPointRange<Double>.toGrid(numPoints: Int): DoubleArray { fun ClosedFloatingPointRange<Double>.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 } return DoubleArray(numPoints) { i -> start + (endInclusive - start) / (numPoints - 1) * i }
} }

View File

@ -1,37 +1,7 @@
package scientifik.kmath.operations package scientifik.kmath.operations
/** interface SpaceOperations<T> {
* 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<T, S> {
/**
* 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<T> {
/**
* Neutral element for sum operation
*/
val zero: T
/** /**
* Addition operation for two context elements * Addition operation for two context elements
*/ */
@ -40,7 +10,7 @@ interface Space<T> {
/** /**
* Multiplication operation for context element and real number * 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 //Operation to be performed in this context
operator fun T.unaryMinus(): T = multiply(this, -1.0) operator fun T.unaryMinus(): T = multiply(this, -1.0)
@ -50,69 +20,61 @@ interface Space<T> {
operator fun T.times(k: Number) = multiply(this, k.toDouble()) operator fun T.times(k: Number) = multiply(this, k.toDouble())
operator fun T.div(k: Number) = multiply(this, 1.0 / k.toDouble()) operator fun T.div(k: Number) = multiply(this, 1.0 / k.toDouble())
operator fun Number.times(b: T) = b * this 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<T, S : Space<T>> : MathElement<T, S> {
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<T> : Space<T> { interface Space<T> : SpaceOperations<T> {
/** /**
* neutral operation for multiplication * Neutral element for sum operation
*/ */
val one: T val zero: T
}
interface RingOperations<T> : SpaceOperations<T> {
/** /**
* Multiplication for two field elements * Multiplication for two field elements
*/ */
fun multiply(a: T, b: T): T fun multiply(a: T, b: T): T
operator fun T.times(b: T): T = multiply(this, b) operator fun T.times(b: T): T = multiply(this, b)
} }
/** /**
* Ring element * The same as {@link Space} but with additional multiplication operation
*/ */
interface RingElement<T, S : Ring<T>> : SpaceElement<T, S> { interface Ring<T> : Space<T>, RingOperations<T> {
override val context: S /**
* 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<T> : RingOperations<T> {
fun divide(a: T, b: T): T
operator fun T.div(b: T): T = divide(this, b)
} }
/** /**
* Four operations algebra * Four operations algebra
*/ */
interface Field<T> : Ring<T> { interface Field<T> : Ring<T>, FieldOperations<T> {
fun divide(a: T, b: T): T
operator fun T.div(b: T): T = divide(this, b)
operator fun Number.div(b: T) = this * divide(one, b) 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<T, S : Field<T>> : RingElement<T, S> {
override val context: S
operator fun div(b: T): T = context.divide(self, b)
} }

View File

@ -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<C> {
/**
* The context this element belongs to
*/
val context: C
}
interface MathWrapper<T, I> {
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<T, I : SpaceElement<T, I, S>, S : Space<T>> : MathElement<S>, MathWrapper<T, I> {
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<T, I : RingElement<T, I, R>, R : Ring<T>> : SpaceElement<T, I, R> {
operator fun times(b: T) = context.multiply(unwrap(), b).wrap()
}
/**
* Field element
*/
interface FieldElement<T, I : FieldElement<T, I, F>, F : Field<T>> : RingElement<T, I, F> {
override val context: F
operator fun div(b: T) = context.divide(unwrap(), b).wrap()
}

View File

@ -0,0 +1,7 @@
package scientifik.kmath.operations
import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.asSequence
fun <T> Space<T>.sum(data : Iterable<T>): T = data.fold(zero) { left, right -> add(left,right) }
fun <T> Space<T>.sum(data : Sequence<T>): T = data.fold(zero) { left, right -> add(left, right) }

View File

@ -1,45 +1,97 @@
package scientifik.kmath.operations 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 * A field for complex numbers
*/ */
object ComplexField : Field<Complex> { object ComplexField : ExtendedFieldOperations<Complex>, Field<Complex> {
override val zero: Complex = Complex(0.0, 0.0) 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 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 * Complex number class
*/ */
data class Complex(val re: Double, val im: Double) : FieldElement<Complex, ComplexField> { data class Complex(val re: Double, val im: Double) : FieldElement<Complex, Complex, ComplexField> {
override val self: Complex get() = this override fun unwrap(): Complex = this
override val context: ComplexField
get() = ComplexField override fun Complex.wrap(): Complex = this
override val context: ComplexField get() = ComplexField
/** /**
* A complex conjugate * A complex conjugate
*/ */
val conjugate: Complex val conjugate: Complex get() = Complex(re, -im)
get() = Complex(re, -im)
val square: Double val square: Double get() = re * re + im * im
get() = re * re + im * im
val abs: Double val abs: Double get() = sqrt(square)
get() = kotlin.math.sqrt(square)
companion object { val theta: Double get() = atan(im / re)
companion object : MemorySpec<Complex> {
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)
}
} }
}
fun Double.toComplex() = Complex(this, 0.0)
fun Buffer.Companion.complex(size: Int, init: (Int) -> Complex): Buffer<Complex> {
return ObjectBuffer.create(Complex, size, init)
}
fun MutableBuffer.Companion.complex(size: Int, init: (Int) -> Complex): Buffer<Complex> {
return ObjectBuffer.create(Complex, size, init)
} }

View File

@ -1,89 +0,0 @@
package scientifik.kmath.operations
import kotlin.math.pow
/**
* Advanced Number-like field that implements basic operations
*/
interface ExtendedField<N : Any> :
Field<N>,
TrigonometricOperations<N>,
PowerOperations<N>,
ExponentialOperations<N>
/**
* Field for real values
*/
object RealField : ExtendedField<Real>, Norm<Real, Real> {
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<Real, RealField> {
//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<Double>, Norm<Double, Double> {
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<Int>{
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
}

View File

@ -0,0 +1,109 @@
package scientifik.kmath.operations
import kotlin.math.pow
/**
* Advanced Number-like field that implements basic operations
*/
interface ExtendedFieldOperations<T> :
FieldOperations<T>,
TrigonometricOperations<T>,
PowerOperations<T>,
ExponentialOperations<T>
interface ExtendedField<T> : ExtendedFieldOperations<T>, Field<T>
/**
* 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<Double, Real, RealField> {
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<Double>, ExtendedFieldOperations<Double>, Norm<Double, Double> {
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<Int>, Norm<Int, Int> {
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<Short>, Norm<Short, Short> {
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<Byte>, Norm<Byte, Byte> {
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<Long>, Norm<Long, Long> {
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
}

View File

@ -10,7 +10,7 @@ package scientifik.kmath.operations
* It also allows to override behavior for optional operations * It also allows to override behavior for optional operations
* *
*/ */
interface TrigonometricOperations<T> : Field<T> { interface TrigonometricOperations<T> : FieldOperations<T> {
fun sin(arg: T): T fun sin(arg: T): T
fun cos(arg: T): T fun cos(arg: T): T
@ -19,10 +19,10 @@ interface TrigonometricOperations<T> : Field<T> {
fun ctg(arg: T): T = cos(arg) / sin(arg) fun ctg(arg: T): T = cos(arg) / sin(arg)
} }
fun <T : FieldElement<T, out TrigonometricOperations<T>>> sin(arg: T): T = arg.context.sin(arg) fun <T : MathElement<out TrigonometricOperations<T>>> sin(arg: T): T = arg.context.sin(arg)
fun <T : FieldElement<T, out TrigonometricOperations<T>>> cos(arg: T): T = arg.context.cos(arg) fun <T : MathElement<out TrigonometricOperations<T>>> cos(arg: T): T = arg.context.cos(arg)
fun <T : FieldElement<T, out TrigonometricOperations<T>>> tg(arg: T): T = arg.context.tg(arg) fun <T : MathElement<out TrigonometricOperations<T>>> tg(arg: T): T = arg.context.tg(arg)
fun <T : FieldElement<T, out TrigonometricOperations<T>>> ctg(arg: T): T = arg.context.ctg(arg) fun <T : MathElement<out TrigonometricOperations<T>>> ctg(arg: T): T = arg.context.ctg(arg)
/* Power and roots */ /* Power and roots */
@ -30,12 +30,13 @@ fun <T : FieldElement<T, out TrigonometricOperations<T>>> ctg(arg: T): T = arg.c
* A context extension to include power operations like square roots, etc * A context extension to include power operations like square roots, etc
*/ */
interface PowerOperations<T> { interface PowerOperations<T> {
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 : MathElement<T, out PowerOperations<T>>> T.pow(power: Double): T = context.power(this, power) infix fun <T : MathElement<out PowerOperations<T>>> T.pow(power: Double): T = context.power(this, power)
fun <T : MathElement<T, out PowerOperations<T>>> sqrt(arg: T): T = arg pow 0.5 fun <T : MathElement<out PowerOperations<T>>> sqrt(arg: T): T = arg pow 0.5
fun <T : MathElement<T, out PowerOperations<T>>> sqr(arg: T): T = arg pow 2.0 fun <T : MathElement<out PowerOperations<T>>> sqr(arg: T): T = arg pow 2.0
/* Exponential */ /* Exponential */
@ -44,11 +45,11 @@ interface ExponentialOperations<T> {
fun ln(arg: T): T fun ln(arg: T): T
} }
fun <T : MathElement<T, out ExponentialOperations<T>>> exp(arg: T): T = arg.context.exp(arg) fun <T : MathElement<out ExponentialOperations<T>>> exp(arg: T): T = arg.context.exp(arg)
fun <T : MathElement<T, out ExponentialOperations<T>>> ln(arg: T): T = arg.context.ln(arg) fun <T : MathElement<out ExponentialOperations<T>>> ln(arg: T): T = arg.context.ln(arg)
interface Norm<in T, out R> { interface Norm<in T: Any, out R> {
fun norm(arg: T): R fun norm(arg: T): R
} }
fun <T : MathElement<T, out Norm<T, R>>, R> norm(arg: T): R = arg.context.norm(arg) fun <T : MathElement<out Norm<T, R>>, R> norm(arg: T): R = arg.context.norm(arg)

View File

@ -0,0 +1,73 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.Field
import scientifik.kmath.operations.FieldElement
class BoxingNDField<T, F : Field<T>>(
override val shape: IntArray,
override val elementContext: F,
val bufferFactory: BufferFactory<T>
) : BufferedNDField<T, F> {
override val strides: Strides = DefaultStrides(shape)
fun buildBuffer(size: Int, initializer: (Int) -> T): Buffer<T> =
bufferFactory(size, initializer)
override fun check(vararg elements: NDBuffer<T>) {
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<T>, transform: F.(T) -> T): BufferedNDFieldElement<T, F> {
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<T>,
transform: F.(index: IntArray, T) -> T
): BufferedNDFieldElement<T, F> {
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<T>,
b: NDBuffer<T>,
transform: F.(T, T) -> T
): BufferedNDFieldElement<T, F> {
check(a, b)
return BufferedNDFieldElement(
this,
buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) })
}
override fun NDBuffer<T>.toElement(): FieldElement<NDBuffer<T>, *, out BufferedNDField<T, F>> =
BufferedNDFieldElement(this@BoxingNDField, buffer)
}

View File

@ -0,0 +1,43 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.*
interface BufferedNDAlgebra<T, C>: NDAlgebra<T, C, NDBuffer<T>>{
val strides: Strides
override fun check(vararg elements: NDBuffer<T>) {
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<T>.toBuffer(): NDBuffer<T> {
return if (this is NDBuffer<T> && this.strides == this@BufferedNDAlgebra.strides) {
this
} else {
produce { index -> get(index) }
}
}
/**
* Convert a buffer to element of this algebra
*/
fun NDBuffer<T>.toElement(): MathElement<out BufferedNDAlgebra<T, C>>
}
interface BufferedNDSpace<T, S : Space<T>> : NDSpace<T, S, NDBuffer<T>>, BufferedNDAlgebra<T,S> {
override fun NDBuffer<T>.toElement(): SpaceElement<NDBuffer<T>, *, out BufferedNDSpace<T, S>>
}
interface BufferedNDRing<T, R : Ring<T>> : NDRing<T, R, NDBuffer<T>>, BufferedNDSpace<T, R> {
override fun NDBuffer<T>.toElement(): RingElement<NDBuffer<T>, *, out BufferedNDRing<T, R>>
}
interface BufferedNDField<T, F : Field<T>> : NDField<T, F, NDBuffer<T>>, BufferedNDRing<T, F> {
override fun NDBuffer<T>.toElement(): FieldElement<NDBuffer<T>, *, out BufferedNDField<T, F>>
}

View File

@ -0,0 +1,88 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.*
/**
* Base interface for an element with context, containing strides
*/
interface BufferedNDElement<T, C> : NDBuffer<T>, NDElement<T, C, NDBuffer<T>> {
override val context: BufferedNDAlgebra<T, C>
override val strides get() = context.strides
override val shape: IntArray get() = context.shape
}
class BufferedNDSpaceElement<T, S : Space<T>>(
override val context: BufferedNDSpace<T, S>,
override val buffer: Buffer<T>
) : BufferedNDElement<T, S>, SpaceElement<NDBuffer<T>, BufferedNDSpaceElement<T, S>, BufferedNDSpace<T, S>> {
override fun unwrap(): NDBuffer<T> = this
override fun NDBuffer<T>.wrap(): BufferedNDSpaceElement<T, S> {
context.check(this)
return BufferedNDSpaceElement(context, buffer)
}
}
class BufferedNDRingElement<T, R : Ring<T>>(
override val context: BufferedNDRing<T, R>,
override val buffer: Buffer<T>
) : BufferedNDElement<T, R>, RingElement<NDBuffer<T>, BufferedNDRingElement<T, R>, BufferedNDRing<T, R>> {
override fun unwrap(): NDBuffer<T> = this
override fun NDBuffer<T>.wrap(): BufferedNDRingElement<T, R> {
context.check(this)
return BufferedNDRingElement(context, buffer)
}
}
class BufferedNDFieldElement<T, F : Field<T>>(
override val context: BufferedNDField<T, F>,
override val buffer: Buffer<T>
) : BufferedNDElement<T, F>, FieldElement<NDBuffer<T>, BufferedNDFieldElement<T, F>, BufferedNDField<T, F>> {
override fun unwrap(): NDBuffer<T> = this
override fun NDBuffer<T>.wrap(): BufferedNDFieldElement<T, F> {
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 <T : Any, F : Field<T>> Function1<T, T>.invoke(ndElement: BufferedNDElement<T, F>) =
ndElement.context.run { map(ndElement) { invoke(it) }.toElement() }
/* plus and minus */
/**
* Summation operation for [BufferedNDElement] and single element
*/
operator fun <T : Any, F : Space<T>> BufferedNDElement<T, F>.plus(arg: T) =
context.map(this) { it + arg }.wrap()
/**
* Subtraction operation between [BufferedNDElement] and single element
*/
operator fun <T : Any, F : Space<T>> BufferedNDElement<T, F>.minus(arg: T) =
context.map(this) { it - arg }.wrap()
/* prod and div */
/**
* Product operation for [BufferedNDElement] and single element
*/
operator fun <T : Any, F : Ring<T>> BufferedNDElement<T, F>.times(arg: T) =
context.map(this) { it * arg }.wrap()
/**
* Division operation between [BufferedNDElement] and single element
*/
operator fun <T : Any, F : Field<T>> BufferedNDElement<T, F>.div(arg: T) =
context.map(this) { it / arg }.wrap()

View File

@ -1,20 +1,70 @@
package scientifik.kmath.structures package scientifik.kmath.structures
typealias BufferFactory<T> = (Int, (Int) -> T) -> Buffer<T>
typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
/** /**
* A generic random access structure for both primitives and objects * A generic random access structure for both primitives and objects
*/ */
interface Buffer<T> { interface Buffer<T> {
/**
* The size of the buffer
*/
val size: Int val size: Int
/**
* Get element at given index
*/
operator fun get(index: Int): T operator fun get(index: Int): T
/**
* Iterate over all elements
*/
operator fun iterator(): Iterator<T> operator fun iterator(): Iterator<T>
/**
* 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 <T> boxing(size: Int, initializer: (Int) -> T): Buffer<T> = ListBuffer(List(size, initializer))
/**
* Create most appropriate immutable buffer for given type avoiding boxing wherever possible
*/
@Suppress("UNCHECKED_CAST")
inline fun <reified T : Any> auto(size: Int, crossinline initializer: (Int) -> T): Buffer<T> {
return when (T::class) {
Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer<T>
Short::class -> ShortBuffer(ShortArray(size) { initializer(it) as Short }) as Buffer<T>
Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer<T>
Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as Buffer<T>
else -> boxing(size, initializer)
}
}
val DoubleBufferFactory: BufferFactory<Double> =
{ size, initializer -> DoubleBuffer(DoubleArray(size, initializer)) }
val ShortBufferFactory: BufferFactory<Short> =
{ size, initializer -> ShortBuffer(ShortArray(size, initializer)) }
val IntBufferFactory: BufferFactory<Int> = { size, initializer -> IntBuffer(IntArray(size, initializer)) }
val LongBufferFactory: BufferFactory<Long> = { size, initializer -> LongBuffer(LongArray(size, initializer)) }
}
} }
fun <T> Buffer<T>.asSequence(): Sequence<T> = iterator().asSequence() fun <T> Buffer<T>.asSequence(): Sequence<T> = iterator().asSequence()
fun <T> Buffer<T>.asIterable(): Iterable<T> = iterator().asSequence().asIterable()
interface MutableBuffer<T> : Buffer<T> { interface MutableBuffer<T> : Buffer<T> {
operator fun set(index: Int, value: T) operator fun set(index: Int, value: T)
@ -22,10 +72,32 @@ interface MutableBuffer<T> : Buffer<T> {
* A shallow copy of the buffer * A shallow copy of the buffer
*/ */
fun copy(): MutableBuffer<T> fun copy(): MutableBuffer<T>
companion object {
/**
* Create a boxing mutable buffer of given type
*/
inline fun <T> boxing(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
MutableListBuffer(MutableList(size, initializer))
/**
* Create most appropriate mutable buffer for given type avoiding boxing wherever possible
*/
@Suppress("UNCHECKED_CAST")
inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> {
return when (T::class) {
Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer<T>
Short::class -> ShortBuffer(ShortArray(size) { initializer(it) as Short }) as MutableBuffer<T>
Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer<T>
Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as MutableBuffer<T>
else -> boxing(size, initializer)
}
}
}
} }
inline class ListBuffer<T>(private val list: List<T>) : Buffer<T> { inline class ListBuffer<T>(val list: List<T>) : Buffer<T> {
override val size: Int override val size: Int
get() = list.size get() = list.size
@ -35,7 +107,12 @@ inline class ListBuffer<T>(private val list: List<T>) : Buffer<T> {
override fun iterator(): Iterator<T> = list.iterator() override fun iterator(): Iterator<T> = list.iterator()
} }
inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T> { fun <T> List<T>.asBuffer() = ListBuffer(this)
@Suppress("FunctionName")
inline fun <T> ListBuffer(size: Int, init: (Int) -> T) = List(size, init).asBuffer()
inline class MutableListBuffer<T>(val list: MutableList<T>) : MutableBuffer<T> {
override val size: Int override val size: Int
get() = list.size get() = list.size
@ -47,12 +124,11 @@ inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuf
} }
override fun iterator(): Iterator<T> = list.iterator() override fun iterator(): Iterator<T> = list.iterator()
override fun copy(): MutableBuffer<T> = MutableListBuffer(ArrayList(list)) override fun copy(): MutableBuffer<T> = MutableListBuffer(ArrayList(list))
} }
class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> { class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
//Can't inline because array invariant //Can't inline because array is invariant
override val size: Int override val size: Int
get() = array.size get() = array.size
@ -67,9 +143,10 @@ class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf()) override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
} }
inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> { fun <T> Array<T>.asBuffer() = ArrayBuffer(this)
override val size: Int
get() = array.size inline class DoubleBuffer(val array: DoubleArray) : MutableBuffer<Double> {
override val size: Int get() = array.size
override fun get(index: Int): Double = array[index] override fun get(index: Int): Double = array[index]
@ -77,14 +154,45 @@ inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double
array[index] = value array[index] = value
} }
override fun iterator(): Iterator<Double> = array.iterator() override fun iterator() = array.iterator()
override fun copy(): MutableBuffer<Double> = DoubleBuffer(array.copyOf()) override fun copy(): MutableBuffer<Double> = DoubleBuffer(array.copyOf())
} }
inline class IntBuffer(private val array: IntArray) : MutableBuffer<Int> { @Suppress("FunctionName")
override val size: Int inline fun DoubleBuffer(size: Int, init: (Int) -> Double) = DoubleBuffer(DoubleArray(size) { init(it) })
get() = array.size
/**
* Transform buffer of doubles into array for high performance operations
*/
val Buffer<out Double>.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<Short> {
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<Short> = ShortBuffer(array.copyOf())
}
fun ShortArray.asBuffer() = ShortBuffer(this)
inline class IntBuffer(val array: IntArray) : MutableBuffer<Int> {
override val size: Int get() = array.size
override fun get(index: Int): Int = array[index] override fun get(index: Int): Int = array[index]
@ -92,17 +200,55 @@ inline class IntBuffer(private val array: IntArray) : MutableBuffer<Int> {
array[index] = value array[index] = value
} }
override fun iterator(): Iterator<Int> = array.iterator() override fun iterator() = array.iterator()
override fun copy(): MutableBuffer<Int> = IntBuffer(array.copyOf()) override fun copy(): MutableBuffer<Int> = IntBuffer(array.copyOf())
} }
inline class ReadOnlyBuffer<T>(private val buffer: MutableBuffer<T>) : Buffer<T> { fun IntArray.asBuffer() = IntBuffer(this)
inline class LongBuffer(val array: LongArray) : MutableBuffer<Long> {
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<Long> = LongBuffer(array.copyOf())
}
fun LongArray.asBuffer() = LongBuffer(this)
inline class ReadOnlyBuffer<T>(val buffer: MutableBuffer<T>) : Buffer<T> {
override val size: Int get() = buffer.size override val size: Int get() = buffer.size
override fun get(index: Int): T = buffer.get(index) override fun get(index: Int): T = buffer.get(index)
override fun iterator(): Iterator<T> = 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<T>(override val size: Int, private val generator: (Int) -> T) : Buffer<T> {
override fun get(index: Int): T = generator(index)
override fun iterator(): Iterator<T> = (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 <T> Buffer<T>.asReadOnly(): Buffer<T> = if (this is MutableBuffer) {
} }
/** /**
* Create most appropriate immutable buffer for given type avoiding boxing wherever possible * Typealias for buffer transformations
*/ */
@Suppress("UNCHECKED_CAST") typealias BufferTransform<T, R> = (Buffer<T>) -> Buffer<R>
inline fun <reified T : Any> buffer(size: Int, noinline initializer: (Int) -> T): Buffer<T> {
return when (T::class) {
Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer<T>
Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer<T>
else -> ArrayBuffer(Array(size, initializer))
}
}
/**
* Create most appropriate mutable buffer for given type avoiding boxing wherever possible
*/
@Suppress("UNCHECKED_CAST")
inline fun <reified T : Any> mutableBuffer(size: Int, noinline initializer: (Int) -> T): MutableBuffer<T> {
return when (T::class) {
Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer<T>
Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer<T>
else -> ArrayBuffer(Array(size, initializer))
}
}

View File

@ -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<Complex, ComplexField>
/**
* An optimized nd-field for complex numbers
*/
class ComplexNDField(override val shape: IntArray) :
BufferedNDField<Complex, ComplexField>,
ExtendedNDField<Complex, ComplexField, NDBuffer<Complex>> {
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<Complex> =
Buffer.complex(size) { initializer(it) }
/**
* Inline transform an NDStructure to another structure
*/
override fun map(
arg: NDBuffer<Complex>,
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<Complex>,
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<Complex>,
b: NDBuffer<Complex>,
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<Complex>.toElement(): FieldElement<NDBuffer<Complex>, *, out BufferedNDField<Complex, ComplexField>> =
BufferedNDFieldElement(this@ComplexNDField, buffer)
override fun power(arg: NDBuffer<Complex>, pow: Number) = map(arg) { power(it, pow) }
override fun exp(arg: NDBuffer<Complex>) = map(arg) { exp(it) }
override fun ln(arg: NDBuffer<Complex>) = map(arg) { ln(it) }
override fun sin(arg: NDBuffer<Complex>) = map(arg) { sin(it) }
override fun cos(arg: NDBuffer<Complex>) = map(arg) { cos(it) }
}
/**
* Fast element production using function inlining
*/
inline fun BufferedNDField<Complex, ComplexField>.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<Complex, Complex>.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)

View File

@ -1,42 +1,46 @@
package scientifik.kmath.structures package scientifik.kmath.structures
import scientifik.kmath.operations.ExponentialOperations import scientifik.kmath.operations.*
import scientifik.kmath.operations.ExtendedField
import scientifik.kmath.operations.PowerOperations
import scientifik.kmath.operations.TrigonometricOperations
/** interface ExtendedNDField<T : Any, F, N : NDStructure<T>> :
* NDField that supports [ExtendedField] operations on its elements NDField<T, F, N>,
*/ TrigonometricOperations<N>,
class ExtendedNDField<N : Any, F : ExtendedField<N>>(shape: IntArray, field: F) : NDField<N, F>(shape, field), PowerOperations<N>,
TrigonometricOperations<NDElement<N, F>>, ExponentialOperations<N>
PowerOperations<NDElement<N, F>>, where F : ExtendedFieldOperations<T>, F : Field<T>
ExponentialOperations<NDElement<N, F>> {
override fun produceStructure(initializer: F.(IntArray) -> N): NDStructure<N> { ///**
return genericNdStructure(shape) { field.initializer(it) } // * NDField that supports [ExtendedField] operations on its elements
} // */
//class ExtendedNDFieldWrapper<T : Any, F : ExtendedField<T>, N : NDStructure<T>>(private val ndField: NDField<T, F, N>) :
override fun power(arg: NDElement<N, F>, pow: Double): NDElement<N, F> { // ExtendedNDField<T, F, N>, NDField<T, F, N> by ndField {
return arg.transform { d -> with(field) { power(d, pow) } } //
} // override val shape: IntArray get() = ndField.shape
// override val elementContext: F get() = ndField.elementContext
override fun exp(arg: NDElement<N, F>): NDElement<N, F> { //
return arg.transform { d -> with(field) { exp(d) } } // override fun produce(initializer: F.(IntArray) -> T) = ndField.produce(initializer)
} //
// override fun power(arg: N, pow: Double): N {
override fun ln(arg: NDElement<N, F>): NDElement<N, F> { // return produce { with(elementContext) { power(arg[it], pow) } }
return arg.transform { d -> with(field) { ln(d) } } // }
} //
// override fun exp(arg: N): N {
override fun sin(arg: NDElement<N, F>): NDElement<N, F> { // return produce { with(elementContext) { exp(arg[it]) } }
return arg.transform { d -> with(field) { sin(d) } } // }
} //
// override fun ln(arg: N): N {
override fun cos(arg: NDElement<N, F>): NDElement<N, F> { // return produce { with(elementContext) { ln(arg[it]) } }
return arg.transform { d -> with(field) { cos(d) } } // }
} //
} // 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]) } }
// }
//}

View File

@ -1,20 +0,0 @@
package scientifik.kmath.structures
//
//class LazyStructureField<T: Any>(val field: Field<T>): Field<LazyStructure<T>>{
//
//}
//
//class LazyStructure<T : Any> : NDStructure<T> {
//
// 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<Pair<IntArray, T>> {
// TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
// }
//}

View File

@ -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<T, C, N : NDStructure<T>> {
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<T, T>.invoke(structure: N) = map(structure) { value -> this@invoke(value) }
}
/**
* An nd-space over element space
*/
interface NDSpace<T, S : Space<T>, N : NDStructure<T>> : Space<N>, NDAlgebra<T, S, N> {
/**
* 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<T, R : Ring<T>, N : NDStructure<T>> : Ring<N>, NDSpace<T, R, N> {
/**
* 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<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F, N> {
/**
* 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<IntArray, RealNDField>()
/**
* 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 <T : Any, F : Field<T>> buffered(
shape: IntArray,
field: F,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
) =
BoxingNDField(shape, field, bufferFactory)
/**
* Create a most suitable implementation for nd-field using reified class.
*/
@Suppress("UNCHECKED_CAST")
inline fun <reified T : Any, F : Field<T>> auto(field: F, vararg shape: Int): BufferedNDField<T, F> =
when {
T::class == Double::class -> real(*shape) as BufferedNDField<T, F>
else -> BoxingNDField(shape, field, Buffer.Companion::auto)
}
}
}

View File

@ -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<T, C, N : NDStructure<T>> : NDStructure<T> {
val context: NDAlgebra<T, C, N>
fun unwrap(): N
fun N.wrap(): NDElement<T, C, N>
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 <T : Any, F : Field<T>> buffered(
shape: IntArray,
field: F,
initializer: F.(IntArray) -> T
): BufferedNDElement<T, F> {
val ndField = BoxingNDField(shape, field, Buffer.Companion::boxing)
return ndField.produce(initializer)
}
inline fun <reified T : Any, F : Field<T>> auto(
shape: IntArray,
field: F,
noinline initializer: F.(IntArray) -> T
): BufferedNDFieldElement<T, F> {
val ndField = NDField.auto(field, *shape)
return BufferedNDFieldElement(ndField, ndField.produce(initializer).buffer)
}
}
}
fun <T, C, N : NDStructure<T>> NDElement<T, C, N>.mapIndexed(transform: C.(index: IntArray, T) -> T) =
context.mapIndexed(unwrap(), transform).wrap()
fun <T, C, N : NDStructure<T>> NDElement<T, C, N>.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 <T, C, N : NDStructure<T>> Function1<T, T>.invoke(ndElement: NDElement<T, C, N>) =
ndElement.map { value -> this@invoke(value) }
/* plus and minus */
/**
* Summation operation for [NDElement] and single element
*/
operator fun <T, S : Space<T>, N : NDStructure<T>> NDElement<T, S, N>.plus(arg: T) =
map { value -> arg + value }
/**
* Subtraction operation between [NDElement] and single element
*/
operator fun <T, S : Space<T>, N : NDStructure<T>> NDElement<T, S, N>.minus(arg: T) =
map { value -> arg - value }
/* prod and div */
/**
* Product operation for [NDElement] and single element
*/
operator fun <T, R : Ring<T>, N : NDStructure<T>> NDElement<T, R, N>.times(arg: T) =
map { value -> arg * value }
/**
* Division operation between [NDElement] and single element
*/
operator fun <T, F : Field<T>, N : NDStructure<T>> NDElement<T, F, N>.div(arg: T) =
map { value -> arg / value }
// /**
// * Reverse sum operation
// */
// operator fun T.plus(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
// field.run { this@plus + arg[index] }
// }
//
// /**
// * Reverse minus operation
// */
// operator fun T.minus(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
// field.run { this@minus - arg[index] }
// }
//
// /**
// * Reverse product operation
// */
// operator fun T.times(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
// field.run { this@times * arg[index] }
// }
//
// /**
// * Reverse division operation
// */
// operator fun T.div(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
// field.run { this@div / arg[index] }
// }

View File

@ -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<T, F : Field<T>>(val shape: IntArray, val field: F) : Field<NDElement<T, F>> {
abstract fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T>
/**
* 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<T, F> = NDStructureElement(this, produceStructure(initializer))
override val zero: NDElement<T, F> 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<T, F>) {
elements.forEach {
if (!shape.contentEquals(it.shape)) {
throw ShapeMismatchException(shape, it.shape)
}
}
}
/**
* Element-by-element addition
*/
override fun add(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
checkShape(a, b)
return produce { with(field) { a[it] + b[it] } }
}
/**
* Multiply all elements by cinstant
*/
override fun multiply(a: NDElement<T, F>, k: Double): NDElement<T, F> {
checkShape(a)
return produce { with(field) { a[it] * k } }
}
override val one: NDElement<T, F>
get() = produce { one }
/**
* Element-by-element multiplication
*/
override fun multiply(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
checkShape(a)
return produce { with(field) { a[it] * b[it] } }
}
/**
* Element-by-element division
*/
override fun divide(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
checkShape(a)
return produce { with(field) { a[it] / b[it] } }
}
// /**
// * Reverse sum operation
// */
// operator fun T.plus(arg: NDElement<T, F>): NDElement<T, F> = arg + this
//
// /**
// * Reverse minus operation
// */
// operator fun T.minus(arg: NDElement<T, F>): NDElement<T, F> = arg.transformIndexed { _, value ->
// with(arg.context.field) {
// this@minus - value
// }
// }
//
// /**
// * Reverse product operation
// */
// operator fun T.times(arg: NDElement<T, F>): NDElement<T, F> = arg * this
//
// /**
// * Reverse division operation
// */
// operator fun T.div(arg: NDElement<T, F>): NDElement<T, F> = arg.transformIndexed { _, value ->
// with(arg.context.field) {
// this@div / value
// }
// }
}
interface NDElement<T, F : Field<T>>: FieldElement<NDElement<T, F>, NDField<T, F>>, NDStructure<T>
inline fun <T, F : Field<T>> NDElement<T, F>.transformIndexed(crossinline action: F.(IntArray, T) -> T): NDElement<T, F> = context.produce { action(it, get(*it)) }
inline fun <T, F : Field<T>> NDElement<T, F>.transform(crossinline action: F.(T) -> T): NDElement<T, F> = context.produce { action(get(*it)) }
/**
* Read-only [NDStructure] coupled to the context.
*/
class NDStructureElement<T, F : Field<T>>(override val context: NDField<T, F>, private val structure: NDStructure<T>) : NDElement<T,F>, NDStructure<T> by structure {
//TODO ensure structure is immutable
override val self: NDElement<T, F> get() = this
}
/**
* Element by element application of any operation on elements to the whole array. Just like in numpy
*/
operator fun <T, F : Field<T>> Function1<T, T>.invoke(ndElement: NDElement<T, F>): NDElement<T, F> = ndElement.transform {value -> this@invoke(value) }
/* plus and minus */
/**
* Summation operation for [NDElement] and single element
*/
operator fun <T, F : Field<T>> NDElement<T, F>.plus(arg: T): NDElement<T, F> = transform {value ->
with(context.field) {
arg + value
}
}
/**
* Subtraction operation between [NDElement] and single element
*/
operator fun <T, F : Field<T>> NDElement<T, F>.minus(arg: T): NDElement<T, F> = transform {value ->
with(context.field) {
arg - value
}
}
/* prod and div */
/**
* Product operation for [NDElement] and single element
*/
operator fun <T, F : Field<T>> NDElement<T, F>.times(arg: T): NDElement<T, F> = transform { value ->
with(context.field) {
arg * value
}
}
/**
* Division operation between [NDElement] and single element
*/
operator fun <T, F : Field<T>> NDElement<T, F>.div(arg: T): NDElement<T, F> = transform { value ->
with(context.field) {
arg / value
}
}
class GenericNDField<T : Any, F : Field<T>>(shape: IntArray, field: F) : NDField<T, F>(shape, field) {
override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T> = genericNdStructure(shape) { field.initializer(it) }
}
//typealias NDFieldFactory<T> = (IntArray)->NDField<T>
object NDArrays {
/**
* Create a platform-optimized NDArray of doubles
*/
fun realNDArray(shape: IntArray, initializer: DoubleField.(IntArray) -> Double = { 0.0 }): NDElement<Double, DoubleField> {
return ExtendedNDField(shape, DoubleField).produce(initializer)
}
fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement<Double, DoubleField> {
return realNDArray(intArrayOf(dim)) { initializer(it[0]) }
}
fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement<Double, DoubleField> {
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<Double, DoubleField> {
return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
}
inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double, DoubleField>.() -> NDElement<Double, DoubleField>) =
ExtendedNDField(shape, DoubleField).run(block)
// /**
// * Simple boxing NDField
// */
// fun <T : Any> fieldFactory(field: Field<T>): NDFieldFactory<T> = { shape -> GenericNDField(shape, field) }
/**
* Simple boxing NDArray
*/
fun <T : Any, F : Field<T>> create(field: F, shape: IntArray, initializer: (IntArray) -> T): NDElement<T, F> {
return GenericNDField(shape, field).produce { initializer(it) }
}
}

View File

@ -11,6 +11,45 @@ interface NDStructure<T> {
operator fun get(index: IntArray): T operator fun get(index: IntArray): T
fun elements(): Sequence<Pair<IntArray, T>> fun elements(): Sequence<Pair<IntArray, T>>
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 <T> build(
strides: Strides,
bufferFactory: BufferFactory<T> = 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 <reified T : Any> auto(strides: Strides, crossinline initializer: (IntArray) -> T) =
BufferNDStructure(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
fun <T> build(
shape: IntArray,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
initializer: (IntArray) -> T
) = build(DefaultStrides(shape), bufferFactory, initializer)
inline fun <reified T : Any> auto(shape: IntArray, crossinline initializer: (IntArray) -> T) =
auto(DefaultStrides(shape), initializer)
}
} }
operator fun <T> NDStructure<T>.get(vararg index: Int): T = get(index) operator fun <T> NDStructure<T>.get(vararg index: Int): T = get(index)
@ -19,7 +58,7 @@ interface MutableNDStructure<T> : NDStructure<T> {
operator fun set(index: IntArray, value: T) operator fun set(index: IntArray, value: T)
} }
fun <T> MutableNDStructure<T>.transformInPlace(action: (IntArray, T) -> T) { fun <T> MutableNDStructure<T>.mapInPlace(action: (IntArray, T) -> T) {
elements().forEach { (index, oldValue) -> elements().forEach { (index, oldValue) ->
this[index] = action(index, oldValue) this[index] = action(index, oldValue)
} }
@ -49,6 +88,9 @@ interface Strides {
*/ */
fun index(offset: Int): IntArray fun index(offset: Int): IntArray
/**
* The size of linear buffer to accommodate all elements of ND-structure corresponding to strides
*/
val linearSize: Int 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 * Strides for memory access
*/ */
@ -77,7 +119,7 @@ class DefaultStrides(override val shape: IntArray) : Strides {
override fun offset(index: IntArray): Int { override fun offset(index: IntArray): Int {
return index.mapIndexed { i, value -> 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]})") throw RuntimeException("Index $value out of shape bounds: (0,${this.shape[i]})")
} }
value * strides[i] value * strides[i]
@ -98,19 +140,40 @@ class DefaultStrides(override val shape: IntArray) : Strides {
override val linearSize: Int override val linearSize: Int
get() = strides[shape.size] 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<IntArray, Strides>()
/**
* Cached builder for default strides
*/
operator fun invoke(shape: IntArray): Strides = defaultStridesCache.getOrPut(shape) { DefaultStrides(shape) }
}
} }
abstract class GenericNDStructure<T, B : Buffer<T>> : NDStructure<T> { interface NDBuffer<T> : NDStructure<T> {
protected abstract val buffer: B val buffer: Buffer<T>
protected abstract val strides: Strides val strides: Strides
override fun get(index: IntArray): T = buffer[strides.offset(index)] override fun get(index: IntArray): T = buffer[strides.offset(index)]
override val shape: IntArray override val shape: IntArray get() = strides.shape
get() = strides.shape
override fun elements()= override fun elements() = strides.indices().map { it to this[it] }
strides.indices().map { it to this[it] }
} }
/** /**
@ -119,29 +182,58 @@ abstract class GenericNDStructure<T, B : Buffer<T>> : NDStructure<T> {
class BufferNDStructure<T>( class BufferNDStructure<T>(
override val strides: Strides, override val strides: Strides,
override val buffer: Buffer<T> override val buffer: Buffer<T>
) : GenericNDStructure<T, Buffer<T>>() { ) : NDBuffer<T> {
init { init {
if (strides.linearSize != buffer.size) { if (strides.linearSize != buffer.size) {
error("Expected buffer side of ${strides.linearSize}, but found ${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 <reified T : Any> ndStructure(strides: Strides, noinline initializer: (IntArray) -> T) = /**
BufferNDStructure<T>(strides, buffer(strides.linearSize) { i -> initializer(strides.index(i)) }) * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferNDStructure]
*/
inline fun <reified T : Any> ndStructure(shape: IntArray, noinline initializer: (IntArray) -> T) = inline fun <T, reified R : Any> NDStructure<T>.mapToBuffer(
ndStructure(DefaultStrides(shape), initializer) factory: BufferFactory<R> = Buffer.Companion::auto,
crossinline transform: (T) -> R
): BufferNDStructure<R> {
return if (this is BufferNDStructure<T>) {
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<T>( class MutableBufferNDStructure<T>(
override val strides: Strides, override val strides: Strides,
override val buffer: MutableBuffer<T> override val buffer: MutableBuffer<T>
) : GenericNDStructure<T, MutableBuffer<T>>(), MutableNDStructure<T> { ) : NDBuffer<T>, MutableNDStructure<T> {
init { init {
if (strides.linearSize != buffer.size) { if (strides.linearSize != buffer.size) {
@ -152,25 +244,10 @@ class MutableBufferNDStructure<T>(
override fun set(index: IntArray, value: T) = buffer.set(strides.offset(index), value) override fun set(index: IntArray, value: T) = buffer.set(strides.offset(index), value)
} }
/** inline fun <reified T : Any> NDStructure<T>.combine(
* Create optimized mutable structure for given type struct: NDStructure<T>,
*/ crossinline block: (T, T) -> T
inline fun <reified T : Any> mutableNdStructure(strides: Strides, noinline initializer: (IntArray) -> T) = ): NDStructure<T> {
MutableBufferNDStructure(strides, mutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) }) if (!this.shape.contentEquals(struct.shape)) error("Shape mismatch in structure combination")
return NDStructure.auto(shape) { block(this[it], struct[it]) }
inline fun <reified T : Any> mutableNdStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
mutableNdStructure(DefaultStrides(shape), initializer)
/**
* Create universal mutable structure
*/
fun <T> genericNdStructure(shape: IntArray, initializer: (IntArray) -> T): MutableBufferNDStructure<T> {
val strides = DefaultStrides(shape)
val sequence = sequence {
strides.indices().forEach {
yield(initializer(it))
}
}
val buffer = MutableListBuffer(sequence.toMutableList())
return MutableBufferNDStructure(strides, buffer)
} }

View File

@ -0,0 +1,48 @@
package scientifik.kmath.structures
import scientifik.memory.*
/**
* A non-boxing buffer based on [ByteBuffer] storage
*/
open class ObjectBuffer<T : Any>(protected val memory: Memory, protected val spec: MemorySpec<T>) : Buffer<T> {
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<T> = (0 until size).asSequence().map { get(it) }.iterator()
companion object {
fun <T : Any> create(spec: MemorySpec<T>, size: Int) =
ObjectBuffer(Memory.allocate(size * spec.objectSize), spec)
inline fun <T : Any> create(
spec: MemorySpec<T>,
size: Int,
crossinline initializer: (Int) -> T
): ObjectBuffer<T> =
MutableObjectBuffer(Memory.allocate(size * spec.objectSize), spec).also { buffer ->
(0 until size).forEach {
buffer[it] = initializer(it)
}
}
}
}
class MutableObjectBuffer<T : Any>(memory: Memory, spec: MemorySpec<T>) : ObjectBuffer<T>(memory, spec),
MutableBuffer<T> {
private val writer = memory.writer()
override fun set(index: Int, value: T) = writer.write(spec, spec.objectSize * index, value)
override fun copy(): MutableBuffer<T> = MutableObjectBuffer(memory.copy(), spec)
companion object {
}
}

View File

@ -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<Buffer<Double>> {
override fun add(a: Buffer<Double>, b: Buffer<Double>): 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<Double>, 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<Double>, b: Buffer<Double>): 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<Double>, b: Buffer<Double>): 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<Double>): 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<Double>): 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<Double>, 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<Double>): 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<Double>): 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<Buffer<Double>>, ExtendedFieldOperations<Buffer<Double>> {
override val zero: Buffer<Double> by lazy { DoubleBuffer(size) { 0.0 } }
override val one: Buffer<Double> by lazy { DoubleBuffer(size) { 1.0 } }
override fun add(a: Buffer<Double>, b: Buffer<Double>): 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<Double>, 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<Double>, b: Buffer<Double>): 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<Double>, b: Buffer<Double>): 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<Double>): 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<Double>): 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<Double>, 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<Double>): 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<Double>): DoubleBuffer {
require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
return RealBufferFieldOperations.ln(arg)
}
}

View File

@ -0,0 +1,121 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.FieldElement
import scientifik.kmath.operations.RealField
typealias RealNDElement = BufferedNDFieldElement<Double, RealField>
class RealNDField(override val shape: IntArray) :
BufferedNDField<Double, RealField>,
ExtendedNDField<Double, RealField, NDBuffer<Double>> {
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<Double> =
DoubleBuffer(DoubleArray(size) { initializer(it) })
/**
* Inline transform an NDStructure to
*/
override fun map(
arg: NDBuffer<Double>,
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<Double>,
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<Double>,
b: NDBuffer<Double>,
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<Double>.toElement(): FieldElement<NDBuffer<Double>, *, out BufferedNDField<Double, RealField>> =
BufferedNDFieldElement(this@RealNDField, buffer)
override fun power(arg: NDBuffer<Double>, pow: Number) = map(arg) { power(it, pow) }
override fun exp(arg: NDBuffer<Double>) = map(arg) { exp(it) }
override fun ln(arg: NDBuffer<Double>) = map(arg) { ln(it) }
override fun sin(arg: NDBuffer<Double>) = map(arg) { sin(it) }
override fun cos(arg: NDBuffer<Double>) = map(arg) { cos(it) }
}
/**
* Fast element production using function inlining
*/
inline fun BufferedNDField<Double, RealField>.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<Double, Double>.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 }

View File

@ -0,0 +1,96 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.RingElement
import scientifik.kmath.operations.ShortRing
typealias ShortNDElement = BufferedNDRingElement<Short, ShortRing>
class ShortNDRing(override val shape: IntArray) :
BufferedNDRing<Short, ShortRing> {
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<Short> =
ShortBuffer(ShortArray(size) { initializer(it) })
/**
* Inline transform an NDStructure to
*/
override fun map(
arg: NDBuffer<Short>,
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<Short>,
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<Short>,
b: NDBuffer<Short>,
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<Short>.toElement(): RingElement<NDBuffer<Short>, *, out BufferedNDRing<Short, ShortRing>> =
BufferedNDRingElement(this@ShortNDRing, buffer)
}
/**
* Fast element production using function inlining
*/
inline fun BufferedNDRing<Short, ShortRing>.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<Short, Short>.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() }

View File

@ -0,0 +1,94 @@
package scientifik.kmath.structures
/**
* A structure that is guaranteed to be one-dimensional
*/
interface Structure1D<T> : NDStructure<T>, Buffer<T> {
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<T> = (0 until size).asSequence().map { get(it) }.iterator()
}
/**
* A 1D wrapper for nd-structure
*/
private inline class Structure1DWrapper<T>(val structure: NDStructure<T>) : Structure1D<T> {
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<Pair<IntArray, T>> = structure.elements()
}
/**
* Represent a [NDStructure] as [Structure1D]. Throw error in case of dimension mismatch
*/
fun <T> NDStructure<T>.as1D(): Structure1D<T> = if (shape.size == 1) {
Structure1DWrapper(this)
} else {
error("Can't create 1d-structure from ${shape.size}d-structure")
}
fun <T> NDBuffer<T>.as1D(): Structure1D<T> = 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<T>(val buffer: Buffer<T>) : Structure1D<T> {
override val shape: IntArray get() = intArrayOf(buffer.size)
override val size: Int get() = buffer.size
override fun elements(): Sequence<Pair<IntArray, T>> =
asSequence().mapIndexed { index, value -> intArrayOf(index) to value }
override fun get(index: Int): T = buffer.get(index)
}
/**
* Represent this buffer as 1D structure
*/
fun <T> Buffer<T>.asND(): Structure1D<T> = Buffer1DWrapper(this)
/**
* A structure that is guaranteed to be two-dimensional
*/
interface Structure2D<T> : NDStructure<T> {
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<T>(val structure: NDStructure<T>) : Structure2D<T> {
override fun get(i: Int, j: Int): T = structure[i, j]
override val shape: IntArray get() = structure.shape
override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()
}
/**
* Represent a [NDStructure] as [Structure1D]. Throw error in case of dimension mismatch
*/
fun <T> NDStructure<T>.as2D(): Structure2D<T> = if (shape.size == 2) {
Structure2DWrapper(this)
} else {
error("Can't create 2d-structure from ${shape.size}d-structure")
}

View File

@ -2,17 +2,17 @@ package scientifik.kmath.expressions
import scientifik.kmath.operations.Complex import scientifik.kmath.operations.Complex
import scientifik.kmath.operations.ComplexField import scientifik.kmath.operations.ComplexField
import scientifik.kmath.operations.DoubleField import scientifik.kmath.operations.RealField
import kotlin.test.Test import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
class FieldExpressionContextTest { class ExpressionFieldTest {
@Test @Test
fun testExpression() { fun testExpression() {
val context = FieldExpressionContext(DoubleField) val context = ExpressionField(RealField)
val expression = with(context) { val expression = with(context) {
val x = variable("x", 2.0) 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("x" to 1.0), 4.0)
assertEquals(expression(), 9.0) assertEquals(expression(), 9.0)
@ -20,10 +20,10 @@ class FieldExpressionContextTest {
@Test @Test
fun testComplex() { fun testComplex() {
val context = FieldExpressionContext(ComplexField) val context = ExpressionField(ComplexField)
val expression = with(context) { val expression = with(context) {
val x = variable("x", Complex(2.0, 0.0)) 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("x" to Complex(1.0, 0.0)), Complex(4.0, 0.0))
assertEquals(expression(), Complex(9.0, 0.0)) assertEquals(expression(), Complex(9.0, 0.0))
@ -31,23 +31,23 @@ class FieldExpressionContextTest {
@Test @Test
fun separateContext() { fun separateContext() {
fun <T> FieldExpressionContext<T>.expression(): Expression<T>{ fun <T> ExpressionField<T>.expression(): Expression<T> {
val x = variable("x") 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) assertEquals(expression("x" to 1.0), 4.0)
} }
@Test @Test
fun valueExpression() { fun valueExpression() {
val expressionBuilder: FieldExpressionContext<Double>.()->Expression<Double> = { val expressionBuilder: ExpressionField<Double>.() -> Expression<Double> = {
val x = variable("x") 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) assertEquals(expression("x" to 1.0), 4.0)
} }
} }

View File

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

View File

@ -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<Double>(2, 3)(
1.0, 0.0, 0.0,
0.0, 1.0, 2.0
)
assertEquals(2.0, matrix[1, 2])
}
}

View File

@ -1,21 +1,41 @@
package scientifik.kmath.linear package scientifik.kmath.linear
import scientifik.kmath.operations.DoubleField
import kotlin.test.Test import kotlin.test.Test
import kotlin.test.assertTrue import kotlin.test.assertEquals
class RealLUSolverTest { class RealLUSolverTest {
@Test @Test
fun testInvertOne() { fun testInvertOne() {
val matrix = Matrix.diagonal(2, 2, DoubleField) val matrix = MatrixContext.real.one(2, 2)
val inverted = RealLUSolver.inverse(matrix) val inverted = LUSolver.real.inverse(matrix)
assertTrue { Matrix.equals(matrix,inverted) } assertEquals(matrix, inverted)
} }
// @Test @Test
// fun testInvert() { fun testInvert() {
// val matrix = realMatrix(2,2){} val matrix = Matrix.square(
// val inverted = RealLUSolver.inverse(matrix) 3.0, 1.0,
// assertTrue { Matrix.equals(matrix,inverted) } 1.0, 3.0
// } )
val decomposed = LUSolver.real.decompose(matrix)
val decomposition = decomposed.getFeature<LUPDecomposition<Double>>()!!
//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)
}
} }

View File

@ -6,11 +6,9 @@ import kotlin.test.assertEquals
class RealFieldTest { class RealFieldTest {
@Test @Test
fun testSqrt() { fun testSqrt() {
//fails because KT-27586
val sqrt = with(RealField) { val sqrt = with(RealField) {
sqrt( 25 * one) sqrt(25 * one)
} }
assertEquals(5.0, sqrt.value) assertEquals(5.0, sqrt)
} }
} }

View File

@ -1,17 +1,14 @@
package scientifik.kmath.structures package scientifik.kmath.structures
import org.junit.Test
import scientifik.kmath.operations.Complex import scientifik.kmath.operations.Complex
import scientifik.kmath.operations.complex
import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
class ComplexBufferSpecTest { class ComplexBufferSpecTest {
@Test @Test
fun testComplexBuffer() { fun testComplexBuffer() {
val buffer = Complex.createBuffer(20) val buffer = Buffer.complex(20) { Complex(it.toDouble(), -it.toDouble()) }
(0 until 20).forEach {
buffer[it] = Complex(it.toDouble(), -it.toDouble())
}
assertEquals(Complex(5.0, -5.0), buffer[5]) assertEquals(Complex(5.0, -5.0), buffer[5])
} }
} }

View File

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

View File

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

View File

@ -1,16 +1,15 @@
package scientifik.kmath.structures package scientifik.kmath.structures
import scientifik.kmath.operations.Norm import scientifik.kmath.operations.Norm
import scientifik.kmath.structures.NDArrays.produceReal import scientifik.kmath.structures.NDElement.Companion.real2D
import scientifik.kmath.structures.NDArrays.real2DArray
import kotlin.math.abs import kotlin.math.abs
import kotlin.math.pow import kotlin.math.pow
import kotlin.test.Test import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
class NumberNDFieldTest { class NumberNDFieldTest {
val array1 = real2DArray(3, 3) { i, j -> (i + j).toDouble() } val array1 = real2D(3, 3) { i, j -> (i + j).toDouble() }
val array2 = real2DArray(3, 3) { i, j -> (i - j).toDouble() } val array2 = real2D(3, 3) { i, j -> (i - j).toDouble() }
@Test @Test
fun testSum() { fun testSum() {
@ -27,7 +26,7 @@ class NumberNDFieldTest {
@Test @Test
fun testGeneration() { 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 (i in 0..2) {
for (j in 0..2) { for (j in 0..2) {
@ -51,15 +50,20 @@ class NumberNDFieldTest {
assertEquals(2.0, result[0, 2]) assertEquals(2.0, result[0, 2])
} }
object L2Norm : Norm<NDElement<out Number, *>, Double> { @Test
override fun norm(arg: NDElement<out Number, *>): Double { fun combineTest() {
return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() }) val division = array1.combine(array2, Double::div)
}
object L2Norm : Norm<NDStructure<out Number>, Double> {
override fun norm(arg: NDStructure<out Number>): Double {
return kotlin.math.sqrt(arg.elements().sumByDouble { it.second.toDouble() })
} }
} }
@Test @Test
fun testInternalContext() { fun testInternalContext() {
produceReal(array1.shape) { NDField.real(*array1.shape).run {
with(L2Norm) { with(L2Norm) {
1 + norm(array1) + exp(array2) 1 + norm(array1) + exp(array2)
} }

View File

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

View File

@ -1,55 +0,0 @@
package scientifik.kmath.structures
import java.nio.ByteBuffer
/**
* A specification for serialization and deserialization objects to buffer
*/
interface BufferSpec<T : Any> {
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<T : Any> : BufferSpec<T> {
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)
}
}

View File

@ -1,25 +0,0 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.Complex
import java.nio.ByteBuffer
object ComplexBufferSpec : FixedSizeBufferSpec<Complex> {
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)

View File

@ -1,28 +0,0 @@
package scientifik.kmath.structures
import java.nio.ByteBuffer
class ObjectBuffer<T : Any>(private val buffer: ByteBuffer, private val spec: FixedSizeBufferSpec<T>) : MutableBuffer<T> {
override val size: Int
get() = buffer.limit() / spec.unitSize
override fun get(index: Int): T = with(spec) { buffer.readObject(index) }
override fun iterator(): Iterator<T> = (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<T> {
val dup = buffer.duplicate()
val copy = ByteBuffer.allocate(dup.capacity())
dup.rewind()
copy.put(dup)
copy.flip()
return ObjectBuffer(copy, spec)
}
companion object {
fun <T : Any> create(spec: FixedSizeBufferSpec<T>, size: Int) =
ObjectBuffer<T>(ByteBuffer.allocate(size * spec.unitSize), spec)
}
}

View File

@ -1,23 +0,0 @@
package scientifik.kmath.structures
import scientifik.kmath.operations.Real
import java.nio.ByteBuffer
object RealBufferSpec : FixedSizeBufferSpec<Real> {
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<Double> {
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)

View File

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

View File

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

View File

@ -6,6 +6,4 @@ import kotlinx.coroutines.Dispatchers
import kotlin.coroutines.CoroutineContext import kotlin.coroutines.CoroutineContext
import kotlin.coroutines.EmptyCoroutineContext import kotlin.coroutines.EmptyCoroutineContext
expect fun <R> runBlocking(context: CoroutineContext = EmptyCoroutineContext, function: suspend CoroutineScope.()->R): R
val Dispatchers.Math: CoroutineDispatcher get() = Dispatchers.Default val Dispatchers.Math: CoroutineDispatcher get() = Dispatchers.Default

View File

@ -1,79 +0,0 @@
package scientifik.kmath.structures
import kotlinx.coroutines.*
import scientifik.kmath.operations.Field
class LazyNDField<T, F : Field<T>>(shape: IntArray, field: F, val scope: CoroutineScope = GlobalScope) : NDField<T, F>(shape, field) {
override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T> = LazyNDStructure(this) { initializer(field, it) }
override fun add(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
return LazyNDStructure(this) { index ->
val aDeferred = a.deferred(index)
val bDeferred = b.deferred(index)
aDeferred.await() + bDeferred.await()
}
}
override fun multiply(a: NDElement<T, F>, k: Double): NDElement<T, F> {
return LazyNDStructure(this) { index -> a.await(index) * k }
}
override fun multiply(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
return LazyNDStructure(this) { index ->
val aDeferred = a.deferred(index)
val bDeferred = b.deferred(index)
aDeferred.await() * bDeferred.await()
}
}
override fun divide(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
return LazyNDStructure(this) { index ->
val aDeferred = a.deferred(index)
val bDeferred = b.deferred(index)
aDeferred.await() / bDeferred.await()
}
}
}
class LazyNDStructure<T, F : Field<T>>(override val context: LazyNDField<T, F>, val function: suspend F.(IntArray) -> T) : NDElement<T, F>, NDStructure<T> {
override val self: NDElement<T, F> get() = this
override val shape: IntArray get() = context.shape
private val cache = HashMap<IntArray, Deferred<T>>()
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<Pair<IntArray, T>> {
val strides = DefaultStrides(shape)
return strides.indices().map { index -> index to runBlocking { await(index) } }
}
}
fun <T> NDElement<T, *>.deferred(index: IntArray) = if (this is LazyNDStructure<T, *>) this.deferred(index) else CompletableDeferred(get(index))
suspend fun <T> NDElement<T, *>.await(index: IntArray) = if (this is LazyNDStructure<T, *>) this.await(index) else get(index)
fun <T, F : Field<T>> NDElement<T, F>.lazy(scope: CoroutineScope = GlobalScope): LazyNDStructure<T, F> {
return if (this is LazyNDStructure<T, F>) {
this
} else {
val context = LazyNDField(context.shape, context.field)
LazyNDStructure(context) { get(it) }
}
}
inline fun <T, F : Field<T>> LazyNDStructure<T, F>.transformIndexed(crossinline action: suspend F.(IntArray, T) -> T) = LazyNDStructure(context) { index ->
action.invoke(this, index, await(index))
}
inline fun <T, F : Field<T>> LazyNDStructure<T, F>.transform(crossinline action: suspend F.(T) -> T) = LazyNDStructure(context) { index ->
action.invoke(this, await(index))
}

View File

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

View File

@ -0,0 +1,47 @@
package scientifik.kmath.structures
import kotlinx.coroutines.*
class LazyNDStructure<T>(
val scope: CoroutineScope,
override val shape: IntArray,
val function: suspend (IntArray) -> T
) : NDStructure<T> {
private val cache = HashMap<IntArray, Deferred<T>>()
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<Pair<IntArray, T>> {
val strides = DefaultStrides(shape)
val res = runBlocking {
strides.indices().toList().map { index -> index to await(index) }
}
return res.asSequence()
}
}
fun <T> NDStructure<T>.deferred(index: IntArray) =
if (this is LazyNDStructure<T>) this.deferred(index) else CompletableDeferred(get(index))
suspend fun <T> NDStructure<T>.await(index: IntArray) =
if (this is LazyNDStructure<T>) this.await(index) else get(index)
/**
* PENDING would benifit from KEEP-176
*/
fun <T, R> NDStructure<T>.mapAsyncIndexed(scope: CoroutineScope, function: suspend (T, index: IntArray) -> R) =
LazyNDStructure(scope, shape) { index -> function(get(index), index) }
fun <T, R> NDStructure<T>.mapAsync(scope: CoroutineScope, function: suspend (T) -> R) =
LazyNDStructure(scope, shape) { index -> function(get(index)) }

View File

@ -1,6 +0,0 @@
package scientifik.kmath.structures
import kotlinx.coroutines.CoroutineScope
import kotlin.coroutines.CoroutineContext
actual fun <R> runBlocking(context: CoroutineContext, function: suspend CoroutineScope.() -> R): R = kotlinx.coroutines.runBlocking(context, function)

View File

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

View File

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

View File

@ -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<T : Any> {
operator fun contains(vector: Point<out T>): Boolean
val dimension: Int
}
/**
* The bin in the histogram. The histogram is by definition always done in the real space
*/
interface Bin<T : Any> : Domain<T> {
/**
* The value of this bin
*/
val value: Number
val center: Point<T>
}
interface Histogram<T : Any, out B : Bin<T>> : Iterable<B> {
/**
* Find existing bin, corresponding to given coordinates
*/
operator fun get(point: Point<out T>): B?
/**
* Dimension of the histogram
*/
val dimension: Int
}
interface MutableHistogram<T : Any, out B : Bin<T>> : Histogram<T, B> {
/**
* Increment appropriate bin
*/
fun putWithWeight(point: Point<out T>, weight: Double)
fun put(point: Point<out T>) = putWithWeight(point, 1.0)
}
fun <T : Any> MutableHistogram<T, *>.put(vararg point: T) = put(ArrayBuffer(point))
fun MutableHistogram<Double, *>.put(vararg point: Number) =
put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray()))
fun MutableHistogram<Double, *>.put(vararg point: Double) = put(DoubleBuffer(point))
fun <T : Any> MutableHistogram<T, *>.fill(sequence: Iterable<Point<T>>) = sequence.forEach { put(it) }
/**
* Pass a sequence builder into histogram
*/
fun <T : Any> MutableHistogram<T, *>.fill(buider: suspend SequenceScope<Point<T>>.() -> Unit) =
fill(sequence(buider).asIterable())

View File

@ -1,43 +1,66 @@
package scientifik.kmath.histogram package scientifik.kmath.histogram
import scientifik.kmath.linear.Point
import scientifik.kmath.linear.toVector import scientifik.kmath.linear.toVector
import scientifik.kmath.operations.SpaceOperations
import scientifik.kmath.structures.* import scientifik.kmath.structures.*
import kotlin.math.floor import kotlin.math.floor
private operator fun RealPoint.minus(other: RealPoint) = ListBuffer((0 until size).map { get(it) - other[it] })
private inline fun <T> Buffer<out Double>.mapIndexed(crossinline mapper: (Int, Double) -> T): Sequence<T> = (0 until size).asSequence().map { mapper(it, get(it)) } data class BinDef<T : Comparable<T>>(val space: SpaceOperations<Point<T>>, val center: Point<T>, val sizes: Point<T>) {
fun contains(vector: Point<out T>): 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<T : Comparable<T>>(val def: BinDef<T>, override val value: Number) : Bin<T> {
override fun contains(vector: Point<out T>): Boolean = def.contains(vector)
override val dimension: Int
get() = def.center.size
override val center: Point<T>
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. * 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( class RealHistogram(
private val lower: RealPoint, private val lower: Buffer<Double>,
private val upper: RealPoint, private val upper: Buffer<Double>,
private val binNums: IntArray = IntArray(lower.size) { 20 } private val binNums: IntArray = IntArray(lower.size) { 20 }
) : MutableHistogram<Double, PhantomBin<Double>> { ) : MutableHistogram<Double, MultivariateBin<Double>> {
private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 }) private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 })
private val values: NDStructure<LongCounter> = ndStructure(strides) { LongCounter() } private val values: NDStructure<LongCounter> = NDStructure.auto(strides) { LongCounter() }
//private val weight: NDStructure<DoubleCounter?> = ndStructure(strides){null} //private val weight: NDStructure<DoubleCounter?> = 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 { init {
// argument checks // argument checks
if (lower.size != upper.size) error("Dimension mismatch in histogram lower and upper limits.") 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 (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 * Get internal [NDStructure] bin index for given axis
*/ */
@ -59,49 +82,41 @@ class FastHistogram(
return getValue(getIndex(point)) return getValue(getIndex(point))
} }
private fun getTemplate(index: IntArray): BinTemplate<Double> { private fun getDef(index: IntArray): BinDef<Double> {
val center = index.mapIndexed { axis, i -> val center = index.mapIndexed { axis, i ->
when (i) { when (i) {
0 -> Double.NEGATIVE_INFINITY 0 -> Double.NEGATIVE_INFINITY
strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
else -> lower[axis] + (i.toDouble() - 0.5) * binSize[axis] else -> lower[axis] + (i.toDouble() - 0.5) * binSize[axis]
} }
}.toVector() }.asBuffer()
return BinTemplate(center, binSize) return BinDef(RealBufferFieldOperations, center, binSize)
} }
fun getTemplate(point: Buffer<out Double>): BinTemplate<Double> { fun getDef(point: Buffer<out Double>): BinDef<Double> {
return getTemplate(getIndex(point)) return getDef(getIndex(point))
} }
override fun get(point: Buffer<out Double>): PhantomBin<Double>? { override fun get(point: Buffer<out Double>): MultivariateBin<Double>? {
val index = getIndex(point) val index = getIndex(point)
return PhantomBin(getTemplate(index), getValue(index)) return MultivariateBin(getDef(index), getValue(index))
} }
override fun put(point: Buffer<out Double>, weight: Double) { override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
if (weight != 1.0) TODO("Implement weighting") if (weight != 1.0) TODO("Implement weighting")
val index = getIndex(point) val index = getIndex(point)
values[index].increment() values[index].increment()
} }
override fun iterator(): Iterator<PhantomBin<Double>> = values.elements().map { (index, value) -> override fun iterator(): Iterator<MultivariateBin<Double>> = values.elements().map { (index, value) ->
PhantomBin(getTemplate(index), value.sum()) MultivariateBin(getDef(index), value.sum())
}.iterator() }.iterator()
/** /**
* Convert this histogram into NDStructure containing bin values but not bin descriptions * Convert this histogram into NDStructure containing bin values but not bin descriptions
*/ */
fun asNDStructure(): NDStructure<Number> { fun asNDStructure(): NDStructure<Number> {
return ndStructure(this.values.shape) { values[it].sum() } return NDStructure.auto(this.values.shape) { values[it].sum() }
}
/**
* Create a phantom lightweight immutable copy of this histogram
*/
fun asPhantomHistogram(): PhantomHistogram<Double> {
val binTemplates = values.elements().associate { (index, _) -> getTemplate(index) to index }
return PhantomHistogram(binTemplates, asNDStructure())
} }
companion object { companion object {
@ -115,8 +130,11 @@ class FastHistogram(
*) *)
*``` *```
*/ */
fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): FastHistogram { fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram {
return FastHistogram(ranges.map { it.start }.toVector(), ranges.map { it.endInclusive }.toVector()) return RealHistogram(
ranges.map { it.start }.toVector(),
ranges.map { it.endInclusive }.toVector()
)
} }
/** /**
@ -128,8 +146,8 @@ class FastHistogram(
*) *)
*``` *```
*/ */
fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): FastHistogram { fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): RealHistogram {
return FastHistogram( return RealHistogram(
ListBuffer(ranges.map { it.first.start }), ListBuffer(ranges.map { it.first.start }),
ListBuffer(ranges.map { it.first.endInclusive }), ListBuffer(ranges.map { it.first.endInclusive }),
ranges.map { it.second }.toIntArray() ranges.map { it.second }.toIntArray()

View File

@ -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 scientifik.kmath.linear.Vector
import kotlin.random.Random import kotlin.random.Random
import kotlin.test.Test import kotlin.test.Test
@ -10,7 +13,7 @@ import kotlin.test.assertTrue
class MultivariateHistogramTest { class MultivariateHistogramTest {
@Test @Test
fun testSinglePutHistogram() { fun testSinglePutHistogram() {
val histogram = FastHistogram.fromRanges( val histogram = RealHistogram.fromRanges(
(-1.0..1.0), (-1.0..1.0),
(-1.0..1.0) (-1.0..1.0)
) )
@ -22,21 +25,21 @@ class MultivariateHistogramTest {
} }
@Test @Test
fun testSequentialPut(){ fun testSequentialPut() {
val histogram = FastHistogram.fromRanges( val histogram = RealHistogram.fromRanges(
(-1.0..1.0), (-1.0..1.0),
(-1.0..1.0), (-1.0..1.0),
(-1.0..1.0) (-1.0..1.0)
) )
val random = Random(1234) val random = Random(1234)
fun nextDouble() = random.nextDouble(-1.0,1.0) fun nextDouble() = random.nextDouble(-1.0, 1.0)
val n = 10000 val n = 10000
histogram.fill { histogram.fill {
repeat(n){ repeat(n) {
yield(Vector.ofReal(nextDouble(),nextDouble(),nextDouble())) yield(Vector.ofReal(nextDouble(), nextDouble(), nextDouble()))
} }
} }
assertEquals(n, histogram.sumBy { it.value.toInt() }) assertEquals(n, histogram.sumBy { it.value.toInt() })

View File

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

View File

@ -18,7 +18,7 @@ class UnivariateBin(val position: Double, val size: Double, val counter: LongCou
override fun contains(vector: Buffer<out Double>): Boolean = contains(vector[0]) override fun contains(vector: Buffer<out Double>): 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 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 * Univariate histogram with log(n) bin search speed
*/ */
class UnivariateHistogram private constructor(private val factory: (Double) -> UnivariateBin) : MutableHistogram<Double,UnivariateBin> { class UnivariateHistogram private constructor(private val factory: (Double) -> UnivariateBin) :
MutableHistogram<Double, UnivariateBin> {
private val bins: TreeMap<Double, UnivariateBin> = TreeMap() private val bins: TreeMap<Double, UnivariateBin> = TreeMap()
@ -58,7 +59,7 @@ class UnivariateHistogram private constructor(private val factory: (Double) -> U
(get(value) ?: createBin(value)).inc() (get(value) ?: createBin(value)).inc()
} }
override fun put(point: Buffer<out Double>, weight: Double) { override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
if (weight != 1.0) TODO("Implement weighting") if (weight != 1.0) TODO("Implement weighting")
put(point[0]) put(point[0])
} }
@ -75,8 +76,14 @@ class UnivariateHistogram private constructor(private val factory: (Double) -> U
val sorted = borders.sortedArray() val sorted = borders.sortedArray()
return UnivariateHistogram { value -> return UnivariateHistogram { value ->
when { when {
value < sorted.first() -> UnivariateBin(Double.NEGATIVE_INFINITY, Double.MAX_VALUE) value < sorted.first() -> UnivariateBin(
value > sorted.last() -> UnivariateBin(Double.POSITIVE_INFINITY, Double.MAX_VALUE) Double.NEGATIVE_INFINITY,
Double.MAX_VALUE
)
value > sorted.last() -> UnivariateBin(
Double.POSITIVE_INFINITY,
Double.MAX_VALUE
)
else -> { else -> {
val index = (0 until sorted.size).first { value > sorted[it] } val index = (0 until sorted.size).first { value > sorted[it] }
val left = sorted[index] val left = sorted[index]

View File

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

Some files were not shown because too many files have changed in this diff Show More