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

Merge branch 'dev' into publish-to-space-repo

Browse Source 
# Conflicts:
#	.space.kts
#	build.gradle.kts
#	examples/build.gradle.kts
This commit is contained in:
Iaroslav Postovalov 2020-12-22 20:28:50 +07:00
commit 66671964fa
No known key found for this signature in database
GPG Key ID: 46E15E4A31B3BCD7
122 changed files with 3894 additions and 1445 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
.github/workflows/gradle.yml vendored
View File

@ -1,17 +1,101 @@
name: Gradle build
on: [push]
on: [ push ]
jobs:
build:
runs-on: ubuntu-latest
build-ubuntu:
runs-on: ubuntu-20.04
steps:
- uses: actions/checkout@v1
- name: Set up JDK 11
uses: actions/setup-java@v1
with:
java-version: 11
- name: Build with Gradle
run: ./gradlew build
- uses: actions/checkout@v2
- name: Set up JDK 11
uses: actions/setup-java@v1
with:
java-version: 11
- name: Install Chrome
run: |
sudo apt install -y libappindicator1 fonts-liberation
wget https://dl.google.com/linux/direct/google-chrome-stable_current_amd64.deb
sudo dpkg -i google-chrome*.deb
- name: Cache gradle
uses: actions/cache@v2
with:
path: |
.gradle
build
~/.gradle
key: gradle
restore-keys: gradle
- name: Cache konan
uses: actions/cache@v2
with:
path: |
~/.konan/dependencies
~/.konan/kotlin-native-prebuilt-linux-*
key: ${{ runner.os }}-konan
restore-keys: ${{ runner.os }}-konan
- name: Build with Gradle
run: ./gradlew -Dorg.gradle.daemon=false --build-cache build
build-osx:
runs-on: macos-latest
steps:
- uses: actions/checkout@v2
- name: Set up JDK 11
uses: actions/setup-java@v1
with:
java-version: 11
- name: Cache gradle
uses: actions/cache@v2
with:
path: |
.gradle
build
~/.gradle
key: gradle
restore-keys: gradle
- name: Cache konan
uses: actions/cache@v2
with:
path: |
~/.konan/dependencies
~/.konan/kotlin-native-prebuilt-macos-*
key: ${{ runner.os }}-konan
restore-keys: ${{ runner.os }}-konan
- name: Build with Gradle
run: sudo ./gradlew -Dorg.gradle.daemon=false --build-cache build
build-windows:
runs-on: windows-latest
steps:
- uses: actions/checkout@v2
- name: Set up JDK 11
uses: actions/setup-java@v1
with:
java-version: 11
- name: Add msys to path
run: SETX PATH "%PATH%;C:\msys64\mingw64\bin"
- name: Cache gradle
uses: actions/cache@v2
with:
path: |
.gradle
build
~/.gradle
key: ${{ runner.os }}-gradle
restore-keys: ${{ runner.os }}-gradle
- name: Cache konan
uses: actions/cache@v2
with:
path: |
~/.konan/dependencies
~/.konan/kotlin-native-prebuilt-mingw-*
key: ${{ runner.os }}-konan
restore-keys: ${{ runner.os }}-konan
- name: Build with Gradle
run: ./gradlew --build-cache build

24
CHANGELOG.md
View File

@ -7,21 +7,39 @@
- Better trigonometric and hyperbolic functions for `AutoDiffField` (https://github.com/mipt-npm/kmath/pull/140).
- Automatic README generation for features (#139)
- Native support for `memory`, `core` and `dimensions`
- `kmath-ejml` to supply EJML SimpleMatrix wrapper.
- `kmath-ejml` to supply EJML SimpleMatrix wrapper (https://github.com/mipt-npm/kmath/pull/136).
- A separate `Symbol` entity, which is used for global unbound symbol.
- A `Symbol` indexing scope.
- Basic optimization API for Commons-math.
- Chi squared optimization for array-like data in CM
- `Fitting` utility object in prob/stat
- ND4J support module submitting `NDStructure` and `NDAlgebra` over `INDArray`.
- Coroutine-deterministic Monte-Carlo scope with a random number generator.
- Some minor utilities to `kmath-for-real`.
- Generic operation result parameter to `MatrixContext`
### Changed
- Package changed from `scientifik` to `kscience.kmath`.
- Gradle version: 6.6 -> 6.6.1
- Gradle version: 6.6 -> 6.7.1
- Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`)
- `Polynomial` secondary constructor made function.
- Kotlin version: 1.3.72 -> 1.4.20-M1
- Kotlin version: 1.3.72 -> 1.4.20
- `kmath-ast` doesn't depend on heavy `kotlin-reflect` library.
- Full autodiff refactoring based on `Symbol`
- `kmath-prob` renamed to `kmath-stat`
- Grid generators moved to `kmath-for-real`
- Use `Point<Double>` instead of specialized type in `kmath-for-real`
- Optimized dot product for buffer matrices moved to `kmath-for-real`
- EjmlMatrix context is an object
- Matrix LUP `inverse` renamed to `inverseWithLUP`
### Deprecated
### Removed
- `kmath-koma` module because it doesn't support Kotlin 1.4.
- Support of `legacy` JS backend (we will support only IR)
- `toGrid` method.
- Public visibility of `BufferAccessor2D`
### Fixed
- `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140)

120
README.md
View File

@ -8,41 +8,50 @@ Bintray: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience
Bintray-dev: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion)
# KMath
Could be pronounced as `key-math`.
The Kotlin MATHematics library was initially intended as a Kotlin-based analog to Python's `numpy` library. Later we found that kotlin is much more flexible language and allows superior architecture designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could be achieved with [kmath-for-real](/kmath-for-real) extension module.
Could be pronounced as `key-math`. The Kotlin MATHematics library was initially intended as a Kotlin-based analog to
Python's NumPy library. Later we found that kotlin is much more flexible language and allows superior architecture
designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could
be achieved with [kmath-for-real](/kmath-for-real) extension module.
## Publications and talks
* [A conceptual article about context-oriented design](https://proandroiddev.com/an-introduction-context-oriented-programming-in-kotlin-2e79d316b0a2)
* [Another article about context-oriented design](https://proandroiddev.com/diving-deeper-into-context-oriented-programming-in-kotlin-3ecb4ec38814)
* [ACAT 2019 conference paper](https://aip.scitation.org/doi/abs/10.1063/1.5130103)
# Goal
* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM and JS for now and Native in future).
* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS and Native).
* Provide basic multiplatform implementations for those abstractions (without significant performance optimization).
* Provide bindings and wrappers with those abstractions for popular optimized platform libraries.
## Non-goals
* Be like Numpy. It was the idea at the beginning, but we decided that we can do better in terms of API.
* Provide best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
* Be like NumPy. It was the idea at the beginning, but we decided that we can do better in terms of API.
* Provide the best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
* Cover all cases as immediately and in one bundle. We will modularize everything and add new features gradually.
* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better experience for those, who want to work with specific types.
* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like
for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better
experience for those, who want to work with specific types.
## Features
Actual feature list is [here](/docs/features.md)
Current feature list is [here](/docs/features.md)
* **Algebra**
* Algebraic structures like rings, spaces and field (**TODO** add example to wiki)
* Algebraic structures like rings, spaces and fields (**TODO** add example to wiki)
* Basic linear algebra operations (sums, products, etc.), backed by the `Space` API.
* Complex numbers backed by the `Field` API (meaning that they will be usable in any structure like vectors and N-dimensional arrays).
* Complex numbers backed by the `Field` API (meaning they will be usable in any structure like vectors and
N-dimensional arrays).
* Advanced linear algebra operations like matrix inversion and LU decomposition.
* **Array-like structures** Full support of many-dimensional array-like structures
including mixed arithmetic operations and function operations over arrays and numbers (with the added benefit of static type checking).
* **Expressions** By writing a single mathematical expression
once, users will be able to apply different types of objects to the expression by providing a context. Expressions
can be used for a wide variety of purposes from high performance calculations to code generation.
* **Expressions** By writing a single mathematical expression once, users will be able to apply different types of
objects to the expression by providing a context. Expressions can be used for a wide variety of purposes from high
performance calculations to code generation.
* **Histograms** Fast multi-dimensional histograms.
@ -50,11 +59,10 @@ can be used for a wide variety of purposes from high performance calculations to
* **Type-safe dimensions** Type-safe dimensions for matrix operations.
* **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/)
library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free
to submit a feature request if you want something to be done first.
* **EJML wrapper** Provides EJML `SimpleMatrix` wrapper consistent with the core matrix structures.
* **Commons-math wrapper** It is planned to gradually wrap most parts of
[Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some
parts to better suit the Kotlin programming paradigm, however there is no established roadmap for that. Feel free to
submit a feature request if you want something to be implemented first.
## Planned features
@ -101,7 +109,7 @@ can be used for a wide variety of purposes from high performance calculations to
> - [buffers](kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
> - [expressions](kmath-core/src/commonMain/kotlin/kscience/kmath/expressions) : Functional Expressions
> - [domains](kmath-core/src/commonMain/kotlin/kscience/kmath/domains) : Domains
> - [autodif](kmath-core/src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt) : Automatic differentiation
> - [autodif](kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
<hr/>
@ -117,12 +125,26 @@ can be used for a wide variety of purposes from high performance calculations to
> **Maturity**: EXPERIMENTAL
<hr/>
* ### [kmath-for-real](kmath-for-real)
* ### [kmath-ejml](kmath-ejml)
>
>
> **Maturity**: EXPERIMENTAL
<hr/>
* ### [kmath-for-real](kmath-for-real)
> Extension module that should be used to achieve numpy-like behavior.
All operations are specialized to work with `Double` numbers without declaring algebraic contexts.
One can still use generic algebras though.
>
> **Maturity**: EXPERIMENTAL
>
> **Features:**
> - [RealVector](kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt) : Numpy-like operations for Buffers/Points
> - [RealMatrix](kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt) : Numpy-like operations for 2d real structures
> - [grids](kmath-for-real/src/commonMain/kotlin/kscience/kmath/structures/grids.kt) : Uniform grid generators
<hr/>
* ### [kmath-functions](kmath-functions)
>
>
@ -141,13 +163,31 @@ can be used for a wide variety of purposes from high performance calculations to
> **Maturity**: EXPERIMENTAL
<hr/>
* ### [kmath-kotlingrad](kmath-kotlingrad)
>
>
> **Maturity**: EXPERIMENTAL
<hr/>
* ### [kmath-memory](kmath-memory)
>
>
> **Maturity**: EXPERIMENTAL
<hr/>
* ### [kmath-prob](kmath-prob)
* ### [kmath-nd4j](kmath-nd4j)
> ND4J NDStructure implementation and according NDAlgebra classes
>
> **Maturity**: EXPERIMENTAL
>
> **Features:**
> - [nd4jarraystructure](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
> - [nd4jarrayrings](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Rings over Nd4jArrayStructure of Int and Long
> - [nd4jarrayfields](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : Fields over Nd4jArrayStructure of Float and Double
<hr/>
* ### [kmath-stat](kmath-stat)
>
>
> **Maturity**: EXPERIMENTAL
@ -162,39 +202,53 @@ can be used for a wide variety of purposes from high performance calculations to
## Multi-platform support
KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the [common module](/kmath-core/src/commonMain). Implementation is also done in the common module wherever possible. In some cases, features are delegated to platform-specific implementations even if they could be done in the common module for performance reasons. Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome.
KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the
[common source sets](/kmath-core/src/commonMain) and implemented there wherever it is possible. In some cases, features
are delegated to platform-specific implementations even if they could be provided in the common module for performance
reasons. Currently, the Kotlin/JVM is the primary platform, however Kotlin/Native and Kotlin/JS contributions and
feedback are also welcome.
## Performance
Calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve both performance and flexibility. We expect to focus on creating convenient universal API first and then work on increasing performance for specific cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be better than SciPy.
Calculation performance is one of major goals of KMath in the future, but in some cases it is impossible to achieve
both performance and flexibility.
### Dependency
We expect to focus on creating convenient universal API first and then work on increasing performance for specific
cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized
native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be
better than SciPy.
Release artifacts are accessible from bintray with following configuration (see documentation for [kotlin-multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) form more details):
### Repositories
Release artifacts are accessible from bintray with following configuration (see documentation of
[Kotlin Multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) for more details):
```kotlin
repositories{
repositories {
maven("https://dl.bintray.com/mipt-npm/kscience")
}
dependencies{
api("kscience.kmath:kmath-core:0.2.0-dev-1")
//api("kscience.kmath:kmath-core-jvm:0.2.0-dev-1") for jvm-specific version
dependencies {
api("kscience.kmath:kmath-core:0.2.0-dev-4")
// api("kscience.kmath:kmath-core-jvm:0.2.0-dev-4") for jvm-specific version
}
```
Gradle `6.0+` is required for multiplatform artifacts.
### Development
#### Development
Development builds are uploaded to the separate repository:
Development builds are accessible from the reposirtory
```kotlin
repositories{
repositories {
maven("https://dl.bintray.com/mipt-npm/dev")
}
```
with the same artifact names.
## Contributing
The project requires a lot of additional work. Please feel free to contribute in any way and propose new features.
The project requires a lot of additional work. The most important thing we need is a feedback about what features are
required the most. Feel free to create feature requests. We are also welcome to code contributions,
especially in issues marked with
[waiting for a hero](https://github.com/mipt-npm/kmath/labels/waiting%20for%20a%20hero) label.

28
build.gradle.kts
View File

@ -5,16 +5,23 @@ plugins {
id("ru.mipt.npm.publish") apply false
}
internal val kmathVersion: String by extra("0.2.0-dev-2")
internal val kmathVersion: String by extra("0.2.0-dev-4")
internal val bintrayRepo: String by extra("kscience")
internal val githubProject: String by extra("kmath")
allprojects {
repositories {
jcenter()
maven(url = "https://dl.bintray.com/kotlin/kotlin-eap")
maven(url = "https://dl.bintray.com/kotlin/kotlinx")
maven(url = "https://dl.bintray.com/hotkeytlt/maven")
maven("https://clojars.org/repo")
maven("https://dl.bintray.com/egor-bogomolov/astminer/")
maven("https://dl.bintray.com/hotkeytlt/maven")
maven("https://dl.bintray.com/kotlin/kotlin-eap")
maven("https://dl.bintray.com/kotlin/kotlinx")
maven("https://dl.bintray.com/mipt-npm/dev")
maven("https://dl.bintray.com/mipt-npm/kscience")
maven("https://jitpack.io")
maven("http://logicrunch.research.it.uu.se/maven/")
mavenCentral()
}
group = "kscience.kmath"
@ -22,8 +29,15 @@ allprojects {
}
subprojects {
if (!name.startsWith("kmath")) return@subprojects
apply<KSciencePublishPlugin>()
if (name.startsWith("kmath")) apply<KSciencePublishPlugin>()
}
readme {
readmeTemplate = file("docs/templates/README-TEMPLATE.md")
}
apiValidation {
validationDisabled = true
}
ksciencePublish {
@ -31,5 +45,3 @@ ksciencePublish {
spaceUser = System.getenv("SPACE_USER")
spaceToken = System.getenv("SPACE_TOKEN")
}
readme.readmeTemplate = file("docs/templates/README-TEMPLATE.md")

78
docs/templates/README-TEMPLATE.md vendored
View File

@ -8,41 +8,50 @@ Bintray: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience
Bintray-dev: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion)
# KMath
Could be pronounced as `key-math`.
The Kotlin MATHematics library was initially intended as a Kotlin-based analog to Python's `numpy` library. Later we found that kotlin is much more flexible language and allows superior architecture designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could be achieved with [kmath-for-real](/kmath-for-real) extension module.
Could be pronounced as `key-math`. The Kotlin MATHematics library was initially intended as a Kotlin-based analog to
Python's NumPy library. Later we found that kotlin is much more flexible language and allows superior architecture
designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could
be achieved with [kmath-for-real](/kmath-for-real) extension module.
## Publications and talks
* [A conceptual article about context-oriented design](https://proandroiddev.com/an-introduction-context-oriented-programming-in-kotlin-2e79d316b0a2)
* [Another article about context-oriented design](https://proandroiddev.com/diving-deeper-into-context-oriented-programming-in-kotlin-3ecb4ec38814)
* [ACAT 2019 conference paper](https://aip.scitation.org/doi/abs/10.1063/1.5130103)
# Goal
* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM and JS for now and Native in future).
* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS and Native).
* Provide basic multiplatform implementations for those abstractions (without significant performance optimization).
* Provide bindings and wrappers with those abstractions for popular optimized platform libraries.
## Non-goals
* Be like Numpy. It was the idea at the beginning, but we decided that we can do better in terms of API.
* Provide best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
* Be like NumPy. It was the idea at the beginning, but we decided that we can do better in terms of API.
* Provide the best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
* Cover all cases as immediately and in one bundle. We will modularize everything and add new features gradually.
* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better experience for those, who want to work with specific types.
* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like
for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better
experience for those, who want to work with specific types.
## Features
Actual feature list is [here](/docs/features.md)
Current feature list is [here](/docs/features.md)
* **Algebra**
* Algebraic structures like rings, spaces and field (**TODO** add example to wiki)
* Algebraic structures like rings, spaces and fields (**TODO** add example to wiki)
* Basic linear algebra operations (sums, products, etc.), backed by the `Space` API.
* Complex numbers backed by the `Field` API (meaning that they will be usable in any structure like vectors and N-dimensional arrays).
* Complex numbers backed by the `Field` API (meaning they will be usable in any structure like vectors and
N-dimensional arrays).
* Advanced linear algebra operations like matrix inversion and LU decomposition.
* **Array-like structures** Full support of many-dimensional array-like structures
including mixed arithmetic operations and function operations over arrays and numbers (with the added benefit of static type checking).
* **Expressions** By writing a single mathematical expression
once, users will be able to apply different types of objects to the expression by providing a context. Expressions
can be used for a wide variety of purposes from high performance calculations to code generation.
* **Expressions** By writing a single mathematical expression once, users will be able to apply different types of
objects to the expression by providing a context. Expressions can be used for a wide variety of purposes from high
performance calculations to code generation.
* **Histograms** Fast multi-dimensional histograms.
@ -50,9 +59,10 @@ can be used for a wide variety of purposes from high performance calculations to
* **Type-safe dimensions** Type-safe dimensions for matrix operations.
* **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/)
library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free
to submit a feature request if you want something to be done first.
* **Commons-math wrapper** It is planned to gradually wrap most parts of
[Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some
parts to better suit the Kotlin programming paradigm, however there is no established roadmap for that. Feel free to
submit a feature request if you want something to be implemented first.
## Planned features
@ -72,39 +82,53 @@ $modules
## Multi-platform support
KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the [common module](/kmath-core/src/commonMain). Implementation is also done in the common module wherever possible. In some cases, features are delegated to platform-specific implementations even if they could be done in the common module for performance reasons. Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome.
KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the
[common source sets](/kmath-core/src/commonMain) and implemented there wherever it is possible. In some cases, features
are delegated to platform-specific implementations even if they could be provided in the common module for performance
reasons. Currently, the Kotlin/JVM is the primary platform, however Kotlin/Native and Kotlin/JS contributions and
feedback are also welcome.
## Performance
Calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve both performance and flexibility. We expect to focus on creating convenient universal API first and then work on increasing performance for specific cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be better than SciPy.
Calculation performance is one of major goals of KMath in the future, but in some cases it is impossible to achieve
both performance and flexibility.
### Dependency
We expect to focus on creating convenient universal API first and then work on increasing performance for specific
cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized
native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be
better than SciPy.
Release artifacts are accessible from bintray with following configuration (see documentation for [kotlin-multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) form more details):
### Repositories
Release artifacts are accessible from bintray with following configuration (see documentation of
[Kotlin Multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) for more details):
```kotlin
repositories{
repositories {
maven("https://dl.bintray.com/mipt-npm/kscience")
}
dependencies{
dependencies {
api("kscience.kmath:kmath-core:$version")
//api("kscience.kmath:kmath-core-jvm:$version") for jvm-specific version
// api("kscience.kmath:kmath-core-jvm:$version") for jvm-specific version
}
```
Gradle `6.0+` is required for multiplatform artifacts.
### Development
#### Development
Development builds are uploaded to the separate repository:
Development builds are accessible from the reposirtory
```kotlin
repositories{
repositories {
maven("https://dl.bintray.com/mipt-npm/dev")
}
```
with the same artifact names.
## Contributing
The project requires a lot of additional work. Please feel free to contribute in any way and propose new features.
The project requires a lot of additional work. The most important thing we need is a feedback about what features are
required the most. Feel free to create feature requests. We are also welcome to code contributions,
especially in issues marked with
[waiting for a hero](https://github.com/mipt-npm/kmath/labels/waiting%20for%20a%20hero) label.

45
examples/build.gradle.kts
View File

@ -8,30 +8,57 @@ plugins {
}
allOpen.annotation("org.openjdk.jmh.annotations.State")
sourceSets.register("benchmarks")
repositories {
jcenter()
maven("https://dl.bintray.com/kotlin/kotlin-eap/")
maven("https://clojars.org/repo")
maven("https://dl.bintray.com/egor-bogomolov/astminer/")
maven("https://dl.bintray.com/hotkeytlt/maven")
maven("https://dl.bintray.com/kotlin/kotlin-eap")
maven("https://dl.bintray.com/kotlin/kotlinx")
maven("https://dl.bintray.com/mipt-npm/dev")
maven("https://dl.bintray.com/mipt-npm/kscience")
maven("https://jitpack.io")
maven("http://logicrunch.research.it.uu.se/maven/")
mavenCentral()
}
sourceSets.register("benchmarks")
dependencies {
// implementation(project(":kmath-ast"))
implementation(project(":kmath-ast"))
implementation(project(":kmath-kotlingrad"))
implementation(project(":kmath-core"))
implementation(project(":kmath-coroutines"))
implementation(project(":kmath-commons"))
implementation(project(":kmath-prob"))
implementation(project(":kmath-stat"))
implementation(project(":kmath-viktor"))
implementation(project(":kmath-dimensions"))
implementation(project(":kmath-ejml"))
implementation(project(":kmath-nd4j"))
implementation(project(":kmath-for-real"))
implementation("org.deeplearning4j:deeplearning4j-core:1.0.0-beta7")
implementation("org.nd4j:nd4j-native:1.0.0-beta7")
// uncomment if your system supports AVX2
// val os = System.getProperty("os.name")
//
// if (System.getProperty("os.arch") in arrayOf("x86_64", "amd64")) when {
// os.startsWith("Windows") -> implementation("org.nd4j:nd4j-native:1.0.0-beta7:windows-x86_64-avx2")
// os == "Linux" -> implementation("org.nd4j:nd4j-native:1.0.0-beta7:linux-x86_64-avx2")
// os == "Mac OS X" -> implementation("org.nd4j:nd4j-native:1.0.0-beta7:macosx-x86_64-avx2")
// } else
implementation("org.nd4j:nd4j-native-platform:1.0.0-beta7")
implementation("org.jetbrains.kotlinx:kotlinx-io:0.2.0-npm-dev-11")
implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-20")
implementation("org.slf4j:slf4j-simple:1.7.30")
"benchmarksImplementation"("org.jetbrains.kotlinx:kotlinx.benchmark.runtime-jvm:0.2.0-dev-8")
// plotting
implementation("kscience.plotlykt:plotlykt-server:0.3.1-dev")
"benchmarksImplementation"("org.jetbrains.kotlinx:kotlinx.benchmark.runtime-jvm:0.2.0-dev-20")
"benchmarksImplementation"(sourceSets.main.get().output + sourceSets.main.get().runtimeClasspath)
}
@ -42,7 +69,7 @@ benchmark {
// This one matches sourceSet name above
configurations.register("fast") {
warmups = 5 // number of warmup iterations
warmups = 1 // number of warmup iterations
iterations = 3 // number of iterations
iterationTime = 500 // time in seconds per iteration
iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
@ -56,4 +83,6 @@ kotlin.sourceSets.all {
}
}
tasks.withType<KotlinCompile> { kotlinOptions.jvmTarget = "11" }
tasks.withType<KotlinCompile> {
kotlinOptions.jvmTarget = "11"
}

2
examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt → examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt
View File

@ -1,4 +1,4 @@
package kscience.kmath.structures
package kscience.kmath.benchmarks
import org.openjdk.jmh.annotations.Benchmark
import org.openjdk.jmh.annotations.Scope

4
examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt → examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt
View File

@ -1,7 +1,9 @@
package kscience.kmath.structures
package kscience.kmath.benchmarks
import kscience.kmath.operations.Complex
import kscience.kmath.operations.complex
import kscience.kmath.structures.MutableBuffer
import kscience.kmath.structures.RealBuffer
import org.openjdk.jmh.annotations.Benchmark
import org.openjdk.jmh.annotations.Scope
import org.openjdk.jmh.annotations.State

50
examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt Normal file
View File

@ -0,0 +1,50 @@
package kscience.kmath.linear
import kotlinx.benchmark.Benchmark
import kscience.kmath.commons.linear.CMMatrixContext
import kscience.kmath.commons.linear.CMMatrixContext.dot
import kscience.kmath.commons.linear.inverse
import kscience.kmath.commons.linear.toCM
import kscience.kmath.ejml.EjmlMatrixContext
import kscience.kmath.ejml.inverse
import kscience.kmath.ejml.toEjml
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
import org.openjdk.jmh.annotations.Scope
import org.openjdk.jmh.annotations.State
import kotlin.random.Random
@State(Scope.Benchmark)
class LinearAlgebraBenchmark {
companion object {
val random = Random(1224)
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
}
@Benchmark
fun kmathLUPInversion() {
MatrixContext.real.inverseWithLUP(matrix)
}
@Benchmark
fun cmLUPInversion() {
CMMatrixContext {
val cm = matrix.toCM() //avoid overhead on conversion
inverse(cm)
}
}
@Benchmark
fun ejmlInverse() {
EjmlMatrixContext {
val km = matrix.toEjml() //avoid overhead on conversion
inverse(km)
}
}
}

60
examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt Normal file
View File

@ -0,0 +1,60 @@
package kscience.kmath.benchmarks
import kotlinx.benchmark.Benchmark
import kscience.kmath.commons.linear.CMMatrixContext
import kscience.kmath.commons.linear.CMMatrixContext.dot
import kscience.kmath.commons.linear.toCM
import kscience.kmath.ejml.EjmlMatrixContext
import kscience.kmath.ejml.toEjml
import kscience.kmath.linear.real
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
import org.openjdk.jmh.annotations.Scope
import org.openjdk.jmh.annotations.State
import kotlin.random.Random
@State(Scope.Benchmark)
class MultiplicationBenchmark {
companion object {
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 }
val cmMatrix1 = matrix1.toCM()
val cmMatrix2 = matrix2.toCM()
val ejmlMatrix1 = matrix1.toEjml()
val ejmlMatrix2 = matrix2.toEjml()
}
@Benchmark
fun commonsMathMultiplication() {
CMMatrixContext.invoke {
cmMatrix1 dot cmMatrix2
}
}
@Benchmark
fun ejmlMultiplication() {
EjmlMatrixContext.invoke {
ejmlMatrix1 dot ejmlMatrix2
}
}
@Benchmark
fun ejmlMultiplicationwithConversion() {
val ejmlMatrix1 = matrix1.toEjml()
val ejmlMatrix2 = matrix2.toEjml()
EjmlMatrixContext.invoke {
ejmlMatrix1 dot ejmlMatrix2
}
}
@Benchmark
fun bufferedMultiplication() {
matrix1 dot matrix2
}
}

3
examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt → examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
View File

@ -1,7 +1,8 @@
package kscience.kmath.structures
package kscience.kmath.benchmarks
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.structures.*
import org.openjdk.jmh.annotations.Benchmark
import org.openjdk.jmh.annotations.Scope
import org.openjdk.jmh.annotations.State

7
examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt → examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
View File

@ -1,7 +1,10 @@
package kscience.kmath.structures
package kscience.kmath.benchmarks
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.structures.BufferedNDField
import kscience.kmath.structures.NDField
import kscience.kmath.structures.RealNDField
import kscience.kmath.viktor.ViktorNDField
import org.jetbrains.bio.viktor.F64Array
import org.openjdk.jmh.annotations.Benchmark
@ -36,7 +39,7 @@ internal class ViktorBenchmark {
@Benchmark
fun rawViktor() {
val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim))
val one = F64Array.full(init = 1.0, shape = intArrayOf(dim, dim))
var res = one
repeat(n) { res = res + one }
}

11
examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt
View File

@ -1,11 +0,0 @@
package kscience.kmath.utils
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
import kotlin.system.measureTimeMillis
internal inline fun measureAndPrint(title: String, block: () -> Unit) {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
val time = measureTimeMillis(block)
println("$title completed in $time millis")
}

148
examples/src/main/kotlin/kscience/kmath/ast/ExpressionsInterpretersBenchmark.kt
View File

@ -1,70 +1,80 @@
package kscience.kmath.ast
//
//import kscience.kmath.asm.compile
//import kscience.kmath.expressions.Expression
//import kscience.kmath.expressions.expressionInField
//import kscience.kmath.expressions.invoke
//import kscience.kmath.operations.Field
//import kscience.kmath.operations.RealField
//import kotlin.random.Random
//import kotlin.system.measureTimeMillis
//
//class ExpressionsInterpretersBenchmark {
// private val algebra: Field<Double> = RealField
// fun functionalExpression() {
// val expr = algebra.expressionInField {
// variable("x") * const(2.0) + const(2.0) / variable("x") - const(16.0)
// }
//
// invokeAndSum(expr)
// }
//
// fun mstExpression() {
// val expr = algebra.mstInField {
// symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
// }
//
// invokeAndSum(expr)
// }
//
// fun asmExpression() {
// val expr = algebra.mstInField {
// symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
// }.compile()
//
// invokeAndSum(expr)
// }
//
// private fun invokeAndSum(expr: Expression<Double>) {
// val random = Random(0)
// var sum = 0.0
//
// repeat(1000000) {
// sum += expr("x" to random.nextDouble())
// }
//
// println(sum)
// }
//}
//
//fun main() {
// val benchmark = ExpressionsInterpretersBenchmark()
//
// val fe = measureTimeMillis {
// benchmark.functionalExpression()
// }
//
// println("fe=$fe")
//
// val mst = measureTimeMillis {
// benchmark.mstExpression()
// }
//
// println("mst=$mst")
//
// val asm = measureTimeMillis {
// benchmark.asmExpression()
// }
//
// println("asm=$asm")
//}
import kscience.kmath.asm.compile
import kscience.kmath.expressions.Expression
import kscience.kmath.expressions.expressionInField
import kscience.kmath.expressions.invoke
import kscience.kmath.operations.Field
import kscience.kmath.operations.RealField
import kotlin.random.Random
import kotlin.system.measureTimeMillis
internal class ExpressionsInterpretersBenchmark {
private val algebra: Field<Double> = RealField
fun functionalExpression() {
val expr = algebra.expressionInField {
symbol("x") * const(2.0) + const(2.0) / symbol("x") - const(16.0)
}
invokeAndSum(expr)
}
fun mstExpression() {
val expr = algebra.mstInField {
symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
}
invokeAndSum(expr)
}
fun asmExpression() {
val expr = algebra.mstInField {
symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
}.compile()
invokeAndSum(expr)
}
private fun invokeAndSum(expr: Expression<Double>) {
val random = Random(0)
var sum = 0.0
repeat(1000000) {
sum += expr("x" to random.nextDouble())
}
println(sum)
}
}
/**
* This benchmark compares basically evaluation of simple function with MstExpression interpreter, ASM backend and
* core FunctionalExpressions API.
*
* The expected rating is:
*
* 1. ASM.
* 2. MST.
* 3. FE.
*/
fun main() {
val benchmark = ExpressionsInterpretersBenchmark()
val fe = measureTimeMillis {
benchmark.functionalExpression()
}
println("fe=$fe")
val mst = measureTimeMillis {
benchmark.mstExpression()
}
println("mst=$mst")
val asm = measureTimeMillis {
benchmark.asmExpression()
}
println("asm=$asm")
}

24
examples/src/main/kotlin/kscience/kmath/ast/KotlingradSupport.kt Normal file
View File

@ -0,0 +1,24 @@
package kscience.kmath.ast
import kscience.kmath.asm.compile
import kscience.kmath.expressions.derivative
import kscience.kmath.expressions.invoke
import kscience.kmath.expressions.symbol
import kscience.kmath.kotlingrad.differentiable
import kscience.kmath.operations.RealField
/**
* In this example, x^2-4*x-44 function is differentiated with Kotlin∇, and the autodiff result is compared with
* valid derivative.
*/
fun main() {
val x by symbol
val actualDerivative = MstExpression(RealField, "x^2-4*x-44".parseMath())
.differentiable()
.derivative(x)
.compile()
val expectedDerivative = MstExpression(RealField, "2*x-4".parseMath()).compile()
assert(actualDerivative("x" to 123.0) == expectedDerivative("x" to 123.0))
}

102
examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt Normal file
View File

@ -0,0 +1,102 @@
package kscience.kmath.commons.fit
import kotlinx.html.br
import kotlinx.html.h3
import kscience.kmath.commons.optimization.chiSquared
import kscience.kmath.commons.optimization.minimize
import kscience.kmath.expressions.symbol
import kscience.kmath.real.RealVector
import kscience.kmath.real.map
import kscience.kmath.real.step
import kscience.kmath.stat.*
import kscience.kmath.structures.asIterable
import kscience.kmath.structures.toList
import kscience.plotly.*
import kscience.plotly.models.ScatterMode
import kscience.plotly.models.TraceValues
import kotlin.math.pow
import kotlin.math.sqrt
//Forward declaration of symbols that will be used in expressions.
// This declaration is required for
private val a by symbol
private val b by symbol
private val c by symbol
/**
* Shortcut to use buffers in plotly
*/
operator fun TraceValues.invoke(vector: RealVector) {
numbers = vector.asIterable()
}
/**
* Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules.
*/
fun main() {
//A generator for a normally distributed values
val generator = Distribution.normal()
//A chain/flow of random values with the given seed
val chain = generator.sample(RandomGenerator.default(112667))
//Create a uniformly distributed x values like numpy.arrange
val x = 1.0..100.0 step 1.0
//Perform an operation on each x value (much more effective, than numpy)
val y = x.map {
val value = it.pow(2) + it + 1
value + chain.nextDouble() * sqrt(value)
}
// this will also work, but less effective:
// val y = x.pow(2)+ x + 1 + chain.nextDouble()
// create same errors for all xs
val yErr = y.map { sqrt(it) }//RealVector.same(x.size, sigma)
// compute differentiable chi^2 sum for given model ax^2 + bx + c
val chi2 = Fitting.chiSquared(x, y, yErr) { x1 ->
//bind variables to autodiff context
val a = bind(a)
val b = bind(b)
//Include default value for c if it is not provided as a parameter
val c = bindOrNull(c) ?: one
a * x1.pow(2) + b * x1 + c
}
//minimize the chi^2 in given starting point. Derivatives are not required, they are already included.
val result: OptimizationResult<Double> = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
//display a page with plot and numerical results
val page = Plotly.page {
plot {
scatter {
mode = ScatterMode.markers
x(x)
y(y)
error_y {
array = yErr.toList()
}
name = "data"
}
scatter {
mode = ScatterMode.lines
x(x)
y(x.map { result.point[a]!! * it.pow(2) + result.point[b]!! * it + 1 })
name = "fit"
}
}
br()
h3{
+"Fit result: $result"
}
h3{
+"Chi2/dof = ${result.value / (x.size - 3)}"
}
}
page.makeFile()
}

50
examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt
View File

@ -1,50 +0,0 @@
package kscience.kmath.linear
import kscience.kmath.commons.linear.CMMatrixContext
import kscience.kmath.commons.linear.inverse
import kscience.kmath.commons.linear.toCM
import kscience.kmath.ejml.EjmlMatrixContext
import kscience.kmath.ejml.inverse
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
import kotlin.random.Random
import kotlin.system.measureTimeMillis
fun main() {
val random = Random(1224)
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 = 5000 // iterations
MatrixContext.real {
repeat(50) { inverse(matrix) }
val inverseTime = measureTimeMillis { repeat(n) { inverse(matrix) } }
println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis")
}
//commons-math
val commonsTime = measureTimeMillis {
CMMatrixContext {
val cm = matrix.toCM() //avoid overhead on conversion
repeat(n) { inverse(cm) }
}
}
println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis")
val ejmlTime = measureTimeMillis {
(EjmlMatrixContext(RealField)) {
val km = matrix.toEjml() //avoid overhead on conversion
repeat(n) { inverse(km) }
}
}
println("[ejml] Inversion of $n matrices $dim x $dim finished in $ejmlTime millis")
}

38
examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt
View File

@ -1,38 +0,0 @@
package kscience.kmath.linear
import kscience.kmath.commons.linear.CMMatrixContext
import kscience.kmath.commons.linear.toCM
import kscience.kmath.ejml.EjmlMatrixContext
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
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 {
val cmMatrix1 = matrix1.toCM()
val cmMatrix2 = matrix2.toCM()
val cmTime = measureTimeMillis { cmMatrix1 dot cmMatrix2 }
println("CM implementation time: $cmTime")
}
(EjmlMatrixContext(RealField)) {
val ejmlMatrix1 = matrix1.toEjml()
val ejmlMatrix2 = matrix2.toEjml()
val ejmlTime = measureTimeMillis { ejmlMatrix1 dot ejmlMatrix2 }
println("EJML implementation time: $ejmlTime")
}
val genericTime = measureTimeMillis { val res = matrix1 dot matrix2 }
println("Generic implementation time: $genericTime")
}

4
examples/src/main/kotlin/kscience/kmath/operations/ComplexDemo.kt
View File

@ -6,8 +6,8 @@ import kscience.kmath.structures.complex
fun main() {
// 2d element
val element = NDElement.complex(2, 2) { index: IntArray ->
Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble())
val element = NDElement.complex(2, 2) { (i,j) ->
Complex(i.toDouble() - j.toDouble(), i.toDouble() + j.toDouble())
}
println(element)

2
examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionBenchmark.kt → examples/src/main/kotlin/kscience/kmath/stat/DistributionBenchmark.kt
View File

@ -4,7 +4,7 @@ import kotlinx.coroutines.Dispatchers
import kotlinx.coroutines.async
import kotlinx.coroutines.runBlocking
import kscience.kmath.chains.BlockingRealChain
import kscience.kmath.prob.*
import kscience.kmath.stat.*
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler
import org.apache.commons.rng.simple.RandomSource
import java.time.Duration

11
examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionDemo.kt → examples/src/main/kotlin/kscience/kmath/stat/DistributionDemo.kt
View File

@ -1,14 +1,17 @@
package kscience.kmath.commons.prob
package kscience.kmath.stat
import kotlinx.coroutines.runBlocking
import kscience.kmath.chains.Chain
import kscience.kmath.chains.collectWithState
import kscience.kmath.prob.Distribution
import kscience.kmath.prob.RandomGenerator
import kscience.kmath.prob.normal
/**
* The state of distribution averager
*/
private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
/**
* Averaging
*/
private fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChainState(), { it.copy() }) { chain ->
val next = chain.next()
num++

13
examples/src/main/kotlin/kscience/kmath/structures/NDField.kt
View File

@ -1,8 +1,10 @@
package kscience.kmath.structures
import kotlinx.coroutines.GlobalScope
import kscience.kmath.nd4j.Nd4jArrayField
import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import org.nd4j.linalg.factory.Nd4j
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
import kotlin.system.measureTimeMillis
@ -14,6 +16,8 @@ internal inline fun measureAndPrint(title: String, block: () -> Unit) {
}
fun main() {
// initializing Nd4j
Nd4j.zeros(0)
val dim = 1000
val n = 1000
@ -23,6 +27,8 @@ fun main() {
val specializedField = NDField.real(dim, dim)
//A generic boxing field. It should be used for objects, not primitives.
val genericField = NDField.boxing(RealField, dim, dim)
// Nd4j specialized field.
val nd4jField = Nd4jArrayField.real(dim, dim)
measureAndPrint("Automatic field addition") {
autoField {
@ -43,6 +49,13 @@ fun main() {
}
}
measureAndPrint("Nd4j specialized addition") {
nd4jField {
var res = one
repeat(n) { res += 1.0 as Number }
}
}
measureAndPrint("Lazy addition") {
val res = specializedField.one.mapAsync(GlobalScope) {
var c = 0.0

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

@ -1,5 +1,5 @@
distributionBase=GRADLE_USER_HOME
distributionPath=wrapper/dists
distributionUrl=https\://services.gradle.org/distributions/gradle-6.6.1-bin.zip
distributionUrl=https\://services.gradle.org/distributions/gradle-6.7-bin.zip
zipStoreBase=GRADLE_USER_HOME
zipStorePath=wrapper/dists

4
kmath-ast/build.gradle.kts
View File

@ -17,3 +17,7 @@ kotlin.sourceSets {
}
}
}
readme{
maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
}

113
kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstAlgebra.kt
View File

@ -6,14 +6,14 @@ import kscience.kmath.operations.*
* [Algebra] over [MST] nodes.
*/
public object MstAlgebra : NumericAlgebra<MST> {
override fun number(value: Number): MST = MST.Numeric(value)
override fun number(value: Number): MST.Numeric = MST.Numeric(value)
override fun symbol(value: String): MST = MST.Symbolic(value)
override fun symbol(value: String): MST.Symbolic = MST.Symbolic(value)
override fun unaryOperation(operation: String, arg: MST): MST =
override fun unaryOperation(operation: String, arg: MST): MST.Unary =
MST.Unary(operation, arg)
override fun binaryOperation(operation: String, left: MST, right: MST): MST =
override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
MST.Binary(operation, left, right)
}
@ -21,97 +21,100 @@ public object MstAlgebra : NumericAlgebra<MST> {
* [Space] over [MST] nodes.
*/
public object MstSpace : Space<MST>, NumericAlgebra<MST> {
override val zero: MST = number(0.0)
override val zero: MST.Numeric by lazy { number(0.0) }
override fun number(value: Number): MST = MstAlgebra.number(value)
override fun symbol(value: String): MST = MstAlgebra.symbol(value)
override fun add(a: MST, b: MST): MST = binaryOperation(SpaceOperations.PLUS_OPERATION, a, b)
override fun multiply(a: MST, k: Number): MST = binaryOperation(RingOperations.TIMES_OPERATION, a, number(k))
override fun number(value: Number): MST.Numeric = MstAlgebra.number(value)
override fun symbol(value: String): MST.Symbolic = MstAlgebra.symbol(value)
override fun add(a: MST, b: MST): MST.Binary = binaryOperation(SpaceOperations.PLUS_OPERATION, a, b)
override fun multiply(a: MST, k: Number): MST.Binary = binaryOperation(RingOperations.TIMES_OPERATION, a, number(k))
override fun binaryOperation(operation: String, left: MST, right: MST): MST =
override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
MstAlgebra.binaryOperation(operation, left, right)
override fun unaryOperation(operation: String, arg: MST): MST = MstAlgebra.unaryOperation(operation, arg)
override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstAlgebra.unaryOperation(operation, arg)
}
/**
* [Ring] over [MST] nodes.
*/
public object MstRing : Ring<MST>, NumericAlgebra<MST> {
override val zero: MST
override val zero: MST.Numeric
get() = MstSpace.zero
override val one: MST = number(1.0)
override fun number(value: Number): MST = MstSpace.number(value)
override fun symbol(value: String): MST = MstSpace.symbol(value)
override fun add(a: MST, b: MST): MST = MstSpace.add(a, b)
override val one: MST.Numeric by lazy { number(1.0) }
override fun multiply(a: MST, k: Number): MST = MstSpace.multiply(a, k)
override fun number(value: Number): MST.Numeric = MstSpace.number(value)
override fun symbol(value: String): MST.Symbolic = MstSpace.symbol(value)
override fun add(a: MST, b: MST): MST.Binary = MstSpace.add(a, b)
override fun multiply(a: MST, k: Number): MST.Binary = MstSpace.multiply(a, k)
override fun multiply(a: MST, b: MST): MST.Binary = binaryOperation(RingOperations.TIMES_OPERATION, a, b)
override fun multiply(a: MST, b: MST): MST = binaryOperation(RingOperations.TIMES_OPERATION, a, b)
override fun binaryOperation(operation: String, left: MST, right: MST): MST =
override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
MstSpace.binaryOperation(operation, left, right)
override fun unaryOperation(operation: String, arg: MST): MST = MstAlgebra.unaryOperation(operation, arg)
override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstSpace.unaryOperation(operation, arg)
}
/**
* [Field] over [MST] nodes.
*/
public object MstField : Field<MST> {
public override val zero: MST
public override val zero: MST.Numeric
get() = MstRing.zero
public override val one: MST
public override val one: MST.Numeric
get() = MstRing.one
public override fun symbol(value: String): MST = MstRing.symbol(value)
public override fun number(value: Number): MST = MstRing.number(value)
public override fun add(a: MST, b: MST): MST = MstRing.add(a, b)
public override fun multiply(a: MST, k: Number): MST = MstRing.multiply(a, k)
public override fun multiply(a: MST, b: MST): MST = MstRing.multiply(a, b)
public override fun divide(a: MST, b: MST): MST = binaryOperation(FieldOperations.DIV_OPERATION, a, b)
public override fun symbol(value: String): MST.Symbolic = MstRing.symbol(value)
public override fun number(value: Number): MST.Numeric = MstRing.number(value)
public override fun add(a: MST, b: MST): MST.Binary = MstRing.add(a, b)
public override fun multiply(a: MST, k: Number): MST.Binary = MstRing.multiply(a, k)
public override fun multiply(a: MST, b: MST): MST.Binary = MstRing.multiply(a, b)
public override fun divide(a: MST, b: MST): MST.Binary = binaryOperation(FieldOperations.DIV_OPERATION, a, b)
public override fun binaryOperation(operation: String, left: MST, right: MST): MST =
public override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
MstRing.binaryOperation(operation, left, right)
override fun unaryOperation(operation: String, arg: MST): MST = MstRing.unaryOperation(operation, arg)
override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstRing.unaryOperation(operation, arg)
}
/**
* [ExtendedField] over [MST] nodes.
*/
public object MstExtendedField : ExtendedField<MST> {
override val zero: MST
override val zero: MST.Numeric
get() = MstField.zero
override val one: MST
override val one: MST.Numeric
get() = MstField.one
override fun symbol(value: String): MST = MstField.symbol(value)
override fun sin(arg: MST): MST = unaryOperation(TrigonometricOperations.SIN_OPERATION, arg)
override fun cos(arg: MST): MST = unaryOperation(TrigonometricOperations.COS_OPERATION, arg)
override fun tan(arg: MST): MST = unaryOperation(TrigonometricOperations.TAN_OPERATION, arg)
override fun asin(arg: MST): MST = unaryOperation(TrigonometricOperations.ASIN_OPERATION, arg)
override fun acos(arg: MST): MST = unaryOperation(TrigonometricOperations.ACOS_OPERATION, arg)
override fun atan(arg: MST): MST = unaryOperation(TrigonometricOperations.ATAN_OPERATION, arg)
override fun sinh(arg: MST): MST = unaryOperation(HyperbolicOperations.SINH_OPERATION, arg)
override fun cosh(arg: MST): MST = unaryOperation(HyperbolicOperations.COSH_OPERATION, arg)
override fun tanh(arg: MST): MST = unaryOperation(HyperbolicOperations.TANH_OPERATION, arg)
override fun asinh(arg: MST): MST = unaryOperation(HyperbolicOperations.ASINH_OPERATION, arg)
override fun acosh(arg: MST): MST = unaryOperation(HyperbolicOperations.ACOSH_OPERATION, arg)
override fun atanh(arg: MST): MST = unaryOperation(HyperbolicOperations.ATANH_OPERATION, arg)
override fun add(a: MST, b: MST): MST = MstField.add(a, b)
override fun multiply(a: MST, k: Number): MST = MstField.multiply(a, k)
override fun multiply(a: MST, b: MST): MST = MstField.multiply(a, b)
override fun divide(a: MST, b: MST): MST = MstField.divide(a, b)
override fun power(arg: MST, pow: Number): MST = binaryOperation(PowerOperations.POW_OPERATION, arg, number(pow))
override fun exp(arg: MST): MST = unaryOperation(ExponentialOperations.EXP_OPERATION, arg)
override fun ln(arg: MST): MST = unaryOperation(ExponentialOperations.LN_OPERATION, arg)
override fun symbol(value: String): MST.Symbolic = MstField.symbol(value)
override fun number(value: Number): MST.Numeric = MstField.number(value)
override fun sin(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.SIN_OPERATION, arg)
override fun cos(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.COS_OPERATION, arg)
override fun tan(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.TAN_OPERATION, arg)
override fun asin(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.ASIN_OPERATION, arg)
override fun acos(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.ACOS_OPERATION, arg)
override fun atan(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.ATAN_OPERATION, arg)
override fun sinh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.SINH_OPERATION, arg)
override fun cosh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.COSH_OPERATION, arg)
override fun tanh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.TANH_OPERATION, arg)
override fun asinh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.ASINH_OPERATION, arg)
override fun acosh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.ACOSH_OPERATION, arg)
override fun atanh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.ATANH_OPERATION, arg)
override fun add(a: MST, b: MST): MST.Binary = MstField.add(a, b)
override fun multiply(a: MST, k: Number): MST.Binary = MstField.multiply(a, k)
override fun multiply(a: MST, b: MST): MST.Binary = MstField.multiply(a, b)
override fun divide(a: MST, b: MST): MST.Binary = MstField.divide(a, b)
override fun binaryOperation(operation: String, left: MST, right: MST): MST =
override fun power(arg: MST, pow: Number): MST.Binary =
binaryOperation(PowerOperations.POW_OPERATION, arg, number(pow))
override fun exp(arg: MST): MST.Unary = unaryOperation(ExponentialOperations.EXP_OPERATION, arg)
override fun ln(arg: MST): MST.Unary = unaryOperation(ExponentialOperations.LN_OPERATION, arg)
override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
MstField.binaryOperation(operation, left, right)
override fun unaryOperation(operation: String, arg: MST): MST = MstField.unaryOperation(operation, arg)
override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstField.unaryOperation(operation, arg)
}

35
kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstExpression.kt
View File

@ -13,21 +13,22 @@ import kotlin.contracts.contract
* @property mst the [MST] node.
* @author Alexander Nozik
*/
public class MstExpression<T>(public val algebra: Algebra<T>, public val mst: MST) : Expression<T> {
private inner class InnerAlgebra(val arguments: Map<String, T>) : NumericAlgebra<T> {
override fun symbol(value: String): T = arguments[value] ?: algebra.symbol(value)
public class MstExpression<T, out A : Algebra<T>>(public val algebra: A, public val mst: MST) : Expression<T> {
private inner class InnerAlgebra(val arguments: Map<Symbol, T>) : NumericAlgebra<T> {
override fun symbol(value: String): T = arguments[StringSymbol(value)] ?: algebra.symbol(value)
override fun unaryOperation(operation: String, arg: T): T = algebra.unaryOperation(operation, arg)
override fun binaryOperation(operation: String, left: T, right: T): T =
algebra.binaryOperation(operation, left, right)
override fun number(value: Number): T = if (algebra is NumericAlgebra)
algebra.number(value)
@Suppress("UNCHECKED_CAST")
override fun number(value: Number): T = if (algebra is NumericAlgebra<*>)
(algebra as NumericAlgebra<T>).number(value)
else
error("Numeric nodes are not supported by $this")
}
override operator fun invoke(arguments: Map<String, T>): T = InnerAlgebra(arguments).evaluate(mst)
override operator fun invoke(arguments: Map<Symbol, T>): T = InnerAlgebra(arguments).evaluate(mst)
}
/**
@ -37,15 +38,15 @@ public class MstExpression<T>(public val algebra: Algebra<T>, public val mst: MS
*/
public inline fun <reified T : Any, A : Algebra<T>, E : Algebra<MST>> A.mst(
mstAlgebra: E,
block: E.() -> MST
): MstExpression<T> = MstExpression(this, mstAlgebra.block())
block: E.() -> MST,
): MstExpression<T, A> = MstExpression(this, mstAlgebra.block())
/**
* Builds [MstExpression] over [Space].
*
* @author Alexander Nozik
*/
public inline fun <reified T : Any> Space<T>.mstInSpace(block: MstSpace.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : Space<T>> A.mstInSpace(block: MstSpace.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return MstExpression(this, MstSpace.block())
}
@ -55,7 +56,7 @@ public inline fun <reified T : Any> Space<T>.mstInSpace(block: MstSpace.() -> MS
*
* @author Alexander Nozik
*/
public inline fun <reified T : Any> Ring<T>.mstInRing(block: MstRing.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : Ring<T>> A.mstInRing(block: MstRing.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return MstExpression(this, MstRing.block())
}
@ -65,7 +66,7 @@ public inline fun <reified T : Any> Ring<T>.mstInRing(block: MstRing.() -> MST):
*
* @author Alexander Nozik
*/
public inline fun <reified T : Any> Field<T>.mstInField(block: MstField.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : Field<T>> A.mstInField(block: MstField.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return MstExpression(this, MstField.block())
}
@ -75,7 +76,7 @@ public inline fun <reified T : Any> Field<T>.mstInField(block: MstField.() -> MS
*
* @author Iaroslav Postovalov
*/
public inline fun <reified T : Any> Field<T>.mstInExtendedField(block: MstExtendedField.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : ExtendedField<T>> A.mstInExtendedField(block: MstExtendedField.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return MstExpression(this, MstExtendedField.block())
}
@ -85,7 +86,7 @@ public inline fun <reified T : Any> Field<T>.mstInExtendedField(block: MstExtend
*
* @author Alexander Nozik
*/
public inline fun <reified T : Any, A : Space<T>> FunctionalExpressionSpace<T, A>.mstInSpace(block: MstSpace.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : Space<T>> FunctionalExpressionSpace<T, A>.mstInSpace(block: MstSpace.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return algebra.mstInSpace(block)
}
@ -95,7 +96,7 @@ public inline fun <reified T : Any, A : Space<T>> FunctionalExpressionSpace<T, A
*
* @author Alexander Nozik
*/
public inline fun <reified T : Any, A : Ring<T>> FunctionalExpressionRing<T, A>.mstInRing(block: MstRing.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : Ring<T>> FunctionalExpressionRing<T, A>.mstInRing(block: MstRing.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return algebra.mstInRing(block)
}
@ -105,7 +106,7 @@ public inline fun <reified T : Any, A : Ring<T>> FunctionalExpressionRing<T, A>.
*
* @author Alexander Nozik
*/
public inline fun <reified T : Any, A : Field<T>> FunctionalExpressionField<T, A>.mstInField(block: MstField.() -> MST): MstExpression<T> {
public inline fun <reified T : Any, A : Field<T>> FunctionalExpressionField<T, A>.mstInField(block: MstField.() -> MST): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return algebra.mstInField(block)
}
@ -116,8 +117,8 @@ public inline fun <reified T : Any, A : Field<T>> FunctionalExpressionField<T, A
* @author Iaroslav Postovalov
*/
public inline fun <reified T : Any, A : ExtendedField<T>> FunctionalExpressionExtendedField<T, A>.mstInExtendedField(
block: MstExtendedField.() -> MST
): MstExpression<T> {
block: MstExtendedField.() -> MST,
): MstExpression<T, A> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return algebra.mstInExtendedField(block)
}

3
kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/asm.kt
View File

@ -69,4 +69,5 @@ public inline fun <reified T : Any> Algebra<T>.expression(mst: MST): Expression<
*
* @author Alexander Nozik.
*/
public inline fun <reified T : Any> MstExpression<T>.compile(): Expression<T> = mst.compileWith(T::class.java, algebra)
public inline fun <reified T : Any> MstExpression<T, Algebra<T>>.compile(): Expression<T> =
mst.compileWith(T::class.java, algebra)

16
kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/AsmBuilder.kt
View File

@ -25,7 +25,7 @@ internal class AsmBuilder<T> internal constructor(
private val classOfT: Class<*>,
private val algebra: Algebra<T>,
private val className: String,
private val invokeLabel0Visitor: AsmBuilder<T>.() -> Unit
private val invokeLabel0Visitor: AsmBuilder<T>.() -> Unit,
) {
/**
* Internal classloader of [AsmBuilder] with alias to define class from byte array.
@ -379,22 +379,14 @@ internal class AsmBuilder<T> internal constructor(
* Loads a variable [name] from arguments [Map] parameter of [Expression.invoke]. The [defaultValue] may be
* provided.
*/
internal fun loadVariable(name: String, defaultValue: T? = null): Unit = invokeMethodVisitor.run {
internal fun loadVariable(name: String): Unit = invokeMethodVisitor.run {
load(invokeArgumentsVar, MAP_TYPE)
aconst(name)
if (defaultValue != null)
loadTConstant(defaultValue)
invokestatic(
MAP_INTRINSICS_TYPE.internalName,
"getOrFail",
Type.getMethodDescriptor(
OBJECT_TYPE,
MAP_TYPE,
OBJECT_TYPE,
*OBJECT_TYPE.wrapToArrayIf { defaultValue != null }),
Type.getMethodDescriptor(OBJECT_TYPE, MAP_TYPE, STRING_TYPE),
false
)
@ -429,7 +421,7 @@ internal class AsmBuilder<T> internal constructor(
method: String,
descriptor: String,
expectedArity: Int,
opcode: Int = INVOKEINTERFACE
opcode: Int = INVOKEINTERFACE,
) {
run loop@{
repeat(expectedArity) {

9
kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/mapIntrinsics.kt
View File

@ -2,11 +2,12 @@
package kscience.kmath.asm.internal
import kscience.kmath.expressions.StringSymbol
import kscience.kmath.expressions.Symbol
/**
* Gets value with given [key] or throws [IllegalStateException] whenever it is not present.
* Gets value with given [key] or throws [NoSuchElementException] whenever it is not present.
*
* @author Iaroslav Postovalov
*/
@JvmOverloads
internal fun <K, V> Map<K, V>.getOrFail(key: K, default: V? = null): V =
this[key] ?: default ?: error("Parameter not found: $key")
internal fun <V> Map<Symbol, V>.getOrFail(key: String): V = getValue(StringSymbol(key))

4
kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt
View File

@ -83,11 +83,10 @@ public object ArithmeticsEvaluator : Grammar<MST>() {
}
/**
* Tries to parse the string into [MST].
* Tries to parse the string into [MST]. Returns [ParseResult] representing expression or error.
*
* @receiver the string to parse.
* @return the [MST] node.
* @author Alexander Nozik
*/
public fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryParseToEnd(this)
@ -96,6 +95,5 @@ public fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryPar
*
* @receiver the string to parse.
* @return the [MST] node.
* @author Alexander Nozik
*/
public fun String.parseMath(): MST = ArithmeticsEvaluator.parseToEnd(this)

2
kmath-ast/src/jvmTest/kotlin/kscience/kmath/asm/TestAsmAlgebras.kt
View File

@ -1,6 +1,5 @@
package kscience.kmath.asm
import kscience.kmath.asm.compile
import kscience.kmath.ast.mstInField
import kscience.kmath.ast.mstInRing
import kscience.kmath.ast.mstInSpace
@ -11,6 +10,7 @@ import kotlin.test.Test
import kotlin.test.assertEquals
internal class TestAsmAlgebras {
@Test
fun space() {
val res1 = ByteRing.mstInSpace {

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

@ -6,7 +6,7 @@ description = "Commons math binding for kmath"
dependencies {
api(project(":kmath-core"))
api(project(":kmath-coroutines"))
api(project(":kmath-prob"))
// api(project(":kmath-functions"))
api(project(":kmath-stat"))
api(project(":kmath-functions"))
api("org.apache.commons:commons-math3:3.6.1")
}

121
kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt Normal file
View File

@ -0,0 +1,121 @@
package kscience.kmath.commons.expressions
import kscience.kmath.expressions.*
import kscience.kmath.operations.ExtendedField
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
/**
* A field over commons-math [DerivativeStructure].
*
* @property order The derivation order.
* @property bindings The map of bindings values. All bindings are considered free parameters
*/
public class DerivativeStructureField(
public val order: Int,
bindings: Map<Symbol, Double>,
) : ExtendedField<DerivativeStructure>, ExpressionAlgebra<Double, DerivativeStructure> {
public val numberOfVariables: Int = bindings.size
public override val zero: DerivativeStructure by lazy { DerivativeStructure(numberOfVariables, order) }
public override val one: DerivativeStructure by lazy { DerivativeStructure(numberOfVariables, order, 1.0) }
/**
* A class that implements both [DerivativeStructure] and a [Symbol]
*/
public inner class DerivativeStructureSymbol(
size: Int,
index: Int,
symbol: Symbol,
value: Double,
) : DerivativeStructure(size, order, index, value), Symbol {
override val identity: String = symbol.identity
override fun toString(): String = identity
override fun equals(other: Any?): Boolean = this.identity == (other as? Symbol)?.identity
override fun hashCode(): Int = identity.hashCode()
}
/**
* Identity-based symbol bindings map
*/
private val variables: Map<String, DerivativeStructureSymbol> = bindings.entries.mapIndexed { index, (key, value) ->
key.identity to DerivativeStructureSymbol(numberOfVariables, index, key, value)
}.toMap()
override fun const(value: Double): DerivativeStructure = DerivativeStructure(numberOfVariables, order, value)
public override fun bindOrNull(symbol: Symbol): DerivativeStructureSymbol? = variables[symbol.identity]
public fun bind(symbol: Symbol): DerivativeStructureSymbol = variables.getValue(symbol.identity)
override fun symbol(value: String): DerivativeStructureSymbol = bind(StringSymbol(value))
public fun DerivativeStructure.derivative(symbols: List<Symbol>): Double {
require(symbols.size <= order) { "The order of derivative ${symbols.size} exceeds computed order $order" }
val ordersCount = symbols.map { it.identity }.groupBy { it }.mapValues { it.value.size }
return getPartialDerivative(*variables.keys.map { ordersCount[it] ?: 0 }.toIntArray())
}
public fun DerivativeStructure.derivative(vararg symbols: Symbol): Double = derivative(symbols.toList())
public override fun add(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.add(b)
public 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())
}
public override fun multiply(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.multiply(b)
public override fun divide(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.divide(b)
public override fun sin(arg: DerivativeStructure): DerivativeStructure = arg.sin()
public override fun cos(arg: DerivativeStructure): DerivativeStructure = arg.cos()
public override fun tan(arg: DerivativeStructure): DerivativeStructure = arg.tan()
public override fun asin(arg: DerivativeStructure): DerivativeStructure = arg.asin()
public override fun acos(arg: DerivativeStructure): DerivativeStructure = arg.acos()
public override fun atan(arg: DerivativeStructure): DerivativeStructure = arg.atan()
public override fun sinh(arg: DerivativeStructure): DerivativeStructure = arg.sinh()
public override fun cosh(arg: DerivativeStructure): DerivativeStructure = arg.cosh()
public override fun tanh(arg: DerivativeStructure): DerivativeStructure = arg.tanh()
public override fun asinh(arg: DerivativeStructure): DerivativeStructure = arg.asinh()
public override fun acosh(arg: DerivativeStructure): DerivativeStructure = arg.acosh()
public override fun atanh(arg: DerivativeStructure): DerivativeStructure = arg.atanh()
public 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())
}
public fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure = arg.pow(pow)
public override fun exp(arg: DerivativeStructure): DerivativeStructure = arg.exp()
public override fun ln(arg: DerivativeStructure): DerivativeStructure = arg.log()
public override operator fun DerivativeStructure.plus(b: Number): DerivativeStructure = add(b.toDouble())
public override operator fun DerivativeStructure.minus(b: Number): DerivativeStructure = subtract(b.toDouble())
public override operator fun Number.plus(b: DerivativeStructure): DerivativeStructure = b + this
public override operator fun Number.minus(b: DerivativeStructure): DerivativeStructure = b - this
public companion object :
AutoDiffProcessor<Double, DerivativeStructure, DerivativeStructureField, Expression<Double>> {
public override fun process(function: DerivativeStructureField.() -> DerivativeStructure): DifferentiableExpression<Double, Expression<Double>> =
DerivativeStructureExpression(function)
}
}
/**
* A constructs that creates a derivative structure with required order on-demand
*/
public class DerivativeStructureExpression(
public val function: DerivativeStructureField.() -> DerivativeStructure,
) : DifferentiableExpression<Double, Expression<Double>> {
public override operator fun invoke(arguments: Map<Symbol, Double>): Double =
DerivativeStructureField(0, arguments).function().value
/**
* Get the derivative expression with given orders
*/
public override fun derivativeOrNull(symbols: List<Symbol>): Expression<Double> = Expression { arguments ->
with(DerivativeStructureField(symbols.size, arguments)) { function().derivative(symbols) }
}
}

131
kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DiffExpression.kt
View File

@ -1,131 +0,0 @@
package kscience.kmath.commons.expressions
import kscience.kmath.expressions.Expression
import kscience.kmath.expressions.ExpressionAlgebra
import kscience.kmath.operations.ExtendedField
import kscience.kmath.operations.Field
import kscience.kmath.operations.invoke
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
import kotlin.properties.ReadOnlyProperty
/**
* A field over commons-math [DerivativeStructure].
*
* @property order The derivation order.
* @property parameters The map of free parameters.
*/
public class DerivativeStructureField(
public val order: Int,
public val parameters: Map<String, Double>
) : ExtendedField<DerivativeStructure> {
public override val zero: DerivativeStructure by lazy { DerivativeStructure(parameters.size, order) }
public override val one: DerivativeStructure by lazy { DerivativeStructure(parameters.size, order, 1.0) }
private val variables: Map<String, DerivativeStructure> = parameters.mapValues { (key, value) ->
DerivativeStructure(parameters.size, order, parameters.keys.indexOf(key), value)
}
public val variable: ReadOnlyProperty<Any?, DerivativeStructure> = ReadOnlyProperty { _, property ->
variables[property.name] ?: error("A variable with name ${property.name} does not exist")
}
public fun variable(name: String, default: DerivativeStructure? = null): DerivativeStructure =
variables[name] ?: default ?: error("A variable with name $name does not exist")
public fun Number.const(): DerivativeStructure = DerivativeStructure(order, parameters.size, toDouble())
public fun DerivativeStructure.deriv(parName: String, order: Int = 1): Double {
return deriv(mapOf(parName to order))
}
public fun DerivativeStructure.deriv(orders: Map<String, Int>): Double {
return getPartialDerivative(*parameters.keys.map { orders[it] ?: 0 }.toIntArray())
}
public fun DerivativeStructure.deriv(vararg orders: Pair<String, Int>): Double = deriv(mapOf(*orders))
public override fun add(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.add(b)
public 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())
}
public override fun multiply(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.multiply(b)
public override fun divide(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.divide(b)
public override fun sin(arg: DerivativeStructure): DerivativeStructure = arg.sin()
public override fun cos(arg: DerivativeStructure): DerivativeStructure = arg.cos()
public override fun tan(arg: DerivativeStructure): DerivativeStructure = arg.tan()
public override fun asin(arg: DerivativeStructure): DerivativeStructure = arg.asin()
public override fun acos(arg: DerivativeStructure): DerivativeStructure = arg.acos()
public override fun atan(arg: DerivativeStructure): DerivativeStructure = arg.atan()
public override fun sinh(arg: DerivativeStructure): DerivativeStructure = arg.sinh()
public override fun cosh(arg: DerivativeStructure): DerivativeStructure = arg.cosh()
public override fun tanh(arg: DerivativeStructure): DerivativeStructure = arg.tanh()
public override fun asinh(arg: DerivativeStructure): DerivativeStructure = arg.asinh()
public override fun acosh(arg: DerivativeStructure): DerivativeStructure = arg.acosh()
public override fun atanh(arg: DerivativeStructure): DerivativeStructure = arg.atanh()
public 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())
}
public fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure = arg.pow(pow)
public override fun exp(arg: DerivativeStructure): DerivativeStructure = arg.exp()
public override fun ln(arg: DerivativeStructure): DerivativeStructure = arg.log()
public override operator fun DerivativeStructure.plus(b: Number): DerivativeStructure = add(b.toDouble())
public override operator fun DerivativeStructure.minus(b: Number): DerivativeStructure = subtract(b.toDouble())
public override operator fun Number.plus(b: DerivativeStructure): DerivativeStructure = b + this
public override operator fun Number.minus(b: DerivativeStructure): DerivativeStructure = b - this
}
/**
* A constructs that creates a derivative structure with required order on-demand
*/
public class DiffExpression(public val function: DerivativeStructureField.() -> DerivativeStructure) :
Expression<Double> {
public override operator fun invoke(arguments: Map<String, Double>): Double = DerivativeStructureField(
0,
arguments
).function().value
/**
* Get the derivative expression with given orders
* TODO make result [DiffExpression]
*/
public fun derivative(orders: Map<String, Int>): Expression<Double> = Expression { arguments ->
(DerivativeStructureField(orders.values.maxOrNull() ?: 0, arguments)) { function().deriv(orders) }
}
//TODO add gradient and maybe other vector operators
}
public fun DiffExpression.derivative(vararg orders: Pair<String, Int>): Expression<Double> = derivative(mapOf(*orders))
public fun DiffExpression.derivative(name: String): Expression<Double> = derivative(name to 1)
/**
* A context for [DiffExpression] (not to be confused with [DerivativeStructure])
*/
public object DiffExpressionAlgebra : ExpressionAlgebra<Double, DiffExpression>, Field<DiffExpression> {
public override val zero: DiffExpression = DiffExpression { 0.0.const() }
public override val one: DiffExpression = DiffExpression { 1.0.const() }
public override fun variable(name: String, default: Double?): DiffExpression =
DiffExpression { variable(name, default?.const()) }
public override fun const(value: Double): DiffExpression = DiffExpression { value.const() }
public override fun add(a: DiffExpression, b: DiffExpression): DiffExpression =
DiffExpression { a.function(this) + b.function(this) }
public override fun multiply(a: DiffExpression, k: Number): DiffExpression = DiffExpression { a.function(this) * k }
public override fun multiply(a: DiffExpression, b: DiffExpression): DiffExpression =
DiffExpression { a.function(this) * b.function(this) }
public override fun divide(a: DiffExpression, b: DiffExpression): DiffExpression =
DiffExpression { a.function(this) / b.function(this) }
}

7
kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
View File

@ -5,8 +5,7 @@ import kscience.kmath.structures.Matrix
import kscience.kmath.structures.NDStructure
import org.apache.commons.math3.linear.*
public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
FeaturedMatrix<Double> {
public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> {
public override val rowNum: Int get() = origin.rowDimension
public override val colNum: Int get() = origin.columnDimension
@ -55,7 +54,7 @@ public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
public fun RealVector.toPoint(): CMVector = CMVector(this)
public object CMMatrixContext : MatrixContext<Double> {
public object CMMatrixContext : MatrixContext<Double, CMMatrix> {
public 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))
@ -79,7 +78,7 @@ public object CMMatrixContext : MatrixContext<Double> {
public override fun multiply(a: Matrix<Double>, k: Number): CMMatrix =
CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
public override operator fun Matrix<Double>.times(value: Double): Matrix<Double> =
public override operator fun Matrix<Double>.times(value: Double): CMMatrix =
produce(rowNum, colNum) { i, j -> get(i, j) * value }
}

110
kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMOptimizationProblem.kt Normal file
View File

@ -0,0 +1,110 @@
package kscience.kmath.commons.optimization
import kscience.kmath.expressions.*
import kscience.kmath.stat.OptimizationFeature
import kscience.kmath.stat.OptimizationProblem
import kscience.kmath.stat.OptimizationProblemFactory
import kscience.kmath.stat.OptimizationResult
import org.apache.commons.math3.optim.*
import org.apache.commons.math3.optim.nonlinear.scalar.GoalType
import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction
import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunctionGradient
import org.apache.commons.math3.optim.nonlinear.scalar.gradient.NonLinearConjugateGradientOptimizer
import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.AbstractSimplex
import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex
import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer
import kotlin.reflect.KClass
public operator fun PointValuePair.component1(): DoubleArray = point
public operator fun PointValuePair.component2(): Double = value
public class CMOptimizationProblem(override val symbols: List<Symbol>, ) :
OptimizationProblem<Double>, SymbolIndexer, OptimizationFeature {
private val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
private var optimizatorBuilder: (() -> MultivariateOptimizer)? = null
public var convergenceChecker: ConvergenceChecker<PointValuePair> = SimpleValueChecker(DEFAULT_RELATIVE_TOLERANCE,
DEFAULT_ABSOLUTE_TOLERANCE, DEFAULT_MAX_ITER)
public fun addOptimizationData(data: OptimizationData) {
optimizationData[data::class] = data
}
init {
addOptimizationData(MaxEval.unlimited())
}
public fun exportOptimizationData(): List<OptimizationData> = optimizationData.values.toList()
public override fun initialGuess(map: Map<Symbol, Double>): Unit {
addOptimizationData(InitialGuess(map.toDoubleArray()))
}
public override fun expression(expression: Expression<Double>): Unit {
val objectiveFunction = ObjectiveFunction {
val args = it.toMap()
expression(args)
}
addOptimizationData(objectiveFunction)
}
public override fun diffExpression(expression: DifferentiableExpression<Double, Expression<Double>>) {
expression(expression)
val gradientFunction = ObjectiveFunctionGradient {
val args = it.toMap()
DoubleArray(symbols.size) { index ->
expression.derivative(symbols[index])(args)
}
}
addOptimizationData(gradientFunction)
if (optimizatorBuilder == null) {
optimizatorBuilder = {
NonLinearConjugateGradientOptimizer(
NonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES,
convergenceChecker
)
}
}
}
public fun simplex(simplex: AbstractSimplex) {
addOptimizationData(simplex)
//Set optimization builder to simplex if it is not present
if (optimizatorBuilder == null) {
optimizatorBuilder = { SimplexOptimizer(convergenceChecker) }
}
}
public fun simplexSteps(steps: Map<Symbol, Double>) {
simplex(NelderMeadSimplex(steps.toDoubleArray()))
}
public fun goal(goalType: GoalType) {
addOptimizationData(goalType)
}
public fun optimizer(block: () -> MultivariateOptimizer) {
optimizatorBuilder = block
}
override fun update(result: OptimizationResult<Double>) {
initialGuess(result.point)
}
override fun optimize(): OptimizationResult<Double> {
val optimizer = optimizatorBuilder?.invoke() ?: error("Optimizer not defined")
val (point, value) = optimizer.optimize(*optimizationData.values.toTypedArray())
return OptimizationResult(point.toMap(), value, setOf(this))
}
public companion object : OptimizationProblemFactory<Double, CMOptimizationProblem> {
public const val DEFAULT_RELATIVE_TOLERANCE: Double = 1e-4
public const val DEFAULT_ABSOLUTE_TOLERANCE: Double = 1e-4
public const val DEFAULT_MAX_ITER: Int = 1000
override fun build(symbols: List<Symbol>): CMOptimizationProblem = CMOptimizationProblem(symbols)
}
}
public fun CMOptimizationProblem.initialGuess(vararg pairs: Pair<Symbol, Double>): Unit = initialGuess(pairs.toMap())
public fun CMOptimizationProblem.simplexSteps(vararg pairs: Pair<Symbol, Double>): Unit = simplexSteps(pairs.toMap())

67
kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/cmFit.kt Normal file
View File

@ -0,0 +1,67 @@
package kscience.kmath.commons.optimization
import kscience.kmath.commons.expressions.DerivativeStructureField
import kscience.kmath.expressions.DifferentiableExpression
import kscience.kmath.expressions.Expression
import kscience.kmath.expressions.Symbol
import kscience.kmath.stat.Fitting
import kscience.kmath.stat.OptimizationResult
import kscience.kmath.stat.optimizeWith
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.asBuffer
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
import org.apache.commons.math3.optim.nonlinear.scalar.GoalType
/**
* Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
*/
public fun Fitting.chiSquared(
x: Buffer<Double>,
y: Buffer<Double>,
yErr: Buffer<Double>,
model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
): DifferentiableExpression<Double, Expression<Double>> = chiSquared(DerivativeStructureField, x, y, yErr, model)
/**
* Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
*/
public fun Fitting.chiSquared(
x: Iterable<Double>,
y: Iterable<Double>,
yErr: Iterable<Double>,
model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
): DifferentiableExpression<Double, Expression<Double>> = chiSquared(
DerivativeStructureField,
x.toList().asBuffer(),
y.toList().asBuffer(),
yErr.toList().asBuffer(),
model
)
/**
* Optimize expression without derivatives
*/
public fun Expression<Double>.optimize(
vararg symbols: Symbol,
configuration: CMOptimizationProblem.() -> Unit,
): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
/**
* Optimize differentiable expression
*/
public fun DifferentiableExpression<Double, Expression<Double>>.optimize(
vararg symbols: Symbol,
configuration: CMOptimizationProblem.() -> Unit,
): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
public fun DifferentiableExpression<Double, Expression<Double>>.minimize(
vararg startPoint: Pair<Symbol, Double>,
configuration: CMOptimizationProblem.() -> Unit = {},
): OptimizationResult<Double> {
require(startPoint.isNotEmpty()) { "Must provide a list of symbols for optimization" }
val problem = CMOptimizationProblem(startPoint.map { it.first }).apply(configuration)
problem.diffExpression(this)
problem.initialGuess(startPoint.toMap())
problem.goal(GoalType.MINIMIZE)
return problem.optimize()
}

7
kmath-commons/src/main/kotlin/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt
View File

@ -1,9 +1,10 @@
package kscience.kmath.commons.random
import kscience.kmath.prob.RandomGenerator
import kscience.kmath.stat.RandomGenerator
public class CMRandomGeneratorWrapper(public val factory: (IntArray) -> RandomGenerator) :
org.apache.commons.math3.random.RandomGenerator {
public class CMRandomGeneratorWrapper(
public val factory: (IntArray) -> RandomGenerator,
) : org.apache.commons.math3.random.RandomGenerator {
private var generator: RandomGenerator = factory(intArrayOf())
public override fun nextBoolean(): Boolean = generator.nextBoolean()

40
kmath-commons/src/test/kotlin/kscience/kmath/commons/expressions/AutoDiffTest.kt
View File

@ -1,40 +0,0 @@
package kscience.kmath.commons.expressions
import kscience.kmath.expressions.invoke
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
import kotlin.test.Test
import kotlin.test.assertEquals
internal inline fun <R> diff(
order: Int,
vararg parameters: Pair<String, Double>,
block: DerivativeStructureField.() -> R
): R {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return DerivativeStructureField(order, mapOf(*parameters)).run(block)
}
internal 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))
}
}

50
kmath-commons/src/test/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt Normal file
View File

@ -0,0 +1,50 @@
package kscience.kmath.commons.expressions
import kscience.kmath.expressions.*
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertFails
internal inline fun diff(
order: Int,
vararg parameters: Pair<Symbol, Double>,
block: DerivativeStructureField.() -> Unit,
): Unit {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
DerivativeStructureField(order, mapOf(*parameters)).run(block)
}
internal class AutoDiffTest {
private val x by symbol
private val y by symbol
@Test
fun derivativeStructureFieldTest() {
diff(2, x to 1.0, y to 1.0) {
val x = bind(x)//by binding()
val y = symbol("y")
val z = x * (-sin(x * y) + y) + 2.0
println(z.derivative(x))
println(z.derivative(y,x))
assertEquals(z.derivative(x, y), z.derivative(y, x))
//check that improper order cause failure
assertFails { z.derivative(x,x,y) }
}
}
@Test
fun autoDifTest() {
val f = DerivativeStructureExpression {
val x by binding()
val y by binding()
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))
assertEquals(2.0, f.derivative(x, x)(x to 1.234, y to -2.0))
assertEquals(2.0, f.derivative(x, y)(x to 1.0, y to 2.0))
}
}

68
kmath-commons/src/test/kotlin/kscience/kmath/commons/optimization/OptimizeTest.kt Normal file
View File

@ -0,0 +1,68 @@
package kscience.kmath.commons.optimization
import kscience.kmath.commons.expressions.DerivativeStructureExpression
import kscience.kmath.expressions.symbol
import kscience.kmath.stat.Distribution
import kscience.kmath.stat.Fitting
import kscience.kmath.stat.RandomGenerator
import kscience.kmath.stat.normal
import org.junit.jupiter.api.Test
import kotlin.math.pow
internal class OptimizeTest {
val x by symbol
val y by symbol
val normal = DerivativeStructureExpression {
exp(-bind(x).pow(2) / 2) + exp(-bind(y).pow(2) / 2)
}
@Test
fun testGradientOptimization() {
val result = normal.optimize(x, y) {
initialGuess(x to 1.0, y to 1.0)
//no need to select optimizer. Gradient optimizer is used by default because gradients are provided by function
}
println(result.point)
println(result.value)
}
@Test
fun testSimplexOptimization() {
val result = normal.optimize(x, y) {
initialGuess(x to 1.0, y to 1.0)
simplexSteps(x to 2.0, y to 0.5)
//this sets simplex optimizer
}
println(result.point)
println(result.value)
}
@Test
fun testCmFit() {
val a by symbol
val b by symbol
val c by symbol
val sigma = 1.0
val generator = Distribution.normal(0.0, sigma)
val chain = generator.sample(RandomGenerator.default(112667))
val x = (1..100).map(Int::toDouble)
val y = x.map {
it.pow(2) + it + 1 + chain.nextDouble()
}
val yErr = List(x.size) { sigma }
val chi2 = Fitting.chiSquared(x, y, yErr) { x1 ->
val cWithDefault = bindOrNull(c) ?: one
bind(a) * x1.pow(2) + bind(b) * x1 + cWithDefault
}
val result = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
println(result)
println("Chi2/dof = ${result.value / (x.size - 3)}")
}
}

11
kmath-core/README.md
View File

@ -7,12 +7,12 @@ The core features of KMath:
- [buffers](src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
- [expressions](src/commonMain/kotlin/kscience/kmath/expressions) : Functional Expressions
- [domains](src/commonMain/kotlin/kscience/kmath/domains) : Domains
- [autodif](src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt) : Automatic differentiation
- [autodif](src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
> #### Artifact:
>
> This module artifact: `kscience.kmath:kmath-core:0.2.0-dev-1`.
> This module artifact: `kscience.kmath:kmath-core:0.2.0-dev-4`.
>
> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-core/_latestVersion)
>
@ -22,25 +22,28 @@ The core features of KMath:
>
> ```gradle
> repositories {
> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" }
> maven { url 'https://dl.bintray.com/mipt-npm/kscience' }
> maven { url 'https://dl.bintray.com/mipt-npm/dev' }
> maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
> }
>
> dependencies {
> implementation 'kscience.kmath:kmath-core:0.2.0-dev-1'
> implementation 'kscience.kmath:kmath-core:0.2.0-dev-4'
> }
> ```
> **Gradle Kotlin DSL:**
>
> ```kotlin
> repositories {
> maven("https://dl.bintray.com/kotlin/kotlin-eap")
> maven("https://dl.bintray.com/mipt-npm/kscience")
> maven("https://dl.bintray.com/mipt-npm/dev")
> maven("https://dl.bintray.com/hotkeytlt/maven")
> }
>
> dependencies {
> implementation("kscience.kmath:kmath-core:0.2.0-dev-1")
> implementation("kscience.kmath:kmath-core:0.2.0-dev-4")
> }
> ```

8
kmath-core/build.gradle.kts
View File

@ -13,34 +13,40 @@ readme {
description = "Core classes, algebra definitions, basic linear algebra"
maturity = ru.mipt.npm.gradle.Maturity.DEVELOPMENT
propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
feature(
id = "algebras",
description = "Algebraic structures: contexts and elements",
ref = "src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt"
)
feature(
id = "nd",
description = "Many-dimensional structures",
ref = "src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt"
)
feature(
id = "buffers",
description = "One-dimensional structure",
ref = "src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt"
)
feature(
id = "expressions",
description = "Functional Expressions",
ref = "src/commonMain/kotlin/kscience/kmath/expressions"
)
feature(
id = "domains",
description = "Domains",
ref = "src/commonMain/kotlin/kscience/kmath/domains"
)
feature(
id = "autodif",
description = "Automatic differentiation",
ref = "src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt"
ref = "src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt"
)
}

48
kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/DifferentiableExpression.kt Normal file
View File

@ -0,0 +1,48 @@
package kscience.kmath.expressions
/**
* Represents expression which structure can be differentiated.
*
* @param T the type this expression takes as argument and returns.
* @param R the type of expression this expression can be differentiated to.
*/
public interface DifferentiableExpression<T, out R : Expression<T>> : Expression<T> {
/**
* Differentiates this expression by ordered collection of [symbols].
*
* @param symbols the symbols.
* @return the derivative or `null`.
*/
public fun derivativeOrNull(symbols: List<Symbol>): R?
}
public fun <T, R : Expression<T>> DifferentiableExpression<T, R>.derivative(symbols: List<Symbol>): R =
derivativeOrNull(symbols) ?: error("Derivative by symbols $symbols not provided")
public fun <T, R : Expression<T>> DifferentiableExpression<T, R>.derivative(vararg symbols: Symbol): R =
derivative(symbols.toList())
public fun <T, R : Expression<T>> DifferentiableExpression<T, R>.derivative(name: String): R =
derivative(StringSymbol(name))
/**
* A [DifferentiableExpression] that defines only first derivatives
*/
public abstract class FirstDerivativeExpression<T, R : Expression<T>> : DifferentiableExpression<T,R> {
/**
* Returns first derivative of this expression by given [symbol].
*/
public abstract fun derivativeOrNull(symbol: Symbol): R?
public final override fun derivativeOrNull(symbols: List<Symbol>): R? {
val dSymbol = symbols.firstOrNull() ?: return null
return derivativeOrNull(dSymbol)
}
}
/**
* A factory that converts an expression in autodiff variables to a [DifferentiableExpression]
*/
public fun interface AutoDiffProcessor<T : Any, I : Any, A : ExpressionAlgebra<T, I>, out R : Expression<T>> {
public fun process(function: A.() -> I): DifferentiableExpression<T, R>
}

94
kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/Expression.kt
View File

@ -1,9 +1,38 @@
package kscience.kmath.expressions
import kscience.kmath.operations.Algebra
import kotlin.jvm.JvmName
import kotlin.properties.ReadOnlyProperty
import kotlin.reflect.KProperty
/**
* An elementary function that could be invoked on a map of arguments
* A marker interface for a symbol. A symbol mus have an identity
*/
public interface Symbol {
/**
* Identity object for the symbol. Two symbols with the same identity are considered to be the same symbol.
*/
public val identity: String
public companion object : ReadOnlyProperty<Any?, Symbol> {
//TODO deprecate and replace by top level function after fix of https://youtrack.jetbrains.com/issue/KT-40121
override fun getValue(thisRef: Any?, property: KProperty<*>): Symbol {
return StringSymbol(property.name)
}
}
}
/**
* A [Symbol] with a [String] identity
*/
public inline class StringSymbol(override val identity: String) : Symbol {
override fun toString(): String = identity
}
/**
* An elementary function that could be invoked on a map of arguments.
*
* @param T the type this expression takes as argument and returns.
*/
public fun interface Expression<T> {
/**
@ -12,30 +41,75 @@ public fun interface Expression<T> {
* @param arguments the map of arguments.
* @return the value.
*/
public operator fun invoke(arguments: Map<String, T>): T
public companion object
public operator fun invoke(arguments: Map<Symbol, T>): T
}
/**
* Calls this expression without providing any arguments.
*
* @return a value.
*/
public operator fun <T> Expression<T>.invoke(): T = invoke(emptyMap())
/**
* Calls this expression from arguments.
*
* @param pairs the pair of arguments' names to values.
* @return the value.
* @param pairs the pairs of arguments to values.
* @return a value.
*/
public operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs))
@JvmName("callBySymbol")
public operator fun <T> Expression<T>.invoke(vararg pairs: Pair<Symbol, T>): T = invoke(mapOf(*pairs))
/**
* Calls this expression from arguments.
*
* @param pairs the pairs of arguments' names to values.
* @return a value.
*/
@JvmName("callByString")
public operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T =
invoke(mapOf(*pairs).mapKeys { StringSymbol(it.key) })
/**
* A context for expression construction
*
* @param T type of the constants for the expression
* @param E type of the actual expression state
*/
public interface ExpressionAlgebra<T, E> : Algebra<E> {
public interface ExpressionAlgebra<in T, E> : Algebra<E> {
/**
* Introduce a variable into expression context
* Bind a given [Symbol] to this context variable and produce context-specific object. Return null if symbol could not be bound in current context.
*/
public fun variable(name: String, default: T? = null): E
public fun bindOrNull(symbol: Symbol): E?
/**
* Bind a string to a context using [StringSymbol]
*/
override fun symbol(value: String): E = bind(StringSymbol(value))
/**
* A constant expression which does not depend on arguments
*/
public fun const(value: T): E
}
/**
* Bind a given [Symbol] to this context variable and produce context-specific object.
*/
public fun <T, E> ExpressionAlgebra<T, E>.bind(symbol: Symbol): E =
bindOrNull(symbol) ?: error("Symbol $symbol could not be bound to $this")
/**
* A delegate to create a symbol with a string identity in this scope
*/
public val symbol: ReadOnlyProperty<Any?, Symbol> get() = Symbol
//TODO does not work directly on native due to https://youtrack.jetbrains.com/issue/KT-40121
/**
* Bind a symbol by name inside the [ExpressionAlgebra]
*/
public fun <T, E> ExpressionAlgebra<T, E>.binding(): ReadOnlyProperty<Any?, E> = ReadOnlyProperty { _, property ->
bind(StringSymbol(property.name)) ?: error("A variable with name ${property.name} does not exist")
}

67
kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt
View File

@ -2,67 +2,43 @@ package kscience.kmath.expressions
import kscience.kmath.operations.*
internal class FunctionalUnaryOperation<T>(val context: Algebra<T>, val name: String, private val expr: Expression<T>) :
Expression<T> {
override operator fun invoke(arguments: Map<String, T>): T =
context.unaryOperation(name, expr.invoke(arguments))
}
internal class FunctionalBinaryOperation<T>(
val context: Algebra<T>,
val name: String,
val first: Expression<T>,
val second: Expression<T>
) : Expression<T> {
override operator fun invoke(arguments: Map<String, T>): T =
context.binaryOperation(name, first.invoke(arguments), second.invoke(arguments))
}
internal class FunctionalVariableExpression<T>(val name: String, val default: T? = null) : Expression<T> {
override operator fun invoke(arguments: Map<String, T>): T =
arguments[name] ?: default ?: error("Parameter not found: $name")
}
internal class FunctionalConstantExpression<T>(val value: T) : Expression<T> {
override operator fun invoke(arguments: Map<String, T>): T = value
}
internal class FunctionalConstProductExpression<T>(
val context: Space<T>,
private val expr: Expression<T>,
val const: Number
) : Expression<T> {
override operator fun invoke(arguments: Map<String, T>): T = context.multiply(expr.invoke(arguments), const)
}
/**
* A context class for [Expression] construction.
*
* @param algebra The algebra to provide for Expressions built.
*/
public abstract class FunctionalExpressionAlgebra<T, A : Algebra<T>>(public val algebra: A) :
ExpressionAlgebra<T, Expression<T>> {
public abstract class FunctionalExpressionAlgebra<T, A : Algebra<T>>(
public val algebra: A,
) : ExpressionAlgebra<T, Expression<T>> {
/**
* Builds an Expression of constant expression which does not depend on arguments.
*/
public override fun const(value: T): Expression<T> = FunctionalConstantExpression(value)
public override fun const(value: T): Expression<T> = Expression { value }
/**
* Builds an Expression to access a variable.
*/
public override fun variable(name: String, default: T?): Expression<T> = FunctionalVariableExpression(name, default)
public override fun bindOrNull(symbol: Symbol): Expression<T>? = Expression { arguments ->
arguments[symbol] ?: error("Argument not found: $symbol")
}
/**
* Builds an Expression of dynamic call of binary operation [operation] on [left] and [right].
*/
public override fun binaryOperation(operation: String, left: Expression<T>, right: Expression<T>): Expression<T> =
FunctionalBinaryOperation(algebra, operation, left, right)
public override fun binaryOperation(
operation: String,
left: Expression<T>,
right: Expression<T>,
): Expression<T> = Expression { arguments ->
algebra.binaryOperation(operation, left.invoke(arguments), right.invoke(arguments))
}
/**
* Builds an Expression of dynamic call of unary operation with name [operation] on [arg].
*/
public override fun unaryOperation(operation: String, arg: Expression<T>): Expression<T> =
FunctionalUnaryOperation(algebra, operation, arg)
public override fun unaryOperation(operation: String, arg: Expression<T>): Expression<T> = Expression { arguments ->
algebra.unaryOperation(operation, arg.invoke(arguments))
}
}
/**
@ -81,8 +57,9 @@ public open class FunctionalExpressionSpace<T, A : Space<T>>(algebra: A) :
/**
* Builds an Expression of multiplication of expression by number.
*/
public override fun multiply(a: Expression<T>, k: Number): Expression<T> =
FunctionalConstProductExpression(algebra, a, k)
public override fun multiply(a: Expression<T>, k: Number): Expression<T> = Expression { arguments ->
algebra.multiply(a.invoke(arguments), k)
}
public operator fun Expression<T>.plus(arg: T): Expression<T> = this + const(arg)
public operator fun Expression<T>.minus(arg: T): Expression<T> = this - const(arg)
@ -118,8 +95,8 @@ public open class FunctionalExpressionRing<T, A>(algebra: A) : FunctionalExpress
}
public open class FunctionalExpressionField<T, A>(algebra: A) :
FunctionalExpressionRing<T, A>(algebra),
Field<Expression<T>> where A : Field<T>, A : NumericAlgebra<T> {
FunctionalExpressionRing<T, A>(algebra), Field<Expression<T>>
where A : Field<T>, A : NumericAlgebra<T> {
/**
* Builds an Expression of division an expression by another one.
*/

393
kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt Normal file
View File

@ -0,0 +1,393 @@
package kscience.kmath.expressions
import kscience.kmath.linear.Point
import kscience.kmath.operations.*
import kscience.kmath.structures.asBuffer
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
/*
* Implementation of backward-mode automatic differentiation.
* Initial gist by Roman Elizarov: https://gist.github.com/elizarov/1ad3a8583e88cb6ea7a0ad09bb591d3d
*/
public open class AutoDiffValue<out T>(public val value: T)
/**
* Represents result of [simpleAutoDiff] call.
*
* @param T the non-nullable type of value.
* @param value the value of result.
* @property simpleAutoDiff The mapping of differentiated variables to their derivatives.
* @property context The field over [T].
*/
public class DerivationResult<T : Any>(
public val value: T,
private val derivativeValues: Map<String, T>,
public val context: Field<T>,
) {
/**
* Returns derivative of [variable] or returns [Ring.zero] in [context].
*/
public fun derivative(variable: Symbol): T = derivativeValues[variable.identity] ?: context.zero
/**
* Computes the divergence.
*/
public fun div(): T = context { sum(derivativeValues.values) }
}
/**
* Computes the gradient for variables in given order.
*/
public fun <T : Any> DerivationResult<T>.grad(vararg variables: Symbol): Point<T> {
check(variables.isNotEmpty()) { "Variable order is not provided for gradient construction" }
return variables.map(::derivative).asBuffer()
}
/**
* Runs differentiation and establishes [SimpleAutoDiffField] context inside the block of code.
*
* The partial derivatives are placed in argument `d` variable
*
* Example:
* ```
* val x by symbol // define variable(s) and their values
* val y = RealField.withAutoDiff() { sqr(x) + 5 * x + 3 } // write formulate in deriv context
* assertEquals(17.0, y.x) // the value of result (y)
* assertEquals(9.0, x.d) // dy/dx
* ```
*
* @param body the action in [SimpleAutoDiffField] context returning [AutoDiffVariable] to differentiate with respect to.
* @return the result of differentiation.
*/
public fun <T : Any, F : Field<T>> F.simpleAutoDiff(
bindings: Map<Symbol, T>,
body: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>,
): DerivationResult<T> {
contract { callsInPlace(body, InvocationKind.EXACTLY_ONCE) }
return SimpleAutoDiffField(this, bindings).differentiate(body)
}
public fun <T : Any, F : Field<T>> F.simpleAutoDiff(
vararg bindings: Pair<Symbol, T>,
body: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>,
): DerivationResult<T> = simpleAutoDiff(bindings.toMap(), body)
/**
* Represents field in context of which functions can be derived.
*/
public open class SimpleAutoDiffField<T : Any, F : Field<T>>(
public val context: F,
bindings: Map<Symbol, T>,
) : Field<AutoDiffValue<T>>, ExpressionAlgebra<T, AutoDiffValue<T>> {
public override val zero: AutoDiffValue<T>
get() = const(context.zero)
public override val one: AutoDiffValue<T>
get() = const(context.one)
// this stack contains pairs of blocks and values to apply them to
private var stack: Array<Any?> = arrayOfNulls<Any?>(8)
private var sp: Int = 0
private val derivatives: MutableMap<AutoDiffValue<T>, T> = hashMapOf()
private val bindings: Map<String, AutoDiffVariableWithDerivative<T>> = bindings.entries.associate {
it.key.identity to AutoDiffVariableWithDerivative(it.key.identity, it.value, context.zero)
}
/**
* Differentiable variable with value and derivative of differentiation ([simpleAutoDiff]) result
* with respect to this variable.
*
* @param T the non-nullable type of value.
* @property value The value of this variable.
*/
private class AutoDiffVariableWithDerivative<T : Any>(
override val identity: String,
value: T,
var d: T,
) : AutoDiffValue<T>(value), Symbol {
override fun toString(): String = identity
override fun equals(other: Any?): Boolean = this.identity == (other as? Symbol)?.identity
override fun hashCode(): Int = identity.hashCode()
}
public override fun bindOrNull(symbol: Symbol): AutoDiffValue<T>? = bindings[symbol.identity]
private fun getDerivative(variable: AutoDiffValue<T>): T =
(variable as? AutoDiffVariableWithDerivative)?.d ?: derivatives[variable] ?: context.zero
private fun setDerivative(variable: AutoDiffValue<T>, value: T) {
if (variable is AutoDiffVariableWithDerivative) variable.d = value else derivatives[variable] = value
}
@Suppress("UNCHECKED_CAST")
private fun runBackwardPass() {
while (sp > 0) {
val value = stack[--sp]
val block = stack[--sp] as F.(Any?) -> Unit
context.block(value)
}
}
override fun const(value: T): AutoDiffValue<T> = AutoDiffValue(value)
/**
* A variable accessing inner state of derivatives.
* Use this value in inner builders to avoid creating additional derivative bindings.
*/
public var AutoDiffValue<T>.d: T
get() = getDerivative(this)
set(value) = setDerivative(this, value)
public inline fun const(block: F.() -> T): AutoDiffValue<T> = const(context.block())
/**
* Performs update of derivative after the rest of the formula in the back-pass.
*
* For example, implementation of `sin` function is:
*
* ```
* fun AD.sin(x: Variable): Variable = derive(Variable(sin(x.x)) { z -> // call derive with function result
* x.d += z.d * cos(x.x) // update derivative using chain rule and derivative of the function
* }
* ```
*/
@Suppress("UNCHECKED_CAST")
public fun <R> derive(value: R, block: F.(R) -> Unit): R {
// save block to stack for backward pass
if (sp >= stack.size) stack = stack.copyOf(stack.size * 2)
stack[sp++] = block
stack[sp++] = value
return value
}
internal fun differentiate(function: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>): DerivationResult<T> {
val result = function()
result.d = context.one // computing derivative w.r.t result
runBackwardPass()
return DerivationResult(result.value, bindings.mapValues { it.value.d }, context)
}
// Overloads for Double constants
public override operator fun Number.plus(b: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { this@plus.toDouble() * one + b.value }) { z ->
b.d += z.d
}
public override operator fun AutoDiffValue<T>.plus(b: Number): AutoDiffValue<T> = b.plus(this)
public override operator fun Number.minus(b: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { this@minus.toDouble() * one - b.value }) { z -> b.d -= z.d }
public override operator fun AutoDiffValue<T>.minus(b: Number): AutoDiffValue<T> =
derive(const { this@minus.value - one * b.toDouble() }) { z -> this@minus.d += z.d }
// Basic math (+, -, *, /)
public override fun add(a: AutoDiffValue<T>, b: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { a.value + b.value }) { z ->
a.d += z.d
b.d += z.d
}
public override fun multiply(a: AutoDiffValue<T>, b: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { a.value * b.value }) { z ->
a.d += z.d * b.value
b.d += z.d * a.value
}
public override fun divide(a: AutoDiffValue<T>, b: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { a.value / b.value }) { z ->
a.d += z.d / b.value
b.d -= z.d * a.value / (b.value * b.value)
}
public override fun multiply(a: AutoDiffValue<T>, k: Number): AutoDiffValue<T> =
derive(const { k.toDouble() * a.value }) { z ->
a.d += z.d * k.toDouble()
}
}
/**
* A constructs that creates a derivative structure with required order on-demand
*/
public class SimpleAutoDiffExpression<T : Any, F : Field<T>>(
public val field: F,
public val function: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>,
) : FirstDerivativeExpression<T, Expression<T>>() {
public override operator fun invoke(arguments: Map<Symbol, T>): T {
//val bindings = arguments.entries.map { it.key.bind(it.value) }
return SimpleAutoDiffField(field, arguments).function().value
}
public override fun derivativeOrNull(symbol: Symbol): Expression<T> = Expression { arguments ->
//val bindings = arguments.entries.map { it.key.bind(it.value) }
val derivationResult = SimpleAutoDiffField(field, arguments).differentiate(function)
derivationResult.derivative(symbol)
}
}
/**
* Generate [AutoDiffProcessor] for [SimpleAutoDiffExpression]
*/
public fun <T : Any, F : Field<T>> simpleAutoDiff(field: F): AutoDiffProcessor<T, AutoDiffValue<T>, SimpleAutoDiffField<T, F>, Expression<T>> =
AutoDiffProcessor { function ->
SimpleAutoDiffExpression(field, function)
}
// Extensions for differentiation of various basic mathematical functions
// x ^ 2
public fun <T : Any, F : Field<T>> SimpleAutoDiffField<T, F>.sqr(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { x.value * x.value }) { z -> x.d += z.d * 2 * x.value }
// x ^ 1/2
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sqrt(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { sqrt(x.value) }) { z -> x.d += z.d * 0.5 / z.value }
// x ^ y (const)
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
x: AutoDiffValue<T>,
y: Double,
): AutoDiffValue<T> =
derive(const { power(x.value, y) }) { z -> x.d += z.d * y * power(x.value, y - 1) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
x: AutoDiffValue<T>,
y: Int,
): AutoDiffValue<T> = pow(x, y.toDouble())
// exp(x)
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.exp(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { exp(x.value) }) { z -> x.d += z.d * z.value }
// ln(x)
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.ln(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { ln(x.value) }) { z -> x.d += z.d / x.value }
// x ^ y (any)
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
x: AutoDiffValue<T>,
y: AutoDiffValue<T>,
): AutoDiffValue<T> =
exp(y * ln(x))
// sin(x)
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sin(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { sin(x.value) }) { z -> x.d += z.d * cos(x.value) }
// cos(x)
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.cos(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { cos(x.value) }) { z -> x.d -= z.d * sin(x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.tan(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { tan(x.value) }) { z ->
val c = cos(x.value)
x.d += z.d / (c * c)
}
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.asin(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { asin(x.value) }) { z -> x.d += z.d / sqrt(one - x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.acos(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { acos(x.value) }) { z -> x.d -= z.d / sqrt(one - x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.atan(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { atan(x.value) }) { z -> x.d += z.d / (one + x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sinh(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { sinh(x.value) }) { z -> x.d += z.d * cosh(x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.cosh(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { cosh(x.value) }) { z -> x.d += z.d * sinh(x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.tanh(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { tanh(x.value) }) { z ->
val c = cosh(x.value)
x.d += z.d / (c * c)
}
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.asinh(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { asinh(x.value) }) { z -> x.d += z.d / sqrt(one + x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.acosh(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { acosh(x.value) }) { z -> x.d += z.d / (sqrt((x.value - one) * (x.value + one))) }
public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.atanh(x: AutoDiffValue<T>): AutoDiffValue<T> =
derive(const { atanh(x.value) }) { z -> x.d += z.d / (one - x.value * x.value) }
public class SimpleAutoDiffExtendedField<T : Any, F : ExtendedField<T>>(
context: F,
bindings: Map<Symbol, T>,
) : ExtendedField<AutoDiffValue<T>>, SimpleAutoDiffField<T, F>(context, bindings) {
// x ^ 2
public fun sqr(x: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).sqr(x)
// x ^ 1/2
public override fun sqrt(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).sqrt(arg)
// x ^ y (const)
public override fun power(arg: AutoDiffValue<T>, pow: Number): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).pow(arg, pow.toDouble())
// exp(x)
public override fun exp(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).exp(arg)
// ln(x)
public override fun ln(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).ln(arg)
// x ^ y (any)
public fun pow(
x: AutoDiffValue<T>,
y: AutoDiffValue<T>,
): AutoDiffValue<T> = exp(y * ln(x))
// sin(x)
public override fun sin(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).sin(arg)
// cos(x)
public override fun cos(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).cos(arg)
public override fun tan(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).tan(arg)
public override fun asin(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).asin(arg)
public override fun acos(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).acos(arg)
public override fun atan(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).atan(arg)
public override fun sinh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).sinh(arg)
public override fun cosh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).cosh(arg)
public override fun tanh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).tanh(arg)
public override fun asinh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).asinh(arg)
public override fun acosh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).acosh(arg)
public override fun atanh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
(this as SimpleAutoDiffField<T, F>).atanh(arg)
}

61
kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SymbolIndexer.kt Normal file
View File

@ -0,0 +1,61 @@
package kscience.kmath.expressions
import kscience.kmath.linear.Point
import kscience.kmath.structures.BufferFactory
import kscience.kmath.structures.Structure2D
/**
* An environment to easy transform indexed variables to symbols and back.
* TODO requires multi-receivers to be beutiful
*/
public interface SymbolIndexer {
public val symbols: List<Symbol>
public fun indexOf(symbol: Symbol): Int = symbols.indexOf(symbol)
public operator fun <T> List<T>.get(symbol: Symbol): T {
require(size == symbols.size) { "The input list size for indexer should be ${symbols.size} but $size found" }
return get(this@SymbolIndexer.indexOf(symbol))
}
public operator fun <T> Array<T>.get(symbol: Symbol): T {
require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
return get(this@SymbolIndexer.indexOf(symbol))
}
public operator fun DoubleArray.get(symbol: Symbol): Double {
require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
return get(this@SymbolIndexer.indexOf(symbol))
}
public operator fun <T> Point<T>.get(symbol: Symbol): T {
require(size == symbols.size) { "The input buffer size for indexer should be ${symbols.size} but $size found" }
return get(this@SymbolIndexer.indexOf(symbol))
}
public fun DoubleArray.toMap(): Map<Symbol, Double> {
require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
return symbols.indices.associate { symbols[it] to get(it) }
}
public operator fun <T> Structure2D<T>.get(rowSymbol: Symbol, columnSymbol: Symbol): T =
get(indexOf(rowSymbol), indexOf(columnSymbol))
public fun <T> Map<Symbol, T>.toList(): List<T> = symbols.map { getValue(it) }
public fun <T> Map<Symbol, T>.toPoint(bufferFactory: BufferFactory<T>): Point<T> =
bufferFactory(symbols.size) { getValue(symbols[it]) }
public fun Map<Symbol, Double>.toDoubleArray(): DoubleArray = DoubleArray(symbols.size) { getValue(symbols[it]) }
}
public inline class SimpleSymbolIndexer(override val symbols: List<Symbol>) : SymbolIndexer
/**
* Execute the block with symbol indexer based on given symbol order
*/
public inline fun <R> withSymbols(vararg symbols: Symbol, block: SymbolIndexer.() -> R): R =
with(SimpleSymbolIndexer(symbols.toList()), block)
public inline fun <R> withSymbols(symbols: Collection<Symbol>, block: SymbolIndexer.() -> R): R =
with(SimpleSymbolIndexer(symbols.toList()), block)

1
kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/expressionBuilders.kt
View File

@ -7,6 +7,7 @@ import kscience.kmath.operations.Space
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
/**
* Creates a functional expression with this [Space].
*/

41
kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
View File

@ -9,8 +9,8 @@ import kscience.kmath.structures.*
*/
public class BufferMatrixContext<T : Any, R : Ring<T>>(
public override val elementContext: R,
private val bufferFactory: BufferFactory<T>
) : GenericMatrixContext<T, R> {
private val bufferFactory: BufferFactory<T>,
) : GenericMatrixContext<T, R, BufferMatrix<T>> {
public 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)
@ -22,15 +22,15 @@ public class BufferMatrixContext<T : Any, R : Ring<T>>(
}
@Suppress("OVERRIDE_BY_INLINE")
public object RealMatrixContext : GenericMatrixContext<Double, RealField> {
public object RealMatrixContext : GenericMatrixContext<Double, RealField, BufferMatrix<Double>> {
public override val elementContext: RealField
get() = RealField
public override inline fun produce(
rows: Int,
columns: Int,
initializer: (i: Int, j: Int) -> Double
): Matrix<Double> {
initializer: (i: Int, j: Int) -> Double,
): BufferMatrix<Double> {
val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
return BufferMatrix(rows, columns, buffer)
}
@ -43,15 +43,15 @@ public class BufferMatrix<T : Any>(
public override val rowNum: Int,
public override val colNum: Int,
public val buffer: Buffer<out T>,
public override val features: Set<MatrixFeature> = emptySet()
public override val features: Set<MatrixFeature> = emptySet(),
) : FeaturedMatrix<T> {
override val shape: IntArray
get() = intArrayOf(rowNum, colNum)
init {
require(buffer.size == rowNum * colNum) { "Dimension mismatch for matrix structure" }
}
override val shape: IntArray get() = intArrayOf(rowNum, colNum)
public override fun suggestFeature(vararg features: MatrixFeature): BufferMatrix<T> =
BufferMatrix(rowNum, colNum, buffer, this.features + features)
@ -86,28 +86,3 @@ public class BufferMatrix<T : Any>(
else "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)"
}
}
/**
* Optimized dot product for real matrices
*/
public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
val array = DoubleArray(this.rowNum * other.colNum)
//convert to array to insure there is not memory indirection
fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
array
else
DoubleArray(size) { get(it) }
val a = this.buffer.unsafeArray()
val b = other.buffer.unsafeArray()
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 = RealBuffer(array)
return BufferMatrix(rowNum, other.colNum, buffer)
}

11
kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
View File

@ -27,9 +27,8 @@ public interface FeaturedMatrix<T : Any> : Matrix<T> {
public inline fun Structure2D.Companion.real(
rows: Int,
columns: Int,
initializer: (Int, Int) -> Double
): Matrix<Double> =
MatrixContext.real.produce(rows, columns, initializer)
initializer: (Int, Int) -> Double,
): BufferMatrix<Double> = MatrixContext.real.produce(rows, columns, initializer)
/**
* Build a square matrix from given elements.
@ -58,7 +57,7 @@ public inline fun <reified T : Any> Matrix<*>.getFeature(): T? =
/**
* Diagonal matrix of ones. The matrix is virtual no actual matrix is created
*/
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
VirtualMatrix(rows, columns, DiagonalFeature) { i, j ->
if (i == j) elementContext.one else elementContext.zero
}
@ -67,7 +66,7 @@ public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, colu
/**
* A virtual matrix of zeroes
*/
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
VirtualMatrix(rows, columns) { _, _ -> elementContext.zero }
public class TransposedFeature<T : Any>(public val original: Matrix<T>) : MatrixFeature
@ -82,5 +81,3 @@ public fun <T : Any> Matrix<T>.transpose(): Matrix<T> {
setOf(TransposedFeature(this))
) { i, j -> get(j, i) }
}
public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = with(MatrixContext.real) { dot(other) }

86
kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
View File

@ -1,25 +1,18 @@
package kscience.kmath.linear
import kscience.kmath.operations.Field
import kscience.kmath.operations.RealField
import kscience.kmath.operations.Ring
import kscience.kmath.operations.invoke
import kscience.kmath.structures.BufferAccessor2D
import kscience.kmath.structures.Matrix
import kscience.kmath.structures.Structure2D
import kotlin.reflect.KClass
import kscience.kmath.operations.*
import kscience.kmath.structures.*
/**
* Common implementation of [LUPDecompositionFeature]
*/
public class LUPDecomposition<T : Any>(
public val context: GenericMatrixContext<T, out Field<T>>,
public val context: MatrixContext<T, FeaturedMatrix<T>>,
public val elementContext: Field<T>,
public val lu: Structure2D<T>,
public val pivot: IntArray,
private val even: Boolean
private val even: Boolean,
) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
public val elementContext: Field<T>
get() = context.elementContext
/**
* Returns the matrix L of the decomposition.
@ -64,23 +57,25 @@ public class LUPDecomposition<T : Any>(
}
public fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T): T =
@PublishedApi
internal fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, *>.abs(value: T): T =
if (value > elementContext.zero) value else elementContext { -value }
/**
* Create a lup decomposition of generic matrix
* Create a lup decomposition of generic matrix.
*/
public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
type: KClass<T>,
public fun <T : Comparable<T>> MatrixContext<T, FeaturedMatrix<T>>.lup(
factory: MutableBufferFactory<T>,
elementContext: Field<T>,
matrix: Matrix<T>,
checkSingular: (T) -> Boolean
checkSingular: (T) -> Boolean,
): LUPDecomposition<T> {
require(matrix.rowNum == matrix.colNum) { "LU decomposition supports only square matrices" }
val m = matrix.colNum
val pivot = IntArray(matrix.rowNum)
//TODO just waits for KEEP-176
BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
elementContext {
val lu = create(matrix)
@ -112,14 +107,14 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
luRow[col] = sum
// maintain best permutation choice
if (this@lup.abs(sum) > largest) {
largest = this@lup.abs(sum)
if (abs(sum) > largest) {
largest = abs(sum)
max = row
}
}
// Singularity check
check(!checkSingular(this@lup.abs(lu[max, col]))) { "The matrix is singular" }
check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" }
// Pivot if necessary
if (max != col) {
@ -143,23 +138,23 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
for (row in col + 1 until m) lu[row, col] /= luDiag
}
return LUPDecomposition(this@lup, lu.collect(), pivot, even)
return LUPDecomposition(this@lup, elementContext, lu.collect(), pivot, even)
}
}
}
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.lup(
matrix: Matrix<T>,
checkSingular: (T) -> Boolean
): LUPDecomposition<T> = lup(T::class, matrix, checkSingular)
noinline checkSingular: (T) -> Boolean,
): LUPDecomposition<T> = lup(MutableBuffer.Companion::auto, elementContext, matrix, checkSingular)
public fun GenericMatrixContext<Double, RealField>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
lup(Double::class, matrix) { it < 1e-11 }
public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T>): Matrix<T> {
public fun <T : Any> LUPDecomposition<T>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<T> {
require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" }
BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
elementContext {
// Apply permutations to b
val bp = create { _, _ -> zero }
@ -201,27 +196,34 @@ public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T
}
}
public inline fun <reified T : Any> LUPDecomposition<T>.solve(matrix: Matrix<T>): Matrix<T> = solve(T::class, matrix)
public inline fun <reified T : Any> LUPDecomposition<T>.solveWithLUP(matrix: Matrix<T>): Matrix<T> =
solveWithLUP(MutableBuffer.Companion::auto, matrix)
/**
* Solve a linear equation **a*x = b**
* Solve a linear equation **a*x = b** using LUP decomposition
*/
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.solveWithLUP(
a: Matrix<T>,
b: Matrix<T>,
checkSingular: (T) -> Boolean
): Matrix<T> {
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
noinline checkSingular: (T) -> Boolean,
): FeaturedMatrix<T> {
// Use existing decomposition if it is provided by matrix
val decomposition = a.getFeature() ?: lup(T::class, a, checkSingular)
return decomposition.solve(T::class, b)
val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular)
return decomposition.solveWithLUP(bufferFactory, b)
}
public fun RealMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solve(a, b) { it < 1e-11 }
public fun RealMatrixContext.solveWithLUP(a: Matrix<Double>, b: Matrix<Double>): FeaturedMatrix<Double> =
solveWithLUP(a, b) { it < 1e-11 }
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.inverse(
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.inverseWithLUP(
matrix: Matrix<T>,
checkSingular: (T) -> Boolean
): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum), checkSingular)
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
noinline checkSingular: (T) -> Boolean,
): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
public fun RealMatrixContext.inverse(matrix: Matrix<Double>): Matrix<Double> =
solve(matrix, one(matrix.rowNum, matrix.colNum)) { it < 1e-11 }
/**
* Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
*/
public fun RealMatrixContext.inverseWithLUP(matrix: Matrix<Double>): FeaturedMatrix<Double> =
solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), Buffer.Companion::real) { it < 1e-11 }

37
kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
View File

@ -12,15 +12,16 @@ import kscience.kmath.structures.asSequence
/**
* Basic operations on matrices. Operates on [Matrix]
*/
public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Matrix<T>> {
/**
* Produce a matrix with this context and given dimensions
*/
public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
public override fun binaryOperation(operation: String, left: Matrix<T>, right: Matrix<T>): Matrix<T> = when (operation) {
@Suppress("UNCHECKED_CAST")
public override fun binaryOperation(operation: String, left: Matrix<T>, right: Matrix<T>): M = when (operation) {
"dot" -> left dot right
else -> super.binaryOperation(operation, left, right)
else -> super.binaryOperation(operation, left, right) as M
}
/**
@ -30,7 +31,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
* @param other the multiplier.
* @return the dot product.
*/
public infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
public infix fun Matrix<T>.dot(other: Matrix<T>): M
/**
* Computes the dot product of this matrix and a vector.
@ -48,7 +49,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
* @param value the multiplier.
* @receiver the product.
*/
public operator fun Matrix<T>.times(value: T): Matrix<T>
public operator fun Matrix<T>.times(value: T): M
/**
* Multiplies an element by a matrix of it.
@ -57,7 +58,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
* @param value the multiplier.
* @receiver the product.
*/
public operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
public operator fun T.times(m: Matrix<T>): M = m * this
public companion object {
/**
@ -70,18 +71,18 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
*/
public fun <T : Any, R : Ring<T>> buffered(
ring: R,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
): GenericMatrixContext<T, R> = BufferMatrixContext(ring, bufferFactory)
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
): GenericMatrixContext<T, R, BufferMatrix<T>> = BufferMatrixContext(ring, bufferFactory)
/**
* Automatic buffered matrix, unboxed if it is possible
*/
public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R, BufferMatrix<T>> =
buffered(ring, Buffer.Companion::auto)
}
}
public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
public interface GenericMatrixContext<T : Any, R : Ring<T>, out M : Matrix<T>> : MatrixContext<T, M> {
/**
* The ring context for matrix elements
*/
@ -92,7 +93,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
*/
public fun point(size: Int, initializer: (Int) -> T): Point<T>
public override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
public override infix fun Matrix<T>.dot(other: Matrix<T>): M {
//TODO add typed error
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
@ -113,10 +114,10 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
}
}
public override operator fun Matrix<T>.unaryMinus(): Matrix<T> =
public override operator fun Matrix<T>.unaryMinus(): M =
produce(rowNum, colNum) { i, j -> elementContext { -get(i, j) } }
public override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
public override fun add(a: Matrix<T>, b: Matrix<T>): M {
require(a.rowNum == b.rowNum && a.colNum == b.colNum) {
"Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]"
}
@ -124,7 +125,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
return produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] + b[i, j] } }
}
public override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
public override operator fun Matrix<T>.minus(b: Matrix<T>): M {
require(rowNum == b.rowNum && colNum == b.colNum) {
"Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]"
}
@ -132,11 +133,11 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } }
}
public override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
public override fun multiply(a: Matrix<T>, k: Number): M =
produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } }
public operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
public operator fun Number.times(matrix: FeaturedMatrix<T>): M = multiply(matrix, this)
public override operator fun Matrix<T>.times(value: T): Matrix<T> =
public override operator fun Matrix<T>.times(value: T): M =
produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } }
}

8
kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VectorSpace.kt
View File

@ -15,7 +15,7 @@ public interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
public val space: S
override val zero: Point<T> get() = produce { space.zero }
public fun produce(initializer: (Int) -> T): Point<T>
public fun produce(initializer: S.(Int) -> T): Point<T>
/**
* Produce a space-element of this vector space for expressions
@ -48,7 +48,7 @@ public interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
public fun <T : Any, S : Space<T>> buffered(
size: Int,
space: S,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
): BufferVectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
/**
@ -63,8 +63,8 @@ public interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
public class BufferVectorSpace<T : Any, S : Space<T>>(
override val size: Int,
override val space: S,
public val bufferFactory: BufferFactory<T>
public val bufferFactory: BufferFactory<T>,
) : VectorSpace<T, S> {
override fun produce(initializer: (Int) -> T): Buffer<T> = bufferFactory(size, initializer)
override fun produce(initializer: S.(Int) -> T): Buffer<T> = bufferFactory(size) { space.initializer(it) }
//override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
}

266
kmath-core/src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt
View File

@ -1,266 +0,0 @@
package kscience.kmath.misc
import kscience.kmath.linear.Point
import kscience.kmath.operations.*
import kscience.kmath.structures.asBuffer
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
/*
* Implementation of backward-mode automatic differentiation.
* Initial gist by Roman Elizarov: https://gist.github.com/elizarov/1ad3a8583e88cb6ea7a0ad09bb591d3d
*/
/**
* Differentiable variable with value and derivative of differentiation ([deriv]) result
* with respect to this variable.
*
* @param T the non-nullable type of value.
* @property value The value of this variable.
*/
public open class Variable<T : Any>(public val value: T)
/**
* Represents result of [deriv] call.
*
* @param T the non-nullable type of value.
* @param value the value of result.
* @property deriv The mapping of differentiated variables to their derivatives.
* @property context The field over [T].
*/
public class DerivationResult<T : Any>(
value: T,
public val deriv: Map<Variable<T>, T>,
public val context: Field<T>
) : Variable<T>(value) {
/**
* Returns derivative of [variable] or returns [Ring.zero] in [context].
*/
public fun deriv(variable: Variable<T>): T = deriv[variable] ?: context.zero
/**
* Computes the divergence.
*/
public fun div(): T = context { sum(deriv.values) }
/**
* Computes the gradient for variables in given order.
*/
public fun grad(vararg variables: Variable<T>): Point<T> {
check(variables.isNotEmpty()) { "Variable order is not provided for gradient construction" }
return variables.map(::deriv).asBuffer()
}
}
/**
* Runs differentiation and establishes [AutoDiffField] context inside the block of code.
*
* The partial derivatives are placed in argument `d` variable
*
* Example:
* ```
* val x = Variable(2) // define variable(s) and their values
* val y = deriv { sqr(x) + 5 * x + 3 } // write formulate in deriv context
* assertEquals(17.0, y.x) // the value of result (y)
* assertEquals(9.0, x.d) // dy/dx
* ```
*
* @param body the action in [AutoDiffField] context returning [Variable] to differentiate with respect to.
* @return the result of differentiation.
*/
public inline fun <T : Any, F : Field<T>> F.deriv(body: AutoDiffField<T, F>.() -> Variable<T>): DerivationResult<T> {
contract { callsInPlace(body, InvocationKind.EXACTLY_ONCE) }
return (AutoDiffContext(this)) {
val result = body()
result.d = context.one // computing derivative w.r.t result
runBackwardPass()
DerivationResult(result.value, derivatives, this@deriv)
}
}
/**
* Represents field in context of which functions can be derived.
*/
public abstract class AutoDiffField<T : Any, F : Field<T>> : Field<Variable<T>> {
public abstract val context: F
/**
* A variable accessing inner state of derivatives.
* Use this value in inner builders to avoid creating additional derivative bindings.
*/
public abstract var Variable<T>.d: T
/**
* Performs update of derivative after the rest of the formula in the back-pass.
*
* For example, implementation of `sin` function is:
*
* ```
* fun AD.sin(x: Variable): Variable = derive(Variable(sin(x.x)) { z -> // call derive with function result
* x.d += z.d * cos(x.x) // update derivative using chain rule and derivative of the function
* }
* ```
*/
public abstract fun <R> derive(value: R, block: F.(R) -> Unit): R
/**
*
*/
public abstract fun variable(value: T): Variable<T>
public inline fun variable(block: F.() -> T): Variable<T> = variable(context.block())
// Overloads for Double constants
override operator fun Number.plus(b: Variable<T>): Variable<T> =
derive(variable { this@plus.toDouble() * one + b.value }) { z ->
b.d += z.d
}
override operator fun Variable<T>.plus(b: Number): Variable<T> = b.plus(this)
override operator fun Number.minus(b: Variable<T>): Variable<T> =
derive(variable { this@minus.toDouble() * one - b.value }) { z -> b.d -= z.d }
override operator fun Variable<T>.minus(b: Number): Variable<T> =
derive(variable { this@minus.value - one * b.toDouble() }) { z -> this@minus.d += z.d }
}
/**
* Automatic Differentiation context class.
*/
@PublishedApi
internal class AutoDiffContext<T : Any, F : Field<T>>(override val context: F) : AutoDiffField<T, F>() {
// this stack contains pairs of blocks and values to apply them to
private var stack: Array<Any?> = arrayOfNulls<Any?>(8)
private var sp: Int = 0
val derivatives: MutableMap<Variable<T>, T> = hashMapOf()
override val zero: Variable<T> get() = Variable(context.zero)
override val one: Variable<T> get() = Variable(context.one)
/**
* A variable coupled with its derivative. For internal use only
*/
private class VariableWithDeriv<T : Any>(x: T, var d: T) : Variable<T>(x)
override fun variable(value: T): Variable<T> =
VariableWithDeriv(value, context.zero)
override var Variable<T>.d: T
get() = (this as? VariableWithDeriv)?.d ?: derivatives[this] ?: context.zero
set(value) = if (this is VariableWithDeriv) d = value else derivatives[this] = value
@Suppress("UNCHECKED_CAST")
override fun <R> derive(value: R, block: F.(R) -> Unit): R {
// save block to stack for backward pass
if (sp >= stack.size) stack = stack.copyOf(stack.size * 2)
stack[sp++] = block
stack[sp++] = value
return value
}
@Suppress("UNCHECKED_CAST")
fun runBackwardPass() {
while (sp > 0) {
val value = stack[--sp]
val block = stack[--sp] as F.(Any?) -> Unit
context.block(value)
}
}
// Basic math (+, -, *, /)
override fun add(a: Variable<T>, b: Variable<T>): Variable<T> = derive(variable { a.value + b.value }) { z ->
a.d += z.d
b.d += z.d
}
override fun multiply(a: Variable<T>, b: Variable<T>): Variable<T> = derive(variable { a.value * b.value }) { z ->
a.d += z.d * b.value
b.d += z.d * a.value
}
override fun divide(a: Variable<T>, b: Variable<T>): Variable<T> = derive(variable { a.value / b.value }) { z ->
a.d += z.d / b.value
b.d -= z.d * a.value / (b.value * b.value)
}
override fun multiply(a: Variable<T>, k: Number): Variable<T> = derive(variable { k.toDouble() * a.value }) { z ->
a.d += z.d * k.toDouble()
}
}
// Extensions for differentiation of various basic mathematical functions
// x ^ 2
public fun <T : Any, F : Field<T>> AutoDiffField<T, F>.sqr(x: Variable<T>): Variable<T> =
derive(variable { x.value * x.value }) { z -> x.d += z.d * 2 * x.value }
// x ^ 1/2
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.sqrt(x: Variable<T>): Variable<T> =
derive(variable { sqrt(x.value) }) { z -> x.d += z.d * 0.5 / z.value }
// x ^ y (const)
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.pow(x: Variable<T>, y: Double): Variable<T> =
derive(variable { power(x.value, y) }) { z -> x.d += z.d * y * power(x.value, y - 1) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.pow(x: Variable<T>, y: Int): Variable<T> =
pow(x, y.toDouble())
// exp(x)
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.exp(x: Variable<T>): Variable<T> =
derive(variable { exp(x.value) }) { z -> x.d += z.d * z.value }
// ln(x)
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.ln(x: Variable<T>): Variable<T> =
derive(variable { ln(x.value) }) { z -> x.d += z.d / x.value }
// x ^ y (any)
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.pow(x: Variable<T>, y: Variable<T>): Variable<T> =
exp(y * ln(x))
// sin(x)
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.sin(x: Variable<T>): Variable<T> =
derive(variable { sin(x.value) }) { z -> x.d += z.d * cos(x.value) }
// cos(x)
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.cos(x: Variable<T>): Variable<T> =
derive(variable { cos(x.value) }) { z -> x.d -= z.d * sin(x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.tan(x: Variable<T>): Variable<T> =
derive(variable { tan(x.value) }) { z ->
val c = cos(x.value)
x.d += z.d / (c * c)
}
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.asin(x: Variable<T>): Variable<T> =
derive(variable { asin(x.value) }) { z -> x.d += z.d / sqrt(one - x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.acos(x: Variable<T>): Variable<T> =
derive(variable { acos(x.value) }) { z -> x.d -= z.d / sqrt(one - x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.atan(x: Variable<T>): Variable<T> =
derive(variable { atan(x.value) }) { z -> x.d += z.d / (one + x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.sinh(x: Variable<T>): Variable<T> =
derive(variable { sin(x.value) }) { z -> x.d += z.d * cosh(x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.cosh(x: Variable<T>): Variable<T> =
derive(variable { cos(x.value) }) { z -> x.d += z.d * sinh(x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.tanh(x: Variable<T>): Variable<T> =
derive(variable { tan(x.value) }) { z ->
val c = cosh(x.value)
x.d += z.d / (c * c)
}
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.asinh(x: Variable<T>): Variable<T> =
derive(variable { asinh(x.value) }) { z -> x.d += z.d / sqrt(one + x.value * x.value) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.acosh(x: Variable<T>): Variable<T> =
derive(variable { acosh(x.value) }) { z -> x.d += z.d / (sqrt((x.value - one) * (x.value + one))) }
public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.atanh(x: Variable<T>): Variable<T> =
derive(variable { atanh(x.value) }) { z -> x.d += z.d / (one - x.value * x.value) }

4
kmath-core/src/commonMain/kotlin/kscience/kmath/misc/annotations.kt Normal file
View File

@ -0,0 +1,4 @@
package kscience.kmath.misc
@RequiresOptIn("This API is unstable and could change in future", RequiresOptIn.Level.WARNING)
public annotation class UnstableKMathAPI

28
kmath-core/src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt
View File

@ -22,10 +22,22 @@ public interface Algebra<T> {
*/
public fun unaryOperation(operation: String, arg: T): T
/**
* Currying version of [unaryOperation]
*/
public fun unaryOperationFunction(operation: String): (T) -> T = { unaryOperation(operation, it) }
/**
* Dynamic call of binary operation [operation] on [left] and [right]
*/
public fun binaryOperation(operation: String, left: T, right: T): T
/**
* Curring version of [binaryOperation]
*/
public fun binaryOperationFunction(operation: String): (left: T, right: T) -> T = { left, right ->
binaryOperation(operation, left, right)
}
}
/**
@ -45,11 +57,27 @@ public interface NumericAlgebra<T> : Algebra<T> {
public fun leftSideNumberOperation(operation: String, left: Number, right: T): T =
binaryOperation(operation, number(left), right)
/**
* Curring version of [leftSideNumberOperation]
*/
public fun leftSideNumberOperationFunction(operation: String): (left: Number, right: T) -> T =
{ left: Number, right: T ->
leftSideNumberOperation(operation, left, right)
}
/**
* Dynamic call of binary operation [operation] on [left] and [right] where right element is [Number].
*/
public fun rightSideNumberOperation(operation: String, left: T, right: Number): T =
leftSideNumberOperation(operation, right, left)
/**
* Curring version of [rightSideNumberOperation]
*/
public fun rightSideNumberOperationFunction(operation: String): (left: T, right: Number) -> T =
{ left: T, right: Number ->
rightSideNumberOperation(operation, left, right)
}
}
/**

6
kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraElements.kt
View File

@ -74,9 +74,9 @@ public interface SpaceElement<T, I : SpaceElement<T, I, S>, S : Space<T>> : Math
/**
* The element of [Ring].
*
* @param T the type of space operation results.
* @param T the type of ring operation results.
* @param I self type of the element. Needed for static type checking.
* @param R the type of space.
* @param R the type of ring.
*/
public interface RingElement<T, I : RingElement<T, I, R>, R : Ring<T>> : SpaceElement<T, I, R> {
/**
@ -91,7 +91,7 @@ public interface RingElement<T, I : RingElement<T, I, R>, R : Ring<T>> : SpaceEl
/**
* The element of [Field].
*
* @param T the type of space operation results.
* @param T the type of field operation results.
* @param I self type of the element. Needed for static type checking.
* @param F the type of field.
*/

5
kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
View File

@ -38,6 +38,11 @@ public fun <T> Space<T>.average(data: Iterable<T>): T = sum(data) / data.count()
*/
public fun <T> Space<T>.average(data: Sequence<T>): T = sum(data) / data.count()
/**
* Absolute of the comparable [value]
*/
public fun <T : Comparable<T>> Space<T>.abs(value: T): T = if (value > zero) value else -value
/**
* Returns the sum of all elements in the iterable in provided space.
*

1
kmath-core/src/commonMain/kotlin/kscience/kmath/operations/Complex.kt
View File

@ -195,6 +195,7 @@ public data class Complex(val re: Double, val im: Double) : FieldElement<Complex
}
}
/**
* Creates a complex number with real part equal to this real.
*

22
kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
View File

@ -1,25 +1,31 @@
package kscience.kmath.structures
import kotlin.reflect.KClass
/**
* A context that allows to operate on a [MutableBuffer] as on 2d array
*/
public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val rowNum: Int, public val colNum: Int) {
internal class BufferAccessor2D<T : Any>(
public val rowNum: Int,
public val colNum: Int,
val factory: MutableBufferFactory<T>,
) {
public operator fun Buffer<T>.get(i: Int, j: Int): T = get(i + colNum * j)
public operator fun MutableBuffer<T>.set(i: Int, j: Int, value: T) {
set(i + colNum * j, value)
}
public inline fun create(init: (i: Int, j: Int) -> T): MutableBuffer<T> =
MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
public inline fun create(crossinline init: (i: Int, j: Int) -> T): MutableBuffer<T> =
factory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
public fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
//TODO optimize wrapper
public fun MutableBuffer<T>.collect(): Structure2D<T> =
NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D()
public fun MutableBuffer<T>.collect(): Structure2D<T> = NDStructure.build(
DefaultStrides(intArrayOf(rowNum, colNum)),
factory
) { (i, j) ->
get(i, j)
}.as2D()
public inner class Row(public val buffer: MutableBuffer<T>, public val rowIndex: Int) : MutableBuffer<T> {
override val size: Int get() = colNum
@ -30,7 +36,7 @@ public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val ro
buffer[rowIndex, index] = value
}
override fun copy(): MutableBuffer<T> = MutableBuffer.auto(type, colNum) { get(it) }
override fun copy(): MutableBuffer<T> = factory(colNum) { get(it) }
override operator fun iterator(): Iterator<T> = (0 until colNum).map(::get).iterator()
}

5
kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
View File

@ -102,6 +102,11 @@ public fun <T> Buffer<T>.asSequence(): Sequence<T> = Sequence(::iterator)
*/
public fun <T> Buffer<T>.asIterable(): Iterable<T> = Iterable(::iterator)
/**
* Converts this [Buffer] to a new [List]
*/
public fun <T> Buffer<T>.toList(): List<T> = asSequence().toList()
/**
* Returns an [IntRange] of the valid indices for this [Buffer].
*/

2
kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDAlgebra.kt
View File

@ -73,7 +73,7 @@ public interface NDAlgebra<T, C, N : NDStructure<T>> {
public fun check(vararg elements: N): Array<out N> = elements
.map(NDStructure<T>::shape)
.singleOrNull { !shape.contentEquals(it) }
?.let { throw ShapeMismatchException(shape, it) }
?.let<IntArray, Array<out N>> { throw ShapeMismatchException(shape, it) }
?: elements
/**

13
kmath-core/src/commonMain/kotlin/kscience/kmath/structures/RealBuffer.kt
View File

@ -5,19 +5,19 @@ package kscience.kmath.structures
*
* @property array the underlying array.
*/
@Suppress("OVERRIDE_BY_INLINE")
public inline class RealBuffer(public val array: DoubleArray) : MutableBuffer<Double> {
override val size: Int get() = array.size
override operator fun get(index: Int): Double = array[index]
override inline operator fun get(index: Int): Double = array[index]
override operator fun set(index: Int, value: Double) {
override inline operator fun set(index: Int, value: Double) {
array[index] = value
}
override operator fun iterator(): DoubleIterator = array.iterator()
override fun copy(): MutableBuffer<Double> =
RealBuffer(array.copyOf())
override fun copy(): RealBuffer = RealBuffer(array.copyOf())
}
/**
@ -34,6 +34,11 @@ public inline fun RealBuffer(size: Int, init: (Int) -> Double): RealBuffer = Rea
*/
public fun RealBuffer(vararg doubles: Double): RealBuffer = RealBuffer(doubles)
/**
* Simplified [RealBuffer] to array comparison
*/
public fun RealBuffer.contentEquals(vararg doubles: Double): Boolean = array.contentEquals(doubles)
/**
* Returns a [DoubleArray] containing all of the elements of this [MutableBuffer].
*/

22
kmath-core/src/commonTest/kotlin/kscience/kmath/expressions/ExpressionFieldTest.kt
View File

@ -6,19 +6,21 @@ import kscience.kmath.operations.RealField
import kscience.kmath.operations.invoke
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertFails
class ExpressionFieldTest {
val x by symbol
@Test
fun testExpression() {
val context = FunctionalExpressionField(RealField)
val expression = context {
val x = variable("x", 2.0)
val x by binding()
x * x + 2 * x + one
}
assertEquals(expression("x" to 1.0), 4.0)
assertEquals(expression(), 9.0)
assertEquals(expression(x to 1.0), 4.0)
assertFails { expression()}
}
@Test
@ -26,33 +28,33 @@ class ExpressionFieldTest {
val context = FunctionalExpressionField(ComplexField)
val expression = context {
val x = variable("x", Complex(2.0, 0.0))
val x = bind(x)
x * x + 2 * x + one
}
assertEquals(expression("x" to Complex(1.0, 0.0)), Complex(4.0, 0.0))
assertEquals(expression(), Complex(9.0, 0.0))
assertEquals(expression(x to Complex(1.0, 0.0)), Complex(4.0, 0.0))
//assertEquals(expression(), Complex(9.0, 0.0))
}
@Test
fun separateContext() {
fun <T> FunctionalExpressionField<T, *>.expression(): Expression<T> {
val x = variable("x")
val x by binding()
return x * x + 2 * x + one
}
val expression = FunctionalExpressionField(RealField).expression()
assertEquals(expression("x" to 1.0), 4.0)
assertEquals(expression(x to 1.0), 4.0)
}
@Test
fun valueExpression() {
val expressionBuilder: FunctionalExpressionField<Double, *>.() -> Expression<Double> = {
val x = variable("x")
val x by binding()
x * x + 2 * x + one
}
val expression = FunctionalExpressionField(RealField).expressionBuilder()
assertEquals(expression("x" to 1.0), 4.0)
assertEquals(expression(x to 1.0), 4.0)
}
}

285
kmath-core/src/commonTest/kotlin/kscience/kmath/expressions/SimpleAutoDiffTest.kt Normal file
View File

@ -0,0 +1,285 @@
package kscience.kmath.expressions
import kscience.kmath.operations.RealField
import kscience.kmath.structures.asBuffer
import kotlin.math.E
import kotlin.math.PI
import kotlin.math.pow
import kotlin.math.sqrt
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
class SimpleAutoDiffTest {
fun dx(
xBinding: Pair<Symbol, Double>,
body: SimpleAutoDiffField<Double, RealField>.(x: AutoDiffValue<Double>) -> AutoDiffValue<Double>,
): DerivationResult<Double> = RealField.simpleAutoDiff(xBinding) { body(bind(xBinding.first)) }
fun dxy(
xBinding: Pair<Symbol, Double>,
yBinding: Pair<Symbol, Double>,
body: SimpleAutoDiffField<Double, RealField>.(x: AutoDiffValue<Double>, y: AutoDiffValue<Double>) -> AutoDiffValue<Double>,
): DerivationResult<Double> = RealField.simpleAutoDiff(xBinding, yBinding) {
body(bind(xBinding.first), bind(yBinding.first))
}
fun diff(block: SimpleAutoDiffField<Double, RealField>.() -> AutoDiffValue<Double>): SimpleAutoDiffExpression<Double, RealField> {
return SimpleAutoDiffExpression(RealField, block)
}
val x by symbol
val y by symbol
val z by symbol
@Test
fun testPlusX2() {
val y = RealField.simpleAutoDiff(x to 3.0) {
// diff w.r.t this x at 3
val x = bind(x)
x + x
}
assertEquals(6.0, y.value) // y = x + x = 6
assertEquals(2.0, y.derivative(x)) // dy/dx = 2
}
@Test
fun testPlusX2Expr() {
val expr = diff {
val x = bind(x)
x + x
}
assertEquals(6.0, expr(x to 3.0)) // y = x + x = 6
assertEquals(2.0, expr.derivative(x)(x to 3.0)) // dy/dx = 2
}
@Test
fun testPlus() {
// two variables
val z = RealField.simpleAutoDiff(x to 2.0, y to 3.0) {
val x = bind(x)
val y = bind(y)
x + y
}
assertEquals(5.0, z.value) // z = x + y = 5
assertEquals(1.0, z.derivative(x)) // dz/dx = 1
assertEquals(1.0, z.derivative(y)) // dz/dy = 1
}
@Test
fun testMinus() {
// two variables
val z = RealField.simpleAutoDiff(x to 7.0, y to 3.0) {
val x = bind(x)
val y = bind(y)
x - y
}
assertEquals(4.0, z.value) // z = x - y = 4
assertEquals(1.0, z.derivative(x)) // dz/dx = 1
assertEquals(-1.0, z.derivative(y)) // dz/dy = -1
}
@Test
fun testMulX2() {
val y = dx(x to 3.0) { x ->
// diff w.r.t this x at 3
x * x
}
assertEquals(9.0, y.value) // y = x * x = 9
assertEquals(6.0, y.derivative(x)) // dy/dx = 2 * x = 7
}
@Test
fun testSqr() {
val y = dx(x to 3.0) { x -> sqr(x) }
assertEquals(9.0, y.value) // y = x ^ 2 = 9
assertEquals(6.0, y.derivative(x)) // dy/dx = 2 * x = 7
}
@Test
fun testSqrSqr() {
val y = dx(x to 2.0) { x -> sqr(sqr(x)) }
assertEquals(16.0, y.value) // y = x ^ 4 = 16
assertEquals(32.0, y.derivative(x)) // dy/dx = 4 * x^3 = 32
}
@Test
fun testX3() {
val y = dx(x to 2.0) { x ->
// diff w.r.t this x at 2
x * x * x
}
assertEquals(8.0, y.value) // y = x * x * x = 8
assertEquals(12.0, y.derivative(x)) // dy/dx = 3 * x * x = 12
}
@Test
fun testDiv() {
val z = dxy(x to 5.0, y to 2.0) { x, y ->
x / y
}
assertEquals(2.5, z.value) // z = x / y = 2.5
assertEquals(0.5, z.derivative(x)) // dz/dx = 1 / y = 0.5
assertEquals(-1.25, z.derivative(y)) // dz/dy = -x / y^2 = -1.25
}
@Test
fun testPow3() {
val y = dx(x to 2.0) { x ->
// diff w.r.t this x at 2
pow(x, 3)
}
assertEquals(8.0, y.value) // y = x ^ 3 = 8
assertEquals(12.0, y.derivative(x)) // dy/dx = 3 * x ^ 2 = 12
}
@Test
fun testPowFull() {
val z = dxy(x to 2.0, y to 3.0) { x, y ->
pow(x, y)
}
assertApprox(8.0, z.value) // z = x ^ y = 8
assertApprox(12.0, z.derivative(x)) // dz/dx = y * x ^ (y - 1) = 12
assertApprox(8.0 * kotlin.math.ln(2.0), z.derivative(y)) // dz/dy = x ^ y * ln(x)
}
@Test
fun testFromPaper() {
val y = dx(x to 3.0) { x -> 2 * x + x * x * x }
assertEquals(33.0, y.value) // y = 2 * x + x * x * x = 33
assertEquals(29.0, y.derivative(x)) // dy/dx = 2 + 3 * x * x = 29
}
@Test
fun testInnerVariable() {
val y = dx(x to 1.0) { x ->
const(1.0) * x
}
assertEquals(1.0, y.value) // y = x ^ n = 1
assertEquals(1.0, y.derivative(x)) // dy/dx = n * x ^ (n - 1) = n - 1
}
@Test
fun testLongChain() {
val n = 10_000
val y = dx(x to 1.0) { x ->
var res = const(1.0)
for (i in 1..n) res *= x
res
}
assertEquals(1.0, y.value) // y = x ^ n = 1
assertEquals(n.toDouble(), y.derivative(x)) // dy/dx = n * x ^ (n - 1) = n - 1
}
@Test
fun testExample() {
val y = dx(x to 2.0) { x -> sqr(x) + 5 * x + 3 }
assertEquals(17.0, y.value) // the value of result (y)
assertEquals(9.0, y.derivative(x)) // dy/dx
}
@Test
fun testSqrt() {
val y = dx(x to 16.0) { x -> sqrt(x) }
assertEquals(4.0, y.value) // y = x ^ 1/2 = 4
assertEquals(1.0 / 8, y.derivative(x)) // dy/dx = 1/2 / x ^ 1/4 = 1/8
}
@Test
fun testSin() {
val y = dx(x to PI / 6.0) { x -> sin(x) }
assertApprox(0.5, y.value) // y = sin(PI/6) = 0.5
assertApprox(sqrt(3.0) / 2, y.derivative(x)) // dy/dx = cos(pi/6) = sqrt(3)/2
}
@Test
fun testCos() {
val y = dx(x to PI / 6) { x -> cos(x) }
assertApprox(sqrt(3.0) / 2, y.value) //y = cos(pi/6) = sqrt(3)/2
assertApprox(-0.5, y.derivative(x)) // dy/dx = -sin(pi/6) = -0.5
}
@Test
fun testTan() {
val y = dx(x to PI / 6) { x -> tan(x) }
assertApprox(1.0 / sqrt(3.0), y.value) // y = tan(pi/6) = 1/sqrt(3)
assertApprox(4.0 / 3.0, y.derivative(x)) // dy/dx = sec(pi/6)^2 = 4/3
}
@Test
fun testAsin() {
val y = dx(x to PI / 6) { x -> asin(x) }
assertApprox(kotlin.math.asin(PI / 6.0), y.value) // y = asin(pi/6)
assertApprox(6.0 / sqrt(36 - PI * PI), y.derivative(x)) // dy/dx = 6/sqrt(36-pi^2)
}
@Test
fun testAcos() {
val y = dx(x to PI / 6) { x -> acos(x) }
assertApprox(kotlin.math.acos(PI / 6.0), y.value) // y = acos(pi/6)
assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.derivative(x)) // dy/dx = -6/sqrt(36-pi^2)
}
@Test
fun testAtan() {
val y = dx(x to PI / 6) { x -> atan(x) }
assertApprox(kotlin.math.atan(PI / 6.0), y.value) // y = atan(pi/6)
assertApprox(36.0 / (36.0 + PI * PI), y.derivative(x)) // dy/dx = 36/(36+pi^2)
}
@Test
fun testSinh() {
val y = dx(x to 0.0) { x -> sinh(x) }
assertApprox(kotlin.math.sinh(0.0), y.value) // y = sinh(0)
assertApprox(kotlin.math.cosh(0.0), y.derivative(x)) // dy/dx = cosh(0)
}
@Test
fun testCosh() {
val y = dx(x to 0.0) { x -> cosh(x) }
assertApprox(1.0, y.value) //y = cosh(0)
assertApprox(0.0, y.derivative(x)) // dy/dx = sinh(0)
}
@Test
fun testTanh() {
val y = dx(x to 1.0) { x -> tanh(x) }
assertApprox((E * E - 1) / (E * E + 1), y.value) // y = tanh(pi/6)
assertApprox(1.0 / kotlin.math.cosh(1.0).pow(2), y.derivative(x)) // dy/dx = sech(pi/6)^2
}
@Test
fun testAsinh() {
val y = dx(x to PI / 6) { x -> asinh(x) }
assertApprox(kotlin.math.asinh(PI / 6.0), y.value) // y = asinh(pi/6)
assertApprox(6.0 / sqrt(36 + PI * PI), y.derivative(x)) // dy/dx = 6/sqrt(pi^2+36)
}
@Test
fun testAcosh() {
val y = dx(x to PI / 6) { x -> acosh(x) }
assertApprox(kotlin.math.acosh(PI / 6.0), y.value) // y = acosh(pi/6)
assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.derivative(x)) // dy/dx = -6/sqrt(36-pi^2)
}
@Test
fun testAtanh() {
val y = dx(x to PI / 6) { x -> atanh(x) }
assertApprox(kotlin.math.atanh(PI / 6.0), y.value) // y = atanh(pi/6)
assertApprox(-36.0 / (PI * PI - 36.0), y.derivative(x)) // dy/dx = -36/(pi^2-36)
}
@Test
fun testDivGrad() {
val res = dxy(x to 1.0, y to 2.0) { x, y -> x * x + y * y }
assertEquals(6.0, res.div())
assertTrue(res.grad(x, y).contentEquals(doubleArrayOf(2.0, 4.0).asBuffer()))
}
private fun assertApprox(a: Double, b: Double) {
if ((a - b) > 1e-10) assertEquals(a, b)
}
}

3
kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
View File

@ -1,5 +1,6 @@
package kscience.kmath.linear
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
import kscience.kmath.structures.NDStructure
import kscience.kmath.structures.as2D
@ -38,7 +39,7 @@ class MatrixTest {
infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
var res = this
repeat(power - 1) {
res = res dot this
res = RealMatrixContext.invoke { res dot this@pow }
}
return res
}

4
kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
View File

@ -9,7 +9,7 @@ class RealLUSolverTest {
@Test
fun testInvertOne() {
val matrix = MatrixContext.real.one(2, 2)
val inverted = MatrixContext.real.inverse(matrix)
val inverted = MatrixContext.real.inverseWithLUP(matrix)
assertEquals(matrix, inverted)
}
@ -37,7 +37,7 @@ class RealLUSolverTest {
1.0, 3.0
)
val inverted = MatrixContext.real.inverse(matrix)
val inverted = MatrixContext.real.inverseWithLUP(matrix)
val expected = Matrix.square(
0.375, -0.125,

261
kmath-core/src/commonTest/kotlin/kscience/kmath/misc/AutoDiffTest.kt
View File

@ -1,261 +0,0 @@
package kscience.kmath.misc
import kscience.kmath.operations.RealField
import kscience.kmath.structures.asBuffer
import kotlin.math.PI
import kotlin.math.pow
import kotlin.math.sqrt
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
class AutoDiffTest {
inline fun deriv(body: AutoDiffField<Double, RealField>.() -> Variable<Double>): DerivationResult<Double> =
RealField.deriv(body)
@Test
fun testPlusX2() {
val x = Variable(3.0) // diff w.r.t this x at 3
val y = deriv { x + x }
assertEquals(6.0, y.value) // y = x + x = 6
assertEquals(2.0, y.deriv(x)) // dy/dx = 2
}
@Test
fun testPlus() {
// two variables
val x = Variable(2.0)
val y = Variable(3.0)
val z = deriv { x + y }
assertEquals(5.0, z.value) // z = x + y = 5
assertEquals(1.0, z.deriv(x)) // dz/dx = 1
assertEquals(1.0, z.deriv(y)) // dz/dy = 1
}
@Test
fun testMinus() {
// two variables
val x = Variable(7.0)
val y = Variable(3.0)
val z = deriv { x - y }
assertEquals(4.0, z.value) // z = x - y = 4
assertEquals(1.0, z.deriv(x)) // dz/dx = 1
assertEquals(-1.0, z.deriv(y)) // dz/dy = -1
}
@Test
fun testMulX2() {
val x = Variable(3.0) // diff w.r.t this x at 3
val y = deriv { x * x }
assertEquals(9.0, y.value) // y = x * x = 9
assertEquals(6.0, y.deriv(x)) // dy/dx = 2 * x = 7
}
@Test
fun testSqr() {
val x = Variable(3.0)
val y = deriv { sqr(x) }
assertEquals(9.0, y.value) // y = x ^ 2 = 9
assertEquals(6.0, y.deriv(x)) // dy/dx = 2 * x = 7
}
@Test
fun testSqrSqr() {
val x = Variable(2.0)
val y = deriv { sqr(sqr(x)) }
assertEquals(16.0, y.value) // y = x ^ 4 = 16
assertEquals(32.0, y.deriv(x)) // dy/dx = 4 * x^3 = 32
}
@Test
fun testX3() {
val x = Variable(2.0) // diff w.r.t this x at 2
val y = deriv { x * x * x }
assertEquals(8.0, y.value) // y = x * x * x = 8
assertEquals(12.0, y.deriv(x)) // dy/dx = 3 * x * x = 12
}
@Test
fun testDiv() {
val x = Variable(5.0)
val y = Variable(2.0)
val z = deriv { x / y }
assertEquals(2.5, z.value) // z = x / y = 2.5
assertEquals(0.5, z.deriv(x)) // dz/dx = 1 / y = 0.5
assertEquals(-1.25, z.deriv(y)) // dz/dy = -x / y^2 = -1.25
}
@Test
fun testPow3() {
val x = Variable(2.0) // diff w.r.t this x at 2
val y = deriv { pow(x, 3) }
assertEquals(8.0, y.value) // y = x ^ 3 = 8
assertEquals(12.0, y.deriv(x)) // dy/dx = 3 * x ^ 2 = 12
}
@Test
fun testPowFull() {
val x = Variable(2.0)
val y = Variable(3.0)
val z = deriv { pow(x, y) }
assertApprox(8.0, z.value) // z = x ^ y = 8
assertApprox(12.0, z.deriv(x)) // dz/dx = y * x ^ (y - 1) = 12
assertApprox(8.0 * kotlin.math.ln(2.0), z.deriv(y)) // dz/dy = x ^ y * ln(x)
}
@Test
fun testFromPaper() {
val x = Variable(3.0)
val y = deriv { 2 * x + x * x * x }
assertEquals(33.0, y.value) // y = 2 * x + x * x * x = 33
assertEquals(29.0, y.deriv(x)) // dy/dx = 2 + 3 * x * x = 29
}
@Test
fun testInnerVariable() {
val x = Variable(1.0)
val y = deriv {
Variable(1.0) * x
}
assertEquals(1.0, y.value) // y = x ^ n = 1
assertEquals(1.0, y.deriv(x)) // dy/dx = n * x ^ (n - 1) = n - 1
}
@Test
fun testLongChain() {
val n = 10_000
val x = Variable(1.0)
val y = deriv {
var res = Variable(1.0)
for (i in 1..n) res *= x
res
}
assertEquals(1.0, y.value) // y = x ^ n = 1
assertEquals(n.toDouble(), y.deriv(x)) // dy/dx = n * x ^ (n - 1) = n - 1
}
@Test
fun testExample() {
val x = Variable(2.0)
val y = deriv { sqr(x) + 5 * x + 3 }
assertEquals(17.0, y.value) // the value of result (y)
assertEquals(9.0, y.deriv(x)) // dy/dx
}
@Test
fun testSqrt() {
val x = Variable(16.0)
val y = deriv { sqrt(x) }
assertEquals(4.0, y.value) // y = x ^ 1/2 = 4
assertEquals(1.0 / 8, y.deriv(x)) // dy/dx = 1/2 / x ^ 1/4 = 1/8
}
@Test
fun testSin() {
val x = Variable(PI / 6.0)
val y = deriv { sin(x) }
assertApprox(0.5, y.value) // y = sin(PI/6) = 0.5
assertApprox(sqrt(3.0) / 2, y.deriv(x)) // dy/dx = cos(pi/6) = sqrt(3)/2
}
@Test
fun testCos() {
val x = Variable(PI / 6)
val y = deriv { cos(x) }
assertApprox(sqrt(3.0) / 2, y.value) //y = cos(pi/6) = sqrt(3)/2
assertApprox(-0.5, y.deriv(x)) // dy/dx = -sin(pi/6) = -0.5
}
@Test
fun testTan() {
val x = Variable(PI / 6)
val y = deriv { tan(x) }
assertApprox(1.0 / sqrt(3.0), y.value) // y = tan(pi/6) = 1/sqrt(3)
assertApprox(4.0 / 3.0, y.deriv(x)) // dy/dx = sec(pi/6)^2 = 4/3
}
@Test
fun testAsin() {
val x = Variable(PI / 6)
val y = deriv { asin(x) }
assertApprox(kotlin.math.asin(PI / 6.0), y.value) // y = asin(pi/6)
assertApprox(6.0 / sqrt(36 - PI * PI), y.deriv(x)) // dy/dx = 6/sqrt(36-pi^2)
}
@Test
fun testAcos() {
val x = Variable(PI / 6)
val y = deriv { acos(x) }
assertApprox(kotlin.math.acos(PI / 6.0), y.value) // y = acos(pi/6)
assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.deriv(x)) // dy/dx = -6/sqrt(36-pi^2)
}
@Test
fun testAtan() {
val x = Variable(PI / 6)
val y = deriv { atan(x) }
assertApprox(kotlin.math.atan(PI / 6.0), y.value) // y = atan(pi/6)
assertApprox(36.0 / (36.0 + PI * PI), y.deriv(x)) // dy/dx = 36/(36+pi^2)
}
@Test
fun testSinh() {
val x = Variable(0.0)
val y = deriv { sinh(x) }
assertApprox(kotlin.math.sinh(0.0), y.value) // y = sinh(0)
assertApprox(kotlin.math.cosh(0.0), y.deriv(x)) // dy/dx = cosh(0)
}
@Test
fun testCosh() {
val x = Variable(0.0)
val y = deriv { cosh(x) }
assertApprox(1.0, y.value) //y = cosh(0)
assertApprox(0.0, y.deriv(x)) // dy/dx = sinh(0)
}
@Test
fun testTanh() {
val x = Variable(PI / 6)
val y = deriv { tanh(x) }
assertApprox(1.0 / sqrt(3.0), y.value) // y = tanh(pi/6)
assertApprox(1.0 / kotlin.math.cosh(PI / 6.0).pow(2), y.deriv(x)) // dy/dx = sech(pi/6)^2
}
@Test
fun testAsinh() {
val x = Variable(PI / 6)
val y = deriv { asinh(x) }
assertApprox(kotlin.math.asinh(PI / 6.0), y.value) // y = asinh(pi/6)
assertApprox(6.0 / sqrt(36 + PI * PI), y.deriv(x)) // dy/dx = 6/sqrt(pi^2+36)
}
@Test
fun testAcosh() {
val x = Variable(PI / 6)
val y = deriv { acosh(x) }
assertApprox(kotlin.math.acosh(PI / 6.0), y.value) // y = acosh(pi/6)
assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.deriv(x)) // dy/dx = -6/sqrt(36-pi^2)
}
@Test
fun testAtanh() {
val x = Variable(PI / 6.0)
val y = deriv { atanh(x) }
assertApprox(kotlin.math.atanh(PI / 6.0), y.value) // y = atanh(pi/6)
assertApprox(-36.0 / (PI * PI - 36.0), y.deriv(x)) // dy/dx = -36/(pi^2-36)
}
@Test
fun testDivGrad() {
val x = Variable(1.0)
val y = Variable(2.0)
val res = deriv { x * x + y * y }
assertEquals(6.0, res.div())
assertTrue(res.grad(x, y).contentEquals(doubleArrayOf(2.0, 4.0).asBuffer()))
}
private fun assertApprox(a: Double, b: Double) {
if ((a - b) > 1e-10) assertEquals(a, b)
}
}

2
kmath-coroutines/src/commonMain/kotlin/kscience/kmath/chains/Chain.kt
View File

@ -47,7 +47,7 @@ public fun <T> Iterator<T>.asChain(): Chain<T> = SimpleChain { next() }
public fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain()
/**
* A simple chain of independent tokens
* A simple chain of independent tokens. [fork] returns the same chain.
*/
public class SimpleChain<out R>(private val gen: suspend () -> R) : Chain<R> {
public override suspend fun next(): R = gen()

4
kmath-dimensions/build.gradle.kts
View File

@ -18,3 +18,7 @@ kotlin.sourceSets {
}
}
}
readme{
maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
}

4
kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
View File

@ -42,7 +42,7 @@ public interface DMatrix<T, R : Dimension, C : Dimension> : Structure2D<T> {
* An inline wrapper for a Matrix
*/
public inline class DMatrixWrapper<T, R : Dimension, C : Dimension>(
public val structure: Structure2D<T>
private val structure: Structure2D<T>
) : DMatrix<T, R, C> {
override val shape: IntArray get() = structure.shape
override operator fun get(i: Int, j: Int): T = structure[i, j]
@ -81,7 +81,7 @@ public inline class DPointWrapper<T, D : Dimension>(public val point: Point<T>)
/**
* Basic operations on dimension-safe matrices. Operates on [Matrix]
*/
public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri>) {
public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri, Matrix<T>>) {
public inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
require(rowNum == Dimension.dim<R>().toInt()) {
"Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum"

24
kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
View File

@ -1,23 +1,22 @@
package kscience.kmath.ejml
import org.ejml.simple.SimpleMatrix
import kscience.kmath.linear.MatrixContext
import kscience.kmath.linear.Point
import kscience.kmath.operations.Space
import kscience.kmath.operations.invoke
import kscience.kmath.structures.Matrix
import org.ejml.simple.SimpleMatrix
/**
* Converts this matrix to EJML one.
*/
public fun Matrix<Double>.toEjml(): EjmlMatrix =
if (this is EjmlMatrix) this else EjmlMatrixContext.produce(rowNum, colNum) { i, j -> get(i, j) }
/**
* Represents context of basic operations operating with [EjmlMatrix].
*
* @author Iaroslav Postovalov
*/
public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext<Double> {
/**
* Converts this matrix to EJML one.
*/
public fun Matrix<Double>.toEjml(): EjmlMatrix =
if (this is EjmlMatrix) this else produce(rowNum, colNum) { i, j -> get(i, j) }
public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
/**
* Converts this vector to EJML one.
@ -47,11 +46,10 @@ public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext
EjmlMatrix(toEjml().origin - b.toEjml().origin)
public override fun multiply(a: Matrix<Double>, k: Number): EjmlMatrix =
produce(a.rowNum, a.colNum) { i, j -> space { a[i, j] * k } }
produce(a.rowNum, a.colNum) { i, j -> a[i, j] * k.toDouble() }
public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix = EjmlMatrix(toEjml().origin.scale(value))
public companion object
public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
EjmlMatrix(toEjml().origin.scale(value))
}
/**

44
kmath-for-real/README.md Normal file
View File

@ -0,0 +1,44 @@
# Real number specialization module (`kmath-for-real`)
- [RealVector](src/commonMain/kotlin/kscience/kmath/real/RealVector.kt) : Numpy-like operations for Buffers/Points
- [RealMatrix](src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt) : Numpy-like operations for 2d real structures
- [grids](src/commonMain/kotlin/kscience/kmath/structures/grids.kt) : Uniform grid generators
> #### Artifact:
>
> This module artifact: `kscience.kmath:kmath-for-real:0.2.0-dev-4`.
>
> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-for-real/_latestVersion)
>
> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-for-real/_latestVersion)
>
> **Gradle:**
>
> ```gradle
> repositories {
> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" }
> maven { url 'https://dl.bintray.com/mipt-npm/kscience' }
> maven { url 'https://dl.bintray.com/mipt-npm/dev' }
> maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
> }
>
> dependencies {
> implementation 'kscience.kmath:kmath-for-real:0.2.0-dev-4'
> }
> ```
> **Gradle Kotlin DSL:**
>
> ```kotlin
> repositories {
> maven("https://dl.bintray.com/kotlin/kotlin-eap")
> maven("https://dl.bintray.com/mipt-npm/kscience")
> maven("https://dl.bintray.com/mipt-npm/dev")
> maven("https://dl.bintray.com/hotkeytlt/maven")
> }
>
> dependencies {
> implementation("kscience.kmath:kmath-for-real:0.2.0-dev-4")
> }
> ```

28
kmath-for-real/build.gradle.kts
View File

@ -7,3 +7,31 @@ kotlin.sourceSets.commonMain {
api(project(":kmath-core"))
}
}
readme {
description = """
Extension module that should be used to achieve numpy-like behavior.
All operations are specialized to work with `Double` numbers without declaring algebraic contexts.
One can still use generic algebras though.
""".trimIndent()
maturity = ru.mipt.npm.gradle.Maturity.EXPERIMENTAL
propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
feature(
id = "RealVector",
description = "Numpy-like operations for Buffers/Points",
ref = "src/commonMain/kotlin/kscience/kmath/real/RealVector.kt"
)
feature(
id = "RealMatrix",
description = "Numpy-like operations for 2d real structures",
ref = "src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt"
)
feature(
id = "grids",
description = "Uniform grid generators",
ref = "src/commonMain/kotlin/kscience/kmath/structures/grids.kt"
)
}

5
kmath-for-real/docs/README-TEMPLATE.md Normal file
View File

@ -0,0 +1,5 @@
# Real number specialization module (`kmath-for-real`)
${features}
${artifact}

97
kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/realMatrix.kt → kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
View File

@ -1,12 +1,14 @@
package kscience.kmath.real
import kscience.kmath.linear.FeaturedMatrix
import kscience.kmath.linear.MatrixContext
import kscience.kmath.linear.RealMatrixContext.elementContext
import kscience.kmath.linear.VirtualMatrix
import kscience.kmath.linear.inverseWithLUP
import kscience.kmath.misc.UnstableKMathAPI
import kscience.kmath.operations.invoke
import kscience.kmath.operations.sum
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.Matrix
import kscience.kmath.structures.RealBuffer
import kscience.kmath.structures.asIterable
import kotlin.math.pow
@ -23,7 +25,7 @@ import kotlin.math.pow
* Functions that help create a real (Double) matrix
*/
public typealias RealMatrix = Matrix<Double>
public typealias RealMatrix = FeaturedMatrix<Double>
public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum, initializer)
@ -36,7 +38,7 @@ public fun Sequence<DoubleArray>.toMatrix(): RealMatrix = toList().let {
MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] }
}
public fun Matrix<Double>.repeatStackVertical(n: Int): RealMatrix =
public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix =
VirtualMatrix(rowNum * n, colNum) { row, col ->
get(if (row == 0) 0 else row % rowNum, col)
}
@ -45,76 +47,65 @@ public fun Matrix<Double>.repeatStackVertical(n: Int): RealMatrix =
* Operations for matrix and real number
*/
public operator fun Matrix<Double>.times(double: Double): RealMatrix =
public operator fun RealMatrix.times(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col ->
this[row, col] * double
}
public operator fun Matrix<Double>.plus(double: Double): RealMatrix =
public operator fun RealMatrix.plus(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col ->
this[row, col] + double
}
public operator fun Matrix<Double>.minus(double: Double): RealMatrix =
public operator fun RealMatrix.minus(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col ->
this[row, col] - double
}
public operator fun Matrix<Double>.div(double: Double): RealMatrix =
public operator fun RealMatrix.div(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col ->
this[row, col] / double
}
public operator fun Double.times(matrix: Matrix<Double>): RealMatrix =
public operator fun Double.times(matrix: RealMatrix): RealMatrix =
MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
this * matrix[row, col]
}
public operator fun Double.plus(matrix: Matrix<Double>): RealMatrix =
public operator fun Double.plus(matrix: RealMatrix): RealMatrix =
MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
this + matrix[row, col]
}
public operator fun Double.minus(matrix: Matrix<Double>): RealMatrix =
public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
this - matrix[row, col]
}
// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'?
//operator fun Double.div(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
//operator fun Double.div(matrix: RealMatrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
// row, col -> matrix[row, col] / this
//}
/*
* Per-element (!) square and power operations
*/
public fun Matrix<Double>.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col ->
this[row, col].pow(2)
}
public fun Matrix<Double>.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j ->
this[i, j].pow(n)
}
/*
* Operations on two matrices (per-element!)
*/
public operator fun Matrix<Double>.times(other: Matrix<Double>): RealMatrix =
@UnstableKMathAPI
public operator fun RealMatrix.times(other: RealMatrix): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> this[row, col] * other[row, col] }
public operator fun Matrix<Double>.plus(other: Matrix<Double>): RealMatrix =
public operator fun RealMatrix.plus(other: RealMatrix): RealMatrix =
MatrixContext.real.add(this, other)
public operator fun Matrix<Double>.minus(other: Matrix<Double>): RealMatrix =
public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> this[row, col] - other[row, col] }
/*
* Operations on columns
*/
public inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): Matrix<Double> =
public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum + 1) { row, col ->
if (col < colNum)
this[row, col]
@ -122,28 +113,28 @@ public inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double
mapper(rows[row])
}
public fun Matrix<Double>.extractColumns(columnRange: IntRange): RealMatrix =
public fun RealMatrix.extractColumns(columnRange: IntRange): RealMatrix =
MatrixContext.real.produce(rowNum, columnRange.count()) { row, col ->
this[row, columnRange.first + col]
}
public fun Matrix<Double>.extractColumn(columnIndex: Int): RealMatrix =
public fun RealMatrix.extractColumn(columnIndex: Int): RealMatrix =
extractColumns(columnIndex..columnIndex)
public fun Matrix<Double>.sumByColumn(): RealBuffer = RealBuffer(colNum) { j ->
public fun RealMatrix.sumByColumn(): RealBuffer = RealBuffer(colNum) { j ->
val column = columns[j]
elementContext { sum(column.asIterable()) }
}
public fun Matrix<Double>.minByColumn(): RealBuffer = RealBuffer(colNum) { j ->
public fun RealMatrix.minByColumn(): RealBuffer = RealBuffer(colNum) { j ->
columns[j].asIterable().minOrNull() ?: error("Cannot produce min on empty column")
}
public fun Matrix<Double>.maxByColumn(): RealBuffer = RealBuffer(colNum) { j ->
public fun RealMatrix.maxByColumn(): RealBuffer = RealBuffer(colNum) { j ->
columns[j].asIterable().maxOrNull() ?: error("Cannot produce min on empty column")
}
public fun Matrix<Double>.averageByColumn(): RealBuffer = RealBuffer(colNum) { j ->
public fun RealMatrix.averageByColumn(): RealBuffer = RealBuffer(colNum) { j ->
columns[j].asIterable().average()
}
@ -151,7 +142,39 @@ public fun Matrix<Double>.averageByColumn(): RealBuffer = RealBuffer(colNum) { j
* Operations processing all elements
*/
public fun Matrix<Double>.sum(): Double = elements().map { (_, value) -> value }.sum()
public fun Matrix<Double>.min(): Double? = elements().map { (_, value) -> value }.minOrNull()
public fun Matrix<Double>.max(): Double? = elements().map { (_, value) -> value }.maxOrNull()
public fun Matrix<Double>.average(): Double = elements().map { (_, value) -> value }.average()
public fun RealMatrix.sum(): Double = elements().map { (_, value) -> value }.sum()
public fun RealMatrix.min(): Double? = elements().map { (_, value) -> value }.minOrNull()
public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.maxOrNull()
public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average()
public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { i, j ->
transform(get(i, j))
}
/**
* Inverse a square real matrix using LUP decomposition
*/
public fun RealMatrix.inverseWithLUP(): RealMatrix = MatrixContext.real.inverseWithLUP(this)
//extended operations
public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) }
public fun RealMatrix.pow(p: Int): RealMatrix = map { it.pow(p) }
public fun exp(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.exp(it) }
public fun sqrt(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sqrt(it) }
public fun RealMatrix.square(): RealMatrix = map { it.pow(2) }
public fun sin(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sin(it) }
public fun cos(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.cos(it) }
public fun tan(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.tan(it) }
public fun ln(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.ln(it) }
public fun log10(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.log10(it) }

90
kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt
View File

@ -1,46 +1,82 @@
package kscience.kmath.real
import kscience.kmath.linear.BufferVectorSpace
import kscience.kmath.linear.Point
import kscience.kmath.linear.VectorSpace
import kscience.kmath.operations.Norm
import kscience.kmath.operations.RealField
import kscience.kmath.operations.SpaceElement
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.RealBuffer
import kscience.kmath.structures.asBuffer
import kscience.kmath.structures.asIterable
import kotlin.math.pow
import kotlin.math.sqrt
public typealias RealPoint = Point<Double>
public fun DoubleArray.asVector(): RealVector = RealVector(asBuffer())
public fun List<Double>.asVector(): RealVector = RealVector(asBuffer())
public typealias RealVector = Point<Double>
public object VectorL2Norm : Norm<Point<out Number>, Double> {
override fun norm(arg: Point<out Number>): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble))
}
public inline class RealVector(private val point: Point<Double>) :
SpaceElement<RealPoint, RealVector, VectorSpace<Double, RealField>>, RealPoint {
public override val size: Int get() = point.size
public override val context: VectorSpace<Double, RealField> get() = space(point.size)
public operator fun Buffer.Companion.invoke(vararg doubles: Double): RealVector = doubles.asBuffer()
public override fun unwrap(): RealPoint = point
public override fun RealPoint.wrap(): RealVector = RealVector(this)
public override operator fun get(index: Int): Double = point[index]
public override operator fun iterator(): Iterator<Double> = point.iterator()
/**
* Fill the vector of given [size] with given [value]
*/
public fun Buffer.Companion.same(size: Int, value: Number): RealVector = real(size) { value.toDouble() }
public companion object {
private val spaceCache: MutableMap<Int, BufferVectorSpace<Double, RealField>> = hashMapOf()
// Transformation methods
public inline operator fun invoke(dim: Int, initializer: (Int) -> Double): RealVector =
RealVector(RealBuffer(dim, initializer))
public inline fun RealVector.map(transform: (Double) -> Double): RealVector =
Buffer.real(size) { transform(get(it)) }
public operator fun invoke(vararg values: Double): RealVector = values.asVector()
public inline fun RealVector.mapIndexed(transform: (index: Int, value: Double) -> Double): RealVector =
Buffer.real(size) { transform(it, get(it)) }
public fun space(dim: Int): BufferVectorSpace<Double, RealField> = spaceCache.getOrPut(dim) {
BufferVectorSpace(dim, RealField) { size, init -> Buffer.real(size, init) }
}
}
}
public operator fun RealVector.plus(other: RealVector): RealVector =
mapIndexed { index, value -> value + other[index] }
public operator fun RealVector.plus(number: Number): RealVector = map { it + number.toDouble() }
public operator fun Number.plus(vector: RealVector): RealVector = vector + this
public operator fun RealVector.unaryMinus(): Buffer<Double> = map { -it }
public operator fun RealVector.minus(other: RealVector): RealVector =
mapIndexed { index, value -> value - other[index] }
public operator fun RealVector.minus(number: Number): RealVector = map { it - number.toDouble() }
public operator fun Number.minus(vector: RealVector): RealVector = vector.map { toDouble() - it }
public operator fun RealVector.times(other: RealVector): RealVector =
mapIndexed { index, value -> value * other[index] }
public operator fun RealVector.times(number: Number): RealVector = map { it * number.toDouble() }
public operator fun Number.times(vector: RealVector): RealVector = vector * this
public operator fun RealVector.div(other: RealVector): RealVector =
mapIndexed { index, value -> value / other[index] }
public operator fun RealVector.div(number: Number): RealVector = map { it / number.toDouble() }
public operator fun Number.div(vector: RealVector): RealVector = vector.map { toDouble() / it }
//extended operations
public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) }
public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) }
public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) }
public fun sqrt(vector: RealVector): RealVector = vector.map { kotlin.math.sqrt(it) }
public fun RealVector.square(): RealVector = map { it.pow(2) }
public fun sin(vector: RealVector): RealVector = vector.map { kotlin.math.sin(it) }
public fun cos(vector: RealVector): RealVector = vector.map { kotlin.math.cos(it) }
public fun tan(vector: RealVector): RealVector = vector.map { kotlin.math.tan(it) }
public fun ln(vector: RealVector): RealVector = vector.map { kotlin.math.ln(it) }
public fun log10(vector: RealVector): RealVector = vector.map { kotlin.math.log10(it) }

31
kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt Normal file
View File

@ -0,0 +1,31 @@
package kscience.kmath.real
import kscience.kmath.linear.BufferMatrix
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.RealBuffer
/**
* Optimized dot product for real matrices
*/
public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
val resultArray = DoubleArray(this.rowNum * other.colNum)
//convert to array to insure there is no memory indirection
fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
this.array
else
DoubleArray(size) { get(it) }
val a = this.buffer.unsafeArray()
val b = other.buffer.unsafeArray()
for (i in (0 until rowNum))
for (j in (0 until other.colNum))
for (k in (0 until colNum))
resultArray[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
val buffer = RealBuffer(resultArray)
return BufferMatrix(rowNum, other.colNum, buffer)
}

15
kmath-core/src/commonMain/kotlin/kscience/kmath/misc/Grids.kt → kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/grids.kt
View File

@ -1,5 +1,6 @@
package kscience.kmath.misc
package kscience.kmath.real
import kscience.kmath.structures.asBuffer
import kotlin.math.abs
/**
@ -32,6 +33,9 @@ public fun ClosedFloatingPointRange<Double>.toSequenceWithStep(step: Double): Se
}
}
public infix fun ClosedFloatingPointRange<Double>.step(step: Double): RealVector =
toSequenceWithStep(step).toList().asBuffer()
/**
* Convert double range to sequence with the fixed number of points
*/
@ -39,12 +43,3 @@ public fun ClosedFloatingPointRange<Double>.toSequenceWithPoints(numPoints: Int)
require(numPoints > 1) { "The number of points should be more than 2" }
return toSequenceWithStep(abs(endInclusive - start) / (numPoints - 1))
}
/**
* Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints]
*/
@Deprecated("Replace by 'toSequenceWithPoints'")
public fun ClosedFloatingPointRange<Double>.toGrid(numPoints: Int): DoubleArray {
require(numPoints >= 2) { "Can't create generic grid with less than two points" }
return DoubleArray(numPoints) { i -> start + (endInclusive - start) / (numPoints - 1) * i }
}

8
kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/realBuffer.kt
View File

@ -1,8 +0,0 @@
package kscience.kmath.real
import kscience.kmath.structures.RealBuffer
/**
* Simplified [RealBuffer] to array comparison
*/
public fun RealBuffer.contentEquals(vararg doubles: Double): Boolean = array.contentEquals(doubles)

13
kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/GridTest.kt Normal file
View File

@ -0,0 +1,13 @@
package kaceince.kmath.real
import kscience.kmath.real.step
import kotlin.test.Test
import kotlin.test.assertEquals
class GridTest {
@Test
fun testStepGrid(){
val grid = 0.0..1.0 step 0.2
assertEquals(6, grid.size)
}
}

3
kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt → kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealMatrixTest.kt
View File

@ -1,9 +1,10 @@
package scientific.kmath.real
package kaceince.kmath.real
import kscience.kmath.linear.VirtualMatrix
import kscience.kmath.linear.build
import kscience.kmath.real.*
import kscience.kmath.structures.Matrix
import kscience.kmath.structures.contentEquals
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue

17
kmath-for-real/src/commonTest/kotlin/kscience/kmath/linear/VectorTest.kt → kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt
View File

@ -1,30 +1,33 @@
package kscience.kmath.linear
package kaceince.kmath.real
import kscience.kmath.linear.*
import kscience.kmath.operations.invoke
import kscience.kmath.real.RealVector
import kscience.kmath.real.plus
import kscience.kmath.structures.Buffer
import kotlin.test.Test
import kotlin.test.assertEquals
internal class VectorTest {
internal class RealVectorTest {
@Test
fun testSum() {
val vector1 = RealVector(5) { it.toDouble() }
val vector2 = RealVector(5) { 5 - it.toDouble() }
val vector1 = Buffer.real(5) { it.toDouble() }
val vector2 = Buffer.real(5) { 5 - it.toDouble() }
val sum = vector1 + vector2
assertEquals(5.0, sum[2])
}
@Test
fun testVectorToMatrix() {
val vector = RealVector(5) { it.toDouble() }
val vector = Buffer.real(5) { it.toDouble() }
val matrix = vector.asMatrix()
assertEquals(4.0, matrix[4, 0])
}
@Test
fun testDot() {
val vector1 = RealVector(5) { it.toDouble() }
val vector2 = RealVector(5) { 5 - it.toDouble() }
val vector1 = Buffer.real(5) { it.toDouble() }
val vector2 = Buffer.real(5) { 5 - it.toDouble() }
val matrix1 = vector1.asMatrix()
val matrix2 = vector2.asMatrix().transpose()
val product = MatrixContext.real { matrix1 dot matrix2 }

12
kmath-functions/src/commonMain/kotlin/kscience/kmath/functions/Polynomial.kt
View File

@ -8,13 +8,13 @@ import kotlin.contracts.contract
import kotlin.math.max
import kotlin.math.pow
// TODO make `inline`, when KT-41771 gets fixed
/**
* Polynomial coefficients without fixation on specific context they are applied to
* @param coefficients constant is the leftmost coefficient
*/
public inline class Polynomial<T : Any>(public val coefficients: List<T>)
@Suppress("FunctionName")
public fun <T : Any> Polynomial(vararg coefficients: T): Polynomial<T> = Polynomial(coefficients.toList())
public fun Polynomial<Double>.value(): Double = coefficients.reduceIndexed { index, acc, d -> acc + d.pow(index) }
@ -33,14 +33,6 @@ public fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring
res
}
/**
* Represent a polynomial as a context-dependent function
*/
public fun <T : Any, C : Ring<T>> Polynomial<T>.asMathFunction(): MathFunction<T, C, T> =
object : MathFunction<T, C, T> {
override fun C.invoke(arg: T): T = value(this, arg)
}
/**
* Represent the polynomial as a regular context-less function
*/
@ -49,7 +41,7 @@ public fun <T : Any, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T =
/**
* An algebra for polynomials
*/
public class PolynomialSpace<T : Any, C : Ring<T>>(public val ring: C) : Space<Polynomial<T>> {
public class PolynomialSpace<T : Any, C : Ring<T>>(private val ring: C) : Space<Polynomial<T>> {
public override val zero: Polynomial<T> = Polynomial(emptyList())
public override fun add(a: Polynomial<T>, b: Polynomial<T>): Polynomial<T> {

34
kmath-functions/src/commonMain/kotlin/kscience/kmath/functions/functions.kt
View File

@ -1,34 +0,0 @@
package kscience.kmath.functions
import kscience.kmath.operations.Algebra
import kscience.kmath.operations.RealField
// TODO make fun interface when KT-41770 is fixed
/**
* A regular function that could be called only inside specific algebra context
* @param T source type
* @param C source algebra constraint
* @param R result type
*/
public /*fun*/ interface MathFunction<T, C : Algebra<T>, R> {
public operator fun C.invoke(arg: T): R
}
public fun <R> MathFunction<Double, RealField, R>.invoke(arg: Double): R = RealField.invoke(arg)
/**
* A suspendable function defined in algebraic context
*/
// TODO make fun interface, when the new JVM IR is enabled
public interface SuspendableMathFunction<T, C : Algebra<T>, R> {
public suspend operator fun C.invoke(arg: T): R
}
public suspend fun <R> SuspendableMathFunction<Double, RealField, R>.invoke(arg: Double): R = RealField.invoke(arg)
/**
* A parametric function with parameter
*/
public fun interface ParametricFunction<T, P, C : Algebra<T>> {
public operator fun C.invoke(arg: T, parameter: P): T
}

5
kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
View File

@ -3,7 +3,6 @@ package kscience.kmath.histogram
import kscience.kmath.linear.Point
import kscience.kmath.operations.SpaceOperations
import kscience.kmath.operations.invoke
import kscience.kmath.real.asVector
import kscience.kmath.structures.*
import kotlin.math.floor
@ -123,8 +122,8 @@ public class RealHistogram(
*```
*/
public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram = RealHistogram(
ranges.map(ClosedFloatingPointRange<Double>::start).asVector(),
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asVector()
ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
)
/**

2
kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
View File

@ -4,6 +4,8 @@ import kscience.kmath.histogram.RealHistogram
import kscience.kmath.histogram.fill
import kscience.kmath.histogram.put
import kscience.kmath.real.RealVector
import kscience.kmath.real.invoke
import kscience.kmath.structures.Buffer
import kotlin.random.Random
import kotlin.test.*

4
kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
View File

@ -1,8 +1,8 @@
package kscience.kmath.histogram
import kscience.kmath.real.RealVector
import kscience.kmath.real.asVector
import kscience.kmath.structures.Buffer
import kscience.kmath.structures.asBuffer
import java.util.*
import kotlin.math.floor
@ -16,7 +16,7 @@ public class UnivariateBin(
//TODO add weighting
public override val value: Number get() = counter.sum()
public override val center: RealVector get() = doubleArrayOf(position).asVector()
public override val center: RealVector get() = doubleArrayOf(position).asBuffer()
public override val dimension: Int get() = 1
public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2)

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

@ -0,0 +1,9 @@
plugins {
id("ru.mipt.npm.jvm")
}
dependencies {
implementation("com.github.breandan:kaliningraph:0.1.4")
implementation("com.github.breandan:kotlingrad:0.4.0")
api(project(":kmath-ast"))
}

53
kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt Normal file
View File

@ -0,0 +1,53 @@
package kscience.kmath.kotlingrad
import edu.umontreal.kotlingrad.api.SFun
import kscience.kmath.ast.MST
import kscience.kmath.ast.MstAlgebra
import kscience.kmath.ast.MstExpression
import kscience.kmath.expressions.DifferentiableExpression
import kscience.kmath.expressions.Symbol
import kscience.kmath.operations.NumericAlgebra
/**
* Represents wrapper of [MstExpression] implementing [DifferentiableExpression].
*
* The principle of this API is converting the [mst] to an [SFun], differentiating it with Kotlin∇, then converting
* [SFun] back to [MST].
*
* @param T the type of number.
* @param A the [NumericAlgebra] of [T].
* @property expr the underlying [MstExpression].
*/
public inline class DifferentiableMstExpression<T, A>(public val expr: MstExpression<T, A>) :
DifferentiableExpression<T, MstExpression<T, A>> where A : NumericAlgebra<T>, T : Number {
public constructor(algebra: A, mst: MST) : this(MstExpression(algebra, mst))
/**
* The [MstExpression.algebra] of [expr].
*/
public val algebra: A
get() = expr.algebra
/**
* The [MstExpression.mst] of [expr].
*/
public val mst: MST
get() = expr.mst
public override fun invoke(arguments: Map<Symbol, T>): T = expr(arguments)
public override fun derivativeOrNull(symbols: List<Symbol>): MstExpression<T, A> = MstExpression(
algebra,
symbols.map(Symbol::identity)
.map(MstAlgebra::symbol)
.map { it.toSVar<KMathNumber<T, A>>() }
.fold(mst.toSFun(), SFun<KMathNumber<T, A>>::d)
.toMst(),
)
}
/**
* Wraps this [MstExpression] into [DifferentiableMstExpression].
*/
public fun <T : Number, A : NumericAlgebra<T>> MstExpression<T, A>.differentiable(): DifferentiableMstExpression<T, A> =
DifferentiableMstExpression(this)

18
kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/KMathNumber.kt Normal file
View File

@ -0,0 +1,18 @@
package kscience.kmath.kotlingrad
import edu.umontreal.kotlingrad.api.RealNumber
import edu.umontreal.kotlingrad.api.SConst
import kscience.kmath.operations.NumericAlgebra
/**
* Implements [RealNumber] by delegating its functionality to [NumericAlgebra].
*
* @param T the type of number.
* @param A the [NumericAlgebra] of [T].
* @property algebra the algebra.
* @param value the value of this number.
*/
public class KMathNumber<T, A>(public val algebra: A, value: T) :
RealNumber<KMathNumber<T, A>, T>(value) where T : Number, A : NumericAlgebra<T> {
public override fun wrap(number: Number): SConst<KMathNumber<T, A>> = SConst(algebra.number(number))
}

124
kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/ScalarsAdapters.kt Normal file
View File

@ -0,0 +1,124 @@
package kscience.kmath.kotlingrad
import edu.umontreal.kotlingrad.api.*
import kscience.kmath.ast.MST
import kscience.kmath.ast.MstAlgebra
import kscience.kmath.ast.MstExtendedField
import kscience.kmath.ast.MstExtendedField.unaryMinus
import kscience.kmath.operations.*
/**
* Maps [SVar] to [MST.Symbolic] directly.
*
* @receiver the variable.
* @return a node.
*/
public fun <X : SFun<X>> SVar<X>.toMst(): MST.Symbolic = MstAlgebra.symbol(name)
/**
* Maps [SVar] to [MST.Numeric] directly.
*
* @receiver the constant.
* @return a node.
*/
public fun <X : SFun<X>> SConst<X>.toMst(): MST.Numeric = MstAlgebra.number(doubleValue)
/**
* Maps [SFun] objects to [MST]. Some unsupported operations like [Derivative] are bound and converted then.
* [Power] operation is limited to constant right-hand side arguments.
*
* Detailed mapping is:
*
* - [SVar] -> [MstExtendedField.symbol];
* - [SConst] -> [MstExtendedField.number];
* - [Sum] -> [MstExtendedField.add];
* - [Prod] -> [MstExtendedField.multiply];
* - [Power] -> [MstExtendedField.power] (limited to constant exponents only);
* - [Negative] -> [MstExtendedField.unaryMinus];
* - [Log] -> [MstExtendedField.ln] (left) / [MstExtendedField.ln] (right);
* - [Sine] -> [MstExtendedField.sin];
* - [Cosine] -> [MstExtendedField.cos];
* - [Tangent] -> [MstExtendedField.tan];
* - [DProd] is vector operation, and it is requested to be evaluated;
* - [SComposition] is also requested to be evaluated eagerly;
* - [VSumAll] is requested to be evaluated;
* - [Derivative] is requested to be evaluated.
*
* @receiver the scalar function.
* @return a node.
*/
public fun <X : SFun<X>> SFun<X>.toMst(): MST = MstExtendedField {
when (this@toMst) {
is SVar -> toMst()
is SConst -> toMst()
is Sum -> left.toMst() + right.toMst()
is Prod -> left.toMst() * right.toMst()
is Power -> left.toMst() pow ((right as? SConst<*>)?.doubleValue ?: (right() as SConst<*>).doubleValue)
is Negative -> -input.toMst()
is Log -> ln(left.toMst()) / ln(right.toMst())
is Sine -> sin(input.toMst())
is Cosine -> cos(input.toMst())
is Tangent -> tan(input.toMst())
is DProd -> this@toMst().toMst()
is SComposition -> this@toMst().toMst()
is VSumAll<X, *> -> this@toMst().toMst()
is Derivative -> this@toMst().toMst()
}
}
/**
* Maps [MST.Numeric] to [SConst] directly.
*
* @receiver the node.
* @return a new constant.
*/
public fun <X : SFun<X>> MST.Numeric.toSConst(): SConst<X> = SConst(value)
/**
* Maps [MST.Symbolic] to [SVar] directly.
*
* @receiver the node.
* @param proto the prototype instance.
* @return a new variable.
*/
internal fun <X : SFun<X>> MST.Symbolic.toSVar(): SVar<X> = SVar(value)
/**
* Maps [MST] objects to [SFun]. Unsupported operations throw [IllegalStateException].
*
* Detailed mapping is:
*
* - [MST.Numeric] -> [SConst];
* - [MST.Symbolic] -> [SVar];
* - [MST.Unary] -> [Negative], [Sine], [Cosine], [Tangent], [Power], [Log];
* - [MST.Binary] -> [Sum], [Prod], [Power].
*
* @receiver the node.
* @param proto the prototype instance.
* @return a scalar function.
*/
public fun <X : SFun<X>> MST.toSFun(): SFun<X> = when (this) {
is MST.Numeric -> toSConst()
is MST.Symbolic -> toSVar()
is MST.Unary -> when (operation) {
SpaceOperations.PLUS_OPERATION -> +value.toSFun<X>()
SpaceOperations.MINUS_OPERATION -> -value.toSFun<X>()
TrigonometricOperations.SIN_OPERATION -> sin(value.toSFun())
TrigonometricOperations.COS_OPERATION -> cos(value.toSFun())
TrigonometricOperations.TAN_OPERATION -> tan(value.toSFun())
PowerOperations.SQRT_OPERATION -> sqrt(value.toSFun())
ExponentialOperations.EXP_OPERATION -> exp(value.toSFun())
ExponentialOperations.LN_OPERATION -> value.toSFun<X>().ln()
else -> error("Unary operation $operation not defined in $this")
}
is MST.Binary -> when (operation) {
SpaceOperations.PLUS_OPERATION -> left.toSFun<X>() + right.toSFun()
SpaceOperations.MINUS_OPERATION -> left.toSFun<X>() - right.toSFun()
RingOperations.TIMES_OPERATION -> left.toSFun<X>() * right.toSFun()
FieldOperations.DIV_OPERATION -> left.toSFun<X>() / right.toSFun()
PowerOperations.POW_OPERATION -> left.toSFun<X>() pow (right as MST.Numeric).toSConst()
else -> error("Binary operation $operation not defined in $this")
}
}

64
kmath-kotlingrad/src/test/kotlin/kscience/kmath/kotlingrad/AdaptingTests.kt Normal file
View File

@ -0,0 +1,64 @@
package kscience.kmath.kotlingrad
import edu.umontreal.kotlingrad.api.*
import kscience.kmath.asm.compile
import kscience.kmath.ast.MstAlgebra
import kscience.kmath.ast.MstExpression
import kscience.kmath.ast.parseMath
import kscience.kmath.expressions.invoke
import kscience.kmath.operations.RealField
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
import kotlin.test.fail
internal class AdaptingTests {
@Test
fun symbol() {
val c1 = MstAlgebra.symbol("x")
assertTrue(c1.toSVar<KMathNumber<Double, RealField>>().name == "x")
val c2 = "kitten".parseMath().toSFun<KMathNumber<Double, RealField>>()
if (c2 is SVar) assertTrue(c2.name == "kitten") else fail()
}
@Test
fun number() {
val c1 = MstAlgebra.number(12354324)
assertTrue(c1.toSConst<DReal>().doubleValue == 12354324.0)
val c2 = "0.234".parseMath().toSFun<KMathNumber<Double, RealField>>()
if (c2 is SConst) assertTrue(c2.doubleValue == 0.234) else fail()
val c3 = "1e-3".parseMath().toSFun<KMathNumber<Double, RealField>>()
if (c3 is SConst) assertEquals(0.001, c3.value) else fail()
}
@Test
fun simpleFunctionShape() {
val linear = "2*x+16".parseMath().toSFun<KMathNumber<Double, RealField>>()
if (linear !is Sum) fail()
if (linear.left !is Prod) fail()
if (linear.right !is SConst) fail()
}
@Test
fun simpleFunctionDerivative() {
val x = MstAlgebra.symbol("x").toSVar<KMathNumber<Double, RealField>>()
val quadratic = "x^2-4*x-44".parseMath().toSFun<KMathNumber<Double, RealField>>()
val actualDerivative = MstExpression(RealField, quadratic.d(x).toMst()).compile()
val expectedDerivative = MstExpression(RealField, "2*x-4".parseMath()).compile()
assertEquals(actualDerivative("x" to 123.0), expectedDerivative("x" to 123.0))
}
@Test
fun moreComplexDerivative() {
val x = MstAlgebra.symbol("x").toSVar<KMathNumber<Double, RealField>>()
val composition = "-sqrt(sin(x^2)-cos(x)^2-16*x)".parseMath().toSFun<KMathNumber<Double, RealField>>()
val actualDerivative = MstExpression(RealField, composition.d(x).toMst()).compile()
val expectedDerivative = MstExpression(
RealField,
"-(2*x*cos(x^2)+2*sin(x)*cos(x)-16)/(2*sqrt(sin(x^2)-16*x-cos(x)^2))".parseMath()
).compile()
assertEquals(actualDerivative("x" to 0.1), expectedDerivative("x" to 0.1))
}
}

82
kmath-nd4j/README.md Normal file
View File

@ -0,0 +1,82 @@
# ND4J NDStructure implementation (`kmath-nd4j`)
This subproject implements the following features:
- [nd4jarraystructure](src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
- [nd4jarrayrings](src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Rings over Nd4jArrayStructure of Int and Long
- [nd4jarrayfields](src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : Fields over Nd4jArrayStructure of Float and Double
> #### Artifact:
>
> This module artifact: `kscience.kmath:kmath-nd4j:0.2.0-dev-4`.
>
> Bintray release version: [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-nd4j/_latestVersion)
>
> Bintray development version: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-nd4j/_latestVersion)
>
> **Gradle:**
>
> ```gradle
> repositories {
> maven { url "https://dl.bintray.com/kotlin/kotlin-eap" }
> maven { url 'https://dl.bintray.com/mipt-npm/kscience' }
> maven { url 'https://dl.bintray.com/mipt-npm/dev' }
> maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
> }
>
> dependencies {
> implementation 'kscience.kmath:kmath-nd4j:0.2.0-dev-4'
> }
> ```
> **Gradle Kotlin DSL:**
>
> ```kotlin
> repositories {
> maven("https://dl.bintray.com/kotlin/kotlin-eap")
> maven("https://dl.bintray.com/mipt-npm/kscience")
> maven("https://dl.bintray.com/mipt-npm/dev")
> maven("https://dl.bintray.com/hotkeytlt/maven")
> }
>
> dependencies {
> implementation("kscience.kmath:kmath-nd4j:0.2.0-dev-4")
> }
> ```
## Examples
NDStructure wrapper for INDArray:
```kotlin
import org.nd4j.linalg.factory.*
import scientifik.kmath.nd4j.*
import scientifik.kmath.structures.*
val array = Nd4j.ones(2, 2).asRealStructure()
println(array[0, 0]) // 1.0
array[intArrayOf(0, 0)] = 24.0
println(array[0, 0]) // 24.0
```
Fast element-wise and in-place arithmetics for INDArray:
```kotlin
import org.nd4j.linalg.factory.*
import scientifik.kmath.nd4j.*
import scientifik.kmath.operations.*
val field = RealNd4jArrayField(intArrayOf(2, 2))
val array = Nd4j.rand(2, 2).asRealStructure()
val res = field {
(25.0 / array + 20) * 4
}
println(res.ndArray)
// [[ 250.6449, 428.5840],
// [ 269.7913, 202.2077]]
```
Contributed by [Iaroslav Postovalov](https://github.com/CommanderTvis).

37
kmath-nd4j/build.gradle.kts Normal file
View File

@ -0,0 +1,37 @@
import ru.mipt.npm.gradle.Maturity
plugins {
id("ru.mipt.npm.jvm")
}
dependencies {
api(project(":kmath-core"))
api("org.nd4j:nd4j-api:1.0.0-beta7")
testImplementation("org.deeplearning4j:deeplearning4j-core:1.0.0-beta7")
testImplementation("org.nd4j:nd4j-native-platform:1.0.0-beta7")
testImplementation("org.slf4j:slf4j-simple:1.7.30")
}
readme {
description = "ND4J NDStructure implementation and according NDAlgebra classes"
maturity = Maturity.EXPERIMENTAL
propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
feature(
id = "nd4jarraystructure",
description = "NDStructure wrapper for INDArray",
ref = "src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt"
)
feature(
id = "nd4jarrayrings",
description = "Rings over Nd4jArrayStructure of Int and Long",
ref = "src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt"
)
feature(
id = "nd4jarrayfields",
description = "Fields over Nd4jArrayStructure of Float and Double",
ref = "src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt"
)
}

43
kmath-nd4j/docs/README-TEMPLATE.md Normal file
View File

@ -0,0 +1,43 @@
# ND4J NDStructure implementation (`kmath-nd4j`)
This subproject implements the following features:
${features}
${artifact}
## Examples
NDStructure wrapper for INDArray:
```kotlin
import org.nd4j.linalg.factory.*
import scientifik.kmath.nd4j.*
import scientifik.kmath.structures.*
val array = Nd4j.ones(2, 2).asRealStructure()
println(array[0, 0]) // 1.0
array[intArrayOf(0, 0)] = 24.0
println(array[0, 0]) // 24.0
```
Fast element-wise and in-place arithmetics for INDArray:
```kotlin
import org.nd4j.linalg.factory.*
import scientifik.kmath.nd4j.*
import scientifik.kmath.operations.*
val field = RealNd4jArrayField(intArrayOf(2, 2))
val array = Nd4j.rand(2, 2).asRealStructure()
val res = field {
(25.0 / array + 20) * 4
}
println(res.ndArray)
// [[ 250.6449, 428.5840],
// [ 269.7913, 202.2077]]
```
Contributed by [Iaroslav Postovalov](https://github.com/CommanderTvis).

349
kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayAlgebra.kt Normal file
View File

@ -0,0 +1,349 @@
package kscience.kmath.nd4j
import kscience.kmath.operations.*
import kscience.kmath.structures.NDAlgebra
import kscience.kmath.structures.NDField
import kscience.kmath.structures.NDRing
import kscience.kmath.structures.NDSpace
import org.nd4j.linalg.api.ndarray.INDArray
import org.nd4j.linalg.factory.Nd4j
/**
* Represents [NDAlgebra] over [Nd4jArrayAlgebra].
*
* @param T the type of ND-structure element.
* @param C the type of the element context.
*/
public interface Nd4jArrayAlgebra<T, C> : NDAlgebra<T, C, Nd4jArrayStructure<T>> {
/**
* Wraps [INDArray] to [N].
*/
public fun INDArray.wrap(): Nd4jArrayStructure<T>
public override fun produce(initializer: C.(IntArray) -> T): Nd4jArrayStructure<T> {
val struct = Nd4j.create(*shape)!!.wrap()
struct.indicesIterator().forEach { struct[it] = elementContext.initializer(it) }
return struct
}
public override fun map(arg: Nd4jArrayStructure<T>, transform: C.(T) -> T): Nd4jArrayStructure<T> {
check(arg)
val newStruct = arg.ndArray.dup().wrap()
newStruct.elements().forEach { (idx, value) -> newStruct[idx] = elementContext.transform(value) }
return newStruct
}
public override fun mapIndexed(
arg: Nd4jArrayStructure<T>,
transform: C.(index: IntArray, T) -> T
): Nd4jArrayStructure<T> {
check(arg)
val new = Nd4j.create(*shape).wrap()
new.indicesIterator().forEach { idx -> new[idx] = elementContext.transform(idx, arg[idx]) }
return new
}
public override fun combine(
a: Nd4jArrayStructure<T>,
b: Nd4jArrayStructure<T>,
transform: C.(T, T) -> T
): Nd4jArrayStructure<T> {
check(a, b)
val new = Nd4j.create(*shape).wrap()
new.indicesIterator().forEach { idx -> new[idx] = elementContext.transform(a[idx], b[idx]) }
return new
}
}
/**
* Represents [NDSpace] over [Nd4jArrayStructure].
*
* @param T the type of the element contained in ND structure.
* @param S the type of space of structure elements.
*/
public interface Nd4jArraySpace<T, S> : NDSpace<T, S, Nd4jArrayStructure<T>>,
Nd4jArrayAlgebra<T, S> where S : Space<T> {
public override val zero: Nd4jArrayStructure<T>
get() = Nd4j.zeros(*shape).wrap()
public override fun add(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
check(a, b)
return a.ndArray.add(b.ndArray).wrap()
}
public override operator fun Nd4jArrayStructure<T>.minus(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
check(this, b)
return ndArray.sub(b.ndArray).wrap()
}
public override operator fun Nd4jArrayStructure<T>.unaryMinus(): Nd4jArrayStructure<T> {
check(this)
return ndArray.neg().wrap()
}
public override fun multiply(a: Nd4jArrayStructure<T>, k: Number): Nd4jArrayStructure<T> {
check(a)
return a.ndArray.mul(k).wrap()
}
public override operator fun Nd4jArrayStructure<T>.div(k: Number): Nd4jArrayStructure<T> {
check(this)
return ndArray.div(k).wrap()
}
public override operator fun Nd4jArrayStructure<T>.times(k: Number): Nd4jArrayStructure<T> {
check(this)
return ndArray.mul(k).wrap()
}
}
/**
* Represents [NDRing] over [Nd4jArrayStructure].
*
* @param T the type of the element contained in ND structure.
* @param R the type of ring of structure elements.
*/
public interface Nd4jArrayRing<T, R> : NDRing<T, R, Nd4jArrayStructure<T>>, Nd4jArraySpace<T, R> where R : Ring<T> {
public override val one: Nd4jArrayStructure<T>
get() = Nd4j.ones(*shape).wrap()
public override fun multiply(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
check(a, b)
return a.ndArray.mul(b.ndArray).wrap()
}
public override operator fun Nd4jArrayStructure<T>.minus(b: Number): Nd4jArrayStructure<T> {
check(this)
return ndArray.sub(b).wrap()
}
public override operator fun Nd4jArrayStructure<T>.plus(b: Number): Nd4jArrayStructure<T> {
check(this)
return ndArray.add(b).wrap()
}
public override operator fun Number.minus(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
check(b)
return b.ndArray.rsub(this).wrap()
}
public companion object {
private val intNd4jArrayRingCache: ThreadLocal<MutableMap<IntArray, IntNd4jArrayRing>> =
ThreadLocal.withInitial { hashMapOf() }
private val longNd4jArrayRingCache: ThreadLocal<MutableMap<IntArray, LongNd4jArrayRing>> =
ThreadLocal.withInitial { hashMapOf() }
/**
* Creates an [NDRing] for [Int] values or pull it from cache if it was created previously.
*/
public fun int(vararg shape: Int): Nd4jArrayRing<Int, IntRing> =
intNd4jArrayRingCache.get().getOrPut(shape) { IntNd4jArrayRing(shape) }
/**
* Creates an [NDRing] for [Long] values or pull it from cache if it was created previously.
*/
public fun long(vararg shape: Int): Nd4jArrayRing<Long, LongRing> =
longNd4jArrayRingCache.get().getOrPut(shape) { LongNd4jArrayRing(shape) }
/**
* Creates a most suitable implementation of [NDRing] using reified class.
*/
@Suppress("UNCHECKED_CAST")
public inline fun <reified T : Any> auto(vararg shape: Int): Nd4jArrayRing<T, out Ring<T>> = when {
T::class == Int::class -> int(*shape) as Nd4jArrayRing<T, out Ring<T>>
T::class == Long::class -> long(*shape) as Nd4jArrayRing<T, out Ring<T>>
else -> throw UnsupportedOperationException("This factory method only supports Int and Long types.")
}
}
}
/**
* Represents [NDField] over [Nd4jArrayStructure].
*
* @param T the type of the element contained in ND structure.
* @param N the type of ND structure.
* @param F the type field of structure elements.
*/
public interface Nd4jArrayField<T, F> : NDField<T, F, Nd4jArrayStructure<T>>, Nd4jArrayRing<T, F> where F : Field<T> {
public override fun divide(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
check(a, b)
return a.ndArray.div(b.ndArray).wrap()
}
public override operator fun Number.div(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
check(b)
return b.ndArray.rdiv(this).wrap()
}
public companion object {
private val floatNd4jArrayFieldCache: ThreadLocal<MutableMap<IntArray, FloatNd4jArrayField>> =
ThreadLocal.withInitial { hashMapOf() }
private val realNd4jArrayFieldCache: ThreadLocal<MutableMap<IntArray, RealNd4jArrayField>> =
ThreadLocal.withInitial { hashMapOf() }
/**
* Creates an [NDField] for [Float] values or pull it from cache if it was created previously.
*/
public fun float(vararg shape: Int): Nd4jArrayRing<Float, FloatField> =
floatNd4jArrayFieldCache.get().getOrPut(shape) { FloatNd4jArrayField(shape) }
/**
* Creates an [NDField] for [Double] values or pull it from cache if it was created previously.
*/
public fun real(vararg shape: Int): Nd4jArrayRing<Double, RealField> =
realNd4jArrayFieldCache.get().getOrPut(shape) { RealNd4jArrayField(shape) }
/**
* Creates a most suitable implementation of [NDRing] using reified class.
*/
@Suppress("UNCHECKED_CAST")
public inline fun <reified T : Any> auto(vararg shape: Int): Nd4jArrayField<T, out Field<T>> = when {
T::class == Float::class -> float(*shape) as Nd4jArrayField<T, out Field<T>>
T::class == Double::class -> real(*shape) as Nd4jArrayField<T, out Field<T>>
else -> throw UnsupportedOperationException("This factory method only supports Float and Double types.")
}
}
}
/**
* Represents [NDField] over [Nd4jArrayRealStructure].
*/
public class RealNd4jArrayField(public override val shape: IntArray) : Nd4jArrayField<Double, RealField> {
public override val elementContext: RealField
get() = RealField
public override fun INDArray.wrap(): Nd4jArrayStructure<Double> = check(asRealStructure())
public override operator fun Nd4jArrayStructure<Double>.div(arg: Double): Nd4jArrayStructure<Double> {
check(this)
return ndArray.div(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Double>.plus(arg: Double): Nd4jArrayStructure<Double> {
check(this)
return ndArray.add(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Double>.minus(arg: Double): Nd4jArrayStructure<Double> {
check(this)
return ndArray.sub(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Double>.times(arg: Double): Nd4jArrayStructure<Double> {
check(this)
return ndArray.mul(arg).wrap()
}
public override operator fun Double.div(arg: Nd4jArrayStructure<Double>): Nd4jArrayStructure<Double> {
check(arg)
return arg.ndArray.rdiv(this).wrap()
}
public override operator fun Double.minus(arg: Nd4jArrayStructure<Double>): Nd4jArrayStructure<Double> {
check(arg)
return arg.ndArray.rsub(this).wrap()
}
}
/**
* Represents [NDField] over [Nd4jArrayStructure] of [Float].
*/
public class FloatNd4jArrayField(public override val shape: IntArray) : Nd4jArrayField<Float, FloatField> {
public override val elementContext: FloatField
get() = FloatField
public override fun INDArray.wrap(): Nd4jArrayStructure<Float> = check(asFloatStructure())
public override operator fun Nd4jArrayStructure<Float>.div(arg: Float): Nd4jArrayStructure<Float> {
check(this)
return ndArray.div(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Float>.plus(arg: Float): Nd4jArrayStructure<Float> {
check(this)
return ndArray.add(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Float>.minus(arg: Float): Nd4jArrayStructure<Float> {
check(this)
return ndArray.sub(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Float>.times(arg: Float): Nd4jArrayStructure<Float> {
check(this)
return ndArray.mul(arg).wrap()
}
public override operator fun Float.div(arg: Nd4jArrayStructure<Float>): Nd4jArrayStructure<Float> {
check(arg)
return arg.ndArray.rdiv(this).wrap()
}
public override operator fun Float.minus(arg: Nd4jArrayStructure<Float>): Nd4jArrayStructure<Float> {
check(arg)
return arg.ndArray.rsub(this).wrap()
}
}
/**
* Represents [NDRing] over [Nd4jArrayIntStructure].
*/
public class IntNd4jArrayRing(public override val shape: IntArray) : Nd4jArrayRing<Int, IntRing> {
public override val elementContext: IntRing
get() = IntRing
public override fun INDArray.wrap(): Nd4jArrayStructure<Int> = check(asIntStructure())
public override operator fun Nd4jArrayStructure<Int>.plus(arg: Int): Nd4jArrayStructure<Int> {
check(this)
return ndArray.add(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Int>.minus(arg: Int): Nd4jArrayStructure<Int> {
check(this)
return ndArray.sub(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Int>.times(arg: Int): Nd4jArrayStructure<Int> {
check(this)
return ndArray.mul(arg).wrap()
}
public override operator fun Int.minus(arg: Nd4jArrayStructure<Int>): Nd4jArrayStructure<Int> {
check(arg)
return arg.ndArray.rsub(this).wrap()
}
}
/**
* Represents [NDRing] over [Nd4jArrayStructure] of [Long].
*/
public class LongNd4jArrayRing(public override val shape: IntArray) : Nd4jArrayRing<Long, LongRing> {
public override val elementContext: LongRing
get() = LongRing
public override fun INDArray.wrap(): Nd4jArrayStructure<Long> = check(asLongStructure())
public override operator fun Nd4jArrayStructure<Long>.plus(arg: Long): Nd4jArrayStructure<Long> {
check(this)
return ndArray.add(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Long>.minus(arg: Long): Nd4jArrayStructure<Long> {
check(this)
return ndArray.sub(arg).wrap()
}
public override operator fun Nd4jArrayStructure<Long>.times(arg: Long): Nd4jArrayStructure<Long> {
check(this)
return ndArray.mul(arg).wrap()
}
public override operator fun Long.minus(arg: Nd4jArrayStructure<Long>): Nd4jArrayStructure<Long> {
check(arg)
return arg.ndArray.rsub(this).wrap()
}
}

62
kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayIterator.kt Normal file
View File

@ -0,0 +1,62 @@
package kscience.kmath.nd4j
import org.nd4j.linalg.api.ndarray.INDArray
import org.nd4j.linalg.api.shape.Shape
private class Nd4jArrayIndicesIterator(private val iterateOver: INDArray) : Iterator<IntArray> {
private var i: Int = 0
override fun hasNext(): Boolean = i < iterateOver.length()
override fun next(): IntArray {
val la = if (iterateOver.ordering() == 'c')
Shape.ind2subC(iterateOver, i++.toLong())!!
else
Shape.ind2sub(iterateOver, i++.toLong())!!
return la.toIntArray()
}
}
internal fun INDArray.indicesIterator(): Iterator<IntArray> = Nd4jArrayIndicesIterator(this)
private sealed class Nd4jArrayIteratorBase<T>(protected val iterateOver: INDArray) : Iterator<Pair<IntArray, T>> {
private var i: Int = 0
final override fun hasNext(): Boolean = i < iterateOver.length()
abstract fun getSingle(indices: LongArray): T
final override fun next(): Pair<IntArray, T> {
val la = if (iterateOver.ordering() == 'c')
Shape.ind2subC(iterateOver, i++.toLong())!!
else
Shape.ind2sub(iterateOver, i++.toLong())!!
return la.toIntArray() to getSingle(la)
}
}
private class Nd4jArrayRealIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Double>(iterateOver) {
override fun getSingle(indices: LongArray): Double = iterateOver.getDouble(*indices)
}
internal fun INDArray.realIterator(): Iterator<Pair<IntArray, Double>> = Nd4jArrayRealIterator(this)
private class Nd4jArrayLongIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Long>(iterateOver) {
override fun getSingle(indices: LongArray) = iterateOver.getLong(*indices)
}
internal fun INDArray.longIterator(): Iterator<Pair<IntArray, Long>> = Nd4jArrayLongIterator(this)
private class Nd4jArrayIntIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Int>(iterateOver) {
override fun getSingle(indices: LongArray) = iterateOver.getInt(*indices.toIntArray())
}
internal fun INDArray.intIterator(): Iterator<Pair<IntArray, Int>> = Nd4jArrayIntIterator(this)
private class Nd4jArrayFloatIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Float>(iterateOver) {
override fun getSingle(indices: LongArray) = iterateOver.getFloat(*indices)
}
internal fun INDArray.floatIterator(): Iterator<Pair<IntArray, Float>> = Nd4jArrayFloatIterator(this)

68
kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayStructure.kt Normal file
View File

@ -0,0 +1,68 @@
package kscience.kmath.nd4j
import kscience.kmath.structures.MutableNDStructure
import kscience.kmath.structures.NDStructure
import org.nd4j.linalg.api.ndarray.INDArray
/**
* Represents a [NDStructure] wrapping an [INDArray] object.
*
* @param T the type of items.
*/
public sealed class Nd4jArrayStructure<T> : MutableNDStructure<T> {
/**
* The wrapped [INDArray].
*/
public abstract val ndArray: INDArray
public override val shape: IntArray
get() = ndArray.shape().toIntArray()
internal abstract fun elementsIterator(): Iterator<Pair<IntArray, T>>
internal fun indicesIterator(): Iterator<IntArray> = ndArray.indicesIterator()
public override fun elements(): Sequence<Pair<IntArray, T>> = Sequence(::elementsIterator)
}
private data class Nd4jArrayIntStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Int>() {
override fun elementsIterator(): Iterator<Pair<IntArray, Int>> = ndArray.intIterator()
override fun get(index: IntArray): Int = ndArray.getInt(*index)
override fun set(index: IntArray, value: Int): Unit = run { ndArray.putScalar(index, value) }
}
/**
* Wraps this [INDArray] to [Nd4jArrayStructure].
*/
public fun INDArray.asIntStructure(): Nd4jArrayStructure<Int> = Nd4jArrayIntStructure(this)
private data class Nd4jArrayLongStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Long>() {
override fun elementsIterator(): Iterator<Pair<IntArray, Long>> = ndArray.longIterator()
override fun get(index: IntArray): Long = ndArray.getLong(*index.toLongArray())
override fun set(index: IntArray, value: Long): Unit = run { ndArray.putScalar(index, value.toDouble()) }
}
/**
* Wraps this [INDArray] to [Nd4jArrayStructure].
*/
public fun INDArray.asLongStructure(): Nd4jArrayStructure<Long> = Nd4jArrayLongStructure(this)
private data class Nd4jArrayRealStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Double>() {
override fun elementsIterator(): Iterator<Pair<IntArray, Double>> = ndArray.realIterator()
override fun get(index: IntArray): Double = ndArray.getDouble(*index)
override fun set(index: IntArray, value: Double): Unit = run { ndArray.putScalar(index, value) }
}
/**
* Wraps this [INDArray] to [Nd4jArrayStructure].
*/
public fun INDArray.asRealStructure(): Nd4jArrayStructure<Double> = Nd4jArrayRealStructure(this)
private data class Nd4jArrayFloatStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Float>() {
override fun elementsIterator(): Iterator<Pair<IntArray, Float>> = ndArray.floatIterator()
override fun get(index: IntArray): Float = ndArray.getFloat(*index)
override fun set(index: IntArray, value: Float): Unit = run { ndArray.putScalar(index, value) }
}
/**
* Wraps this [INDArray] to [Nd4jArrayStructure].
*/
public fun INDArray.asFloatStructure(): Nd4jArrayStructure<Float> = Nd4jArrayFloatStructure(this)

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