Fix benchmark results script

Cosmetic change Double -> Float64

.github/CODEOWNERS

@@ -0,0 +1,3 @@
+@altavir
+
+/kmath-trajectory @ESchouten

.github/workflows/build.yml

@@ -0,0 +1,24 @@
+name: Gradle build
+
+on:
+  push:
+    branches: [ dev, master ]
+  pull_request:
+
+jobs:
+  build:
+    runs-on: windows-latest
+    timeout-minutes: 20
+    steps:
+      - uses: actions/checkout@v3
+      - uses: actions/setup-java@v3.5.1
+        with:
+          java-version: '11'
+          distribution: 'liberica'
+          cache: 'gradle'
+      - name: Gradle Wrapper Validation
+        uses: gradle/wrapper-validation-action@v1.0.4
+      - name: Gradle Build
+        uses: gradle/gradle-build-action@v2.4.2
+        with:
+          arguments: test jvmTest

.github/workflows/gradle.yml

@@ -1,17 +0,0 @@
-name: Gradle build
-
-on: [push]
-
-jobs:
-  build:
-
-    runs-on: ubuntu-latest
-
-    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

.github/workflows/pages.yml

@@ -0,0 +1,31 @@
+name: Dokka publication
+
+on:
+  workflow_dispatch:
+  release:
+    types: [ created ]
+
+jobs:
+  build:
+    runs-on: ubuntu-20.04
+    timeout-minutes: 40
+    steps:
+      - uses: actions/checkout@v3.0.0
+      - uses: actions/setup-java@v3.0.0
+        with:
+          java-version: 11
+          distribution: liberica
+      - name: Cache konan
+        uses: actions/cache@v3.0.1
+        with:
+          path: ~/.konan
+          key: ${{ runner.os }}-gradle-${{ hashFiles('*.gradle.kts') }}
+          restore-keys: |
+            ${{ runner.os }}-gradle-
+      - uses: gradle/gradle-build-action@v2.4.2
+        with:
+          arguments: dokkaHtmlMultiModule --no-parallel
+      - uses: JamesIves/github-pages-deploy-action@v4.3.0
+        with:
+          branch: gh-pages
+          folder: build/dokka/htmlMultiModule

.github/workflows/publish.yml

@@ -0,0 +1,50 @@
+name: Gradle publish
+
+on:
+  workflow_dispatch:
+  release:
+    types: [ created ]
+
+jobs:
+  publish:
+    environment:
+      name: publish
+    strategy:
+      matrix:
+        os: [ macOS-latest, windows-latest ]
+    runs-on: ${{matrix.os}}
+    steps:
+      - uses: actions/checkout@v3.0.0
+      - uses: actions/setup-java@v3.10.0
+        with:
+          java-version: 11
+          distribution: liberica
+      - name: Cache konan
+        uses: actions/cache@v3.0.1
+        with:
+          path: ~/.konan
+          key: ${{ runner.os }}-gradle-${{ hashFiles('*.gradle.kts') }}
+          restore-keys: |
+            ${{ runner.os }}-gradle-
+      - name: Publish Windows Artifacts
+        if: matrix.os == 'windows-latest'
+        uses: gradle/gradle-build-action@v2.4.2
+        with:
+          arguments: |
+            publishAllPublicationsToSpaceRepository
+            -Ppublishing.targets=all
+            -Ppublishing.space.user=${{ secrets.SPACE_APP_ID }}
+            -Ppublishing.space.token=${{ secrets.SPACE_APP_SECRET }}
+      - name: Publish Mac Artifacts
+        if: matrix.os == 'macOS-latest'
+        uses: gradle/gradle-build-action@v2.4.2
+        with:
+          arguments: |
+            publishMacosX64PublicationToSpaceRepository
+            publishMacosArm64PublicationToSpaceRepository
+            publishIosX64PublicationToSpaceRepository
+            publishIosArm64PublicationToSpaceRepository
+            publishIosSimulatorArm64PublicationToSpaceRepository
+            -Ppublishing.targets=all
+            -Ppublishing.space.user=${{ secrets.SPACE_APP_ID }}
+            -Ppublishing.space.token=${{ secrets.SPACE_APP_SECRET }}

.gitignore

@@ -1,7 +1,12 @@
 .gradle
 build/
 out/
+
 .idea/
+.vscode/
+.fleet/
+.kotlin/
+
 
 # Avoid ignoring Gradle wrapper jar file (.jar files are usually ignored)
 !gradle-wrapper.jar
@@ -9,4 +14,11 @@ out/
 # Cache of project
 .gradletasknamecache
 
-gradle.properties
+# Generated by javac -h and runtime
+*.class
+*.log
+
+!/.idea/copyright/
+!/.idea/scopes/
+/gradle/yarn.lock
+

.idea/copyright/kmath.xml

@@ -0,0 +1,7 @@
+<component name="CopyrightManager">
+  <copyright>
+    <option name="allowReplaceRegexp" value="Copyright \d{4}-\d{4} KMath" />
+    <option name="notice" value="Copyright 2018-&amp;#36;today.year KMath contributors.&#10;Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file." />
+    <option name="myName" value="kmath" />
+  </copyright>
+</component>

.idea/copyright/profiles_settings.xml

@@ -0,0 +1,21 @@
+<component name="CopyrightManager">
+  <settings>
+    <module2copyright>
+      <element module="Apply copyright" copyright="kmath" />
+    </module2copyright>
+    <LanguageOptions name="Groovy">
+      <option name="fileTypeOverride" value="1" />
+    </LanguageOptions>
+    <LanguageOptions name="HTML">
+      <option name="fileTypeOverride" value="1" />
+      <option name="prefixLines" value="false" />
+    </LanguageOptions>
+    <LanguageOptions name="Properties">
+      <option name="fileTypeOverride" value="1" />
+    </LanguageOptions>
+    <LanguageOptions name="XML">
+      <option name="fileTypeOverride" value="1" />
+      <option name="prefixLines" value="false" />
+    </LanguageOptions>
+  </settings>
+</component>

.idea/scopes/Apply_copyright.xml

@@ -0,0 +1,3 @@
+<component name="DependencyValidationManager">
+    <scope name="Apply copyright" pattern="!file[*]:*//testData//*&amp;&amp;!file[*]:testData//*&amp;&amp;!file[*]:*.gradle.kts&amp;&amp;!file[*]:*.gradle&amp;&amp;!file[group:kotlin-ultimate]:*/&amp;&amp;!file[kotlin.libraries]:stdlib/api//*" />
+</component>

.space.kts

@@ -0,0 +1,48 @@
+import kotlin.io.path.readText
+
+val projectName = "kmath"
+
+job("Build") {
+    //Perform only jvm tests
+    gradlew("spc.registry.jetbrains.space/p/sci/containers/kotlin-ci:1.0.3", "test", "jvmTest")
+}
+
+job("Publish") {
+    startOn {
+        gitPush { enabled = false }
+    }
+    container("spc.registry.jetbrains.space/p/sci/containers/kotlin-ci:1.0.3") {
+        env["SPACE_USER"] = "{{ project:space_user }}"
+        env["SPACE_TOKEN"] = "{{ project:space_token }}"
+        kotlinScript { api ->
+
+            val spaceUser = System.getenv("SPACE_USER")
+            val spaceToken = System.getenv("SPACE_TOKEN")
+
+            // write the version to the build directory
+            api.gradlew("version")
+
+            //read the version from build file
+            val version = java.nio.file.Path.of("build/project-version.txt").readText()
+
+            val revisionSuffix = if (version.endsWith("SNAPSHOT")) {
+                "-" + api.gitRevision().take(7)
+            } else {
+                ""
+            }
+
+            api.space().projects.automation.deployments.start(
+                project = api.projectIdentifier(),
+                targetIdentifier = TargetIdentifier.Key(projectName),
+                version = version + revisionSuffix,
+                // automatically update deployment status based on the status of a job
+                syncWithAutomationJob = true
+            )
+            api.gradlew(
+                "publishAllPublicationsToSpaceRepository",
+                "-Ppublishing.space.user=\"$spaceUser\"",
+                "-Ppublishing.space.token=\"$spaceToken\"",
+            )
+        }
+    }
+}

CHANGELOG.md

@@ -1,17 +1,248 @@
 # KMath
 
-## [Unreleased]
+## Unreleased
 
 ### Added
+- Metropolis-Hastings sampler
+
+### Changed
+
+### Deprecated
+
+### Removed
+
+### Fixed
+
+### Security
+
+## 0.4.0-dev-3 - 2024-02-18
+
+### Added
+
+- Reification. Explicit `SafeType` for algebras.
+- Integer division algebras.
+- Float32 geometries.
+- New Attributes-kt module that could be used as stand-alone. It declares. type-safe attributes containers.
+- Explicit `mutableStructureND` builders for mutable structures.
+- `Buffer.asList()` zero-copy transformation.
+- Wasm support.
+- Parallel implementation of `LinearSpace` for Float64
+- Parallel buffer factories
+
+### Changed
+
+- Buffer copy removed from API (added as an extension).
+- Default naming for algebra and buffers now uses IntXX/FloatXX notation instead of Java types.
+- Remove unnecessary inlines in basic algebras.
+- QuaternionField -> QuaternionAlgebra and does not implement `Field` anymore since it is non-commutative
+- kmath-geometry is split into `euclidean2d` and `euclidean3d`
+- Features replaced with Attributes.
+- Transposed refactored.
+- Kmath-memory is moved on top of core.
+
+### Deprecated
+
+- ND4J engine
+
+### Removed
+
+- `asPolynomial` function due to scope pollution
+- Codegend for ejml (450 lines of codegen for 1000 lines of code is too much)
+
+### Fixed
+
+- Median statistics
+- Complex power of negative real numbers
+- Add proper mutability for MutableBufferND rows and columns
+- Generic Float32 and Float64 vectors are used in geometry algebras.
+
+## 0.3.1 - 2023-04-09
+
+### Added
+
+- Wasm support for `memory`, `core`, `complex` and `functions` modules.
+- Generic builders for `BufferND` and `MutableBufferND`
+- `NamedMatrix` - matrix with symbol-based indexing
+- `Expression` with default arguments
+- Type-aliases for numbers like `Float64`
+- Autodiff for generic algebra elements in core!
+- Algebra now has an obligatory `bufferFactory` (#477).
+
+### Changed
+
+- Removed marker `Vector` type for geometry
+- Geometry uses type-safe angles
+- Tensor operations switched to prefix notation
+- Row-wise and column-wise ND shapes in the core
+- Shape is read-only
+- Major refactor of tensors (only minor API changes)
+- Kotlin 1.8.20
+- `LazyStructure` `deffered` -> `async` to comply with coroutines code style
+- Default `dot` operation in tensor algebra no longer support broadcasting. Instead `matmul` operation is added
+  to `DoubleTensorAlgebra`.
+- Multik went MPP
+
+### Removed
+
+- Trajectory moved to https://github.com/SciProgCentre/maps-kt
+- Polynomials moved to https://github.com/SciProgCentre/kmath-polynomial
+
+## 0.3.0
+
+### Added
+
+- `ScaleOperations` interface
+- `Field` extends `ScaleOperations`
+- Basic integration API
+- Basic MPP distributions and samplers
+- `bindSymbolOrNull`
+- Blocking chains and Statistics
+- Multiplatform integration
+- Integration for any Field element
+- Extended operations for ND4J fields
+- Jupyter Notebook integration module (kmath-jupyter)
+- `@PerformancePitfall` annotation to mark possibly slow API
+- Unified architecture for Integration and Optimization using features.
+- `BigInt` operation performance improvement and fixes by @zhelenskiy (#328)
+- Integration between `MST` and Symja `IExpr`
+- Complex power
+- Separate methods for UInt, Int and Number powers. NaN safety.
+- Tensorflow prototype
+- `ValueAndErrorField`
+- MST compilation to WASM: #286
+- Jafama integration: #176
+- `contentEquals` with tolerance: #364
+- Compilation to TeX for MST: #254
+
+### Changed
+
+- Annotations moved to `space.kscience.kmath`
+- Exponential operations merged with hyperbolic functions
+- Space is replaced by Group. Space is reserved for vector spaces.
+- VectorSpace is now a vector space
+- Buffer factories for primitives moved to MutableBuffer.Companion
+- Rename `NDStructure` and `NDAlgebra` to `StructureND` and `AlgebraND` respectively
+- `Real` -> `Double`
+- DataSets are moved from functions to core
+- Redesign advanced Chain API
+- Redesign `MST`. Remove `MstExpression`.
+- Move `MST` to core
+- Separated benchmarks and examples
+- Rewrite `kmath-ejml` without `ejml-simple` artifact, support sparse matrices
+- Promote stability of kmath-ast and kmath-kotlingrad to EXPERIMENTAL.
+- ColumnarData returns nullable column
+- `MST` is made sealed interface
+- Replace `MST.Symbolic` by `Symbol`, `Symbol` now implements MST
+- Remove Any restriction on polynomials
+- Add `out` variance to type parameters of `StructureND` and its implementations where possible
+- Rename `DifferentiableMstExpression` to `KotlingradExpression`
+- `FeatureSet` now accepts only `Feature`. It is possible to override keys and use interfaces.
+- Use `Symbol` factory function instead of `StringSymbol`
+- New discoverability pattern: `<Type>.algebra.<nd/etc>`
+- Adjusted commons-math API for linear solvers to match conventions.
+- Buffer algebra does not require size anymore
+- Operations -> Ops
+- Default Buffer and ND algebras are now Ops and lack neutral elements (0, 1) as well as algebra-level shapes.
+- Tensor algebra takes read-only structures as input and inherits AlgebraND
+- `UnivariateDistribution` renamed to `Distribution1D`
+- Rework of histograms.
+- `UnivariateFunction` -> `Function1D`, `MultivariateFunction` -> `FunctionND`
+
+### Deprecated
+
+- Specialized `DoubleBufferAlgebra`
+
+### Removed
+
+- Nearest in Domain. To be implemented in geometry package.
+- Number multiplication and division in main Algebra chain
+- `contentEquals` from Buffer. It moved to the companion.
+- MSTExpression
+- Expression algebra builders
+- Complex and Quaternion no longer are elements.
+- Second generic from DifferentiableExpression
+- Algebra elements are completely removed. Use algebra contexts instead.
+
+### Fixed
+
+- Ring inherits RingOperations, not GroupOperations
+- Univariate histogram filling
+
+## 0.2.0
+
+### Added
+
+- `fun` annotation for SAM interfaces in library
+- Explicit `public` visibility for all public APIs
+- 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 (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`
+- New `MatrixFeature` interfaces for matrix decompositions
+- Basic Quaternion vector support in `kmath-complex`.
+
+### Changed
+
+- Package changed from `scientifik` to `space.kscience`
+- Gradle version: 6.6 -> 6.8.2
+- Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`)
+- `Polynomial` secondary constructor made function
+- Kotlin version: 1.3.72 -> 1.4.30
+- `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<Float64>` 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`
+- `NumericAlgebra` moved outside of regular algebra chain (`Ring` no longer implements it).
+- Features moved to NDStructure and became transparent.
+- Capitalization of LUP in many names changed to Lup.
+- Refactored `NDStructure` algebra to be more simple, preferring under-the-hood conversion to explicit NDStructure types
+- Refactor histograms. They are marked as prototype
+- `Complex` and related features moved to a separate module `kmath-complex`
+- Refactor AlgebraElement
+- `symbol` method in `Algebra` renamed to `bindSymbol` to avoid ambiguity
+- Add `out` projection to `Buffer` generic
+
+### 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`
+- `Real` class
+- StructureND identity and equals
+
+### Fixed
+
+- `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140)
+
+## 0.1.4
+
+### Added
+
 - Functional Expressions API
 - Mathematical Syntax Tree, its interpreter and API
 - String to MST parser (https://github.com/mipt-npm/kmath/pull/120)
 - MST to JVM bytecode translator (https://github.com/mipt-npm/kmath/pull/94)
 - FloatBuffer (specialized MutableBuffer over FloatArray)
 - FlaggedBuffer to associate primitive numbers buffer with flags (to mark values infinite or missing, etc.)
-- Specialized builder functions for all primitive buffers like `IntBuffer(25) { it + 1 }` (https://github.com/mipt-npm/kmath/pull/125)
+- Specialized builder functions for all primitive buffers
+  like `IntBuffer(25) { it + 1 }` (https://github.com/mipt-npm/kmath/pull/125)
 - Interface `NumericAlgebra` where `number` operation is available to convert numbers to algebraic elements
-- Inverse trigonometric functions support in ExtendedField (`asin`, `acos`, `atan`) (https://github.com/mipt-npm/kmath/pull/114)
+- Inverse trigonometric functions support in
+  ExtendedField (`asin`, `acos`, `atan`) (https://github.com/mipt-npm/kmath/pull/114)
 - New space extensions: `average` and `averageWith`
 - Local coding conventions
 - Geometric Domains API in `kmath-core`
@@ -20,10 +251,12 @@
 - Norm support for `Complex`
 
 ### Changed
+
 - `readAsMemory` now has `throws IOException` in JVM signature.
 - Several functions taking functional types were made `inline`.
 - Several functions taking functional types now have `callsInPlace` contracts.
-- BigInteger and BigDecimal algebra: JBigDecimalField has companion object with default math context; minor optimizations
+- BigInteger and BigDecimal algebra: JBigDecimalField has companion object with default math context; minor
+  optimizations
 - `power(T, Int)` extension function has preconditions and supports `Field<T>`
 - Memory objects have more preconditions (overflow checking)
 - `tg` function is renamed to `tan` (https://github.com/mipt-npm/kmath/pull/114)
@@ -31,6 +264,7 @@
 - Moved probability distributions to commons-rng and to `kmath-prob`
 
 ### Fixed
+
 - Missing copy method in Memory implementation on JS (https://github.com/mipt-npm/kmath/pull/106)
 - D3.dim value in `kmath-dimensions`
 - Multiplication in integer rings in `kmath-core` (https://github.com/mipt-npm/kmath/pull/101)

README.md

@@ -1,109 +1,300 @@
 [![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub)
 [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
-
-![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg)
-
-Bintray:        [ ![Download](https://api.bintray.com/packages/mipt-npm/scientifik/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion)
-
-Bintray-dev:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion)
+![Gradle build](https://github.com/SciProgCentre/kmath/workflows/Gradle%20build/badge.svg)
+[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22)
+[![Space](https://img.shields.io/badge/dynamic/xml?color=orange&label=Space&query=//metadata/versioning/latest&url=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fspace%2Fkscience%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/)
 
 # KMath
-Could be pronounced as `key-math`.
-The Kotlin MATHematics library is intended as a Kotlin-based analog to Python's `numpy` library. In contrast to `numpy` and `scipy` it is modular and has a lightweight core.
 
-## Publications
+Could be pronounced as `key-math`. The **K**otlin **Math**ematics 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.
+
+[Documentation site](https://SciProgCentre.github.io/kmath/)
+
+## 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)
+* [A talk at KotlinConf 2019 about using kotlin for science](https://youtu.be/LI_5TZ7tnOE?si=4LknX41gl_YeUbIe)
+* [A talk on architecture at Joker-2021 (in Russian)](https://youtu.be/1bZ2doHiRRM?si=9w953ro9yu98X_KJ)
+* [The same talk in English](https://youtu.be/yP5DIc2fVwQ?si=louZzQ1dcXV6gP10)
+* [A seminar on tensor API](https://youtu.be/0H99wUs0xTM?si=6c__04jrByFQtVpo)
 
 # 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, Native and
+  Wasm).
 * 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 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 `kmath-for-real`, which will give better
+  experience for those, who want to work with specific types.
 
-## Features
+## Features and stability
 
-Actual feature list is [here](doc/features.md)
+KMath is a modular library. Different modules provide different features with different API stability guarantees. All
+core modules are released with the same version, but with different API change policy. The features are described in
+module definitions below. The module stability could have the following levels:
 
-* **Algebra**
-    * Algebraic structures like rings, spaces and field (**TODO** add example to wiki)
-    * Basic linear algebra operations (sums, products, etc.), backed by the `Space` API.
-    * Complex numbers backed by the `Field` API (meaning that they will be usable in any structure like vectors and N-dimensional arrays).
-    * Advanced linear algebra operations like matrix inversion and LU decomposition.
+* **PROTOTYPE**. On this level there are no compatibility guarantees. All methods and classes form those modules could
+  break any moment. You can still use it, but be sure to fix the specific version.
+* **EXPERIMENTAL**. The general API is decided, but some changes could be made. Volatile API is marked
+  with `@UnstableKMathAPI` or other stability warning annotations.
+* **DEVELOPMENT**. API breaking generally follows semantic versioning ideology. There could be changes in minor
+  versions, but not in patch versions. API is protected
+  with [binary-compatibility-validator](https://github.com/Kotlin/binary-compatibility-validator) tool.
+* **STABLE**. The API stabilized. Breaking changes are allowed only in major releases.
 
-* **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).
+## Modules
 
-* **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.
+### [attributes-kt](attributes-kt)
+> An API and basic implementation for arranging objects in a continuous memory block.
+>
+> **Maturity**: DEVELOPMENT
 
-* **Streaming** Streaming operations on mathematical objects and objects buffers.
+### [benchmarks](benchmarks)
+>
+> **Maturity**: EXPERIMENTAL
 
-* **Type-safe dimensions** Type-safe dimensions for matrix operations.
+### [examples](examples)
+>
+> **Maturity**: EXPERIMENTAL
 
-* **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.
-                           
-* **Koma wrapper** [Koma](https://github.com/kyonifer/koma) is a well established numerics library in Kotlin, specifically linear algebra.
-The plan is to have wrappers for koma implementations for compatibility with kmath API.
+### [kmath-ast](kmath-ast)
+>
+> **Maturity**: EXPERIMENTAL
+>
+> **Features:**
+> - [expression-language](kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/parser.kt) : Expression language and its parser
+> - [mst-jvm-codegen](kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/asm.kt) : Dynamic MST to JVM bytecode compiler
+> - [mst-js-codegen](kmath-ast/src/jsMain/kotlin/space/kscience/kmath/estree/estree.kt) : Dynamic MST to JS compiler
+> - [rendering](kmath-ast/src/commonMain/kotlin/space/kscience/kmath/ast/rendering/MathRenderer.kt) : Extendable MST rendering
 
-## Planned features
 
-* **Messaging** A mathematical notation to support multi-language and multi-node communication for mathematical tasks.
+### [kmath-commons](kmath-commons)
+> Commons math binding for kmath
+>
+> **Maturity**: EXPERIMENTAL
 
-* **Array statistics** 
+### [kmath-complex](kmath-complex)
+> Complex numbers and quaternions.
+>
+> **Maturity**: PROTOTYPE
+>
+> **Features:**
+> - [complex](kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Complex.kt) : Complex numbers operations
+> - [quaternion](kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Quaternion.kt) : Quaternions and their composition
 
-* **Integration** Univariate and multivariate integration framework.
 
-* **Probability and distributions**
+### [kmath-core](kmath-core)
+> Core classes, algebra definitions, basic linear algebra
+>
+> **Maturity**: DEVELOPMENT
+>
+> **Features:**
+> - [algebras](kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt) : Algebraic structures like rings, spaces and fields.
+> - [nd](kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/StructureND.kt) : Many-dimensional structures and operations on them.
+> - [linear](kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt) : Basic linear algebra operations (sums, products, etc.), backed by the `Space` API. Advanced linear algebra operations like matrix inversion and LU decomposition.
+> - [buffers](kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
+> - [expressions](kmath-core/src/commonMain/kotlin/space/kscience/kmath/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.
+> - [domains](kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains) : Domains
+> - [autodiff](kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
+> - [Parallel linear algebra](kmath-core/#) : Parallel implementation for `LinearAlgebra`
+
+
+### [kmath-coroutines](kmath-coroutines)
+>
+> **Maturity**: EXPERIMENTAL
+
+### [kmath-dimensions](kmath-dimensions)
+> A proof of concept module for adding type-safe dimensions to structures
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-ejml](kmath-ejml)
+>
+> **Maturity**: PROTOTYPE
+>
+> **Features:**
+> - [ejml-vector](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt) : Point implementations.
+> - [ejml-matrix](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt) : Matrix implementation.
+> - [ejml-linear-space](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt) : LinearSpace implementations.
+
+
+### [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:**
+> - [DoubleVector](kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt) : Numpy-like operations for Buffers/Points
+> - [DoubleMatrix](kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleMatrix.kt) : Numpy-like operations for 2d real structures
+> - [grids](kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/structures/grids.kt) : Uniform grid generators
+
+
+### [kmath-functions](kmath-functions)
+> Functions, integration and interpolation
+>
+> **Maturity**: EXPERIMENTAL
+>
+> **Features:**
+> - [piecewise](kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt) : Piecewise functions.
+> - [polynomials](kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt) : Polynomial functions.
+> - [linear interpolation](kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt) : Linear XY interpolator.
+> - [spline interpolation](kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt) : Cubic spline XY interpolator.
+> - [integration](kmath-functions/#) : Univariate and multivariate quadratures
+
+
+### [kmath-geometry](kmath-geometry)
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-histograms](kmath-histograms)
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-jafama](kmath-jafama)
+> Jafama integration module
+>
+> **Maturity**: DEPRECATED
+>
+> **Features:**
+> - [jafama-double](kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/) : Double ExtendedField implementations based on Jafama
+
+
+### [kmath-jupyter](kmath-jupyter)
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-kotlingrad](kmath-kotlingrad)
+> Kotlin∇ integration module
+>
+> **Maturity**: EXPERIMENTAL
+>
+> **Features:**
+> - [differentiable-mst-expression](kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/KotlingradExpression.kt) : MST based DifferentiableExpression.
+> - [scalars-adapters](kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/scalarsAdapters.kt) : Conversions between Kotlin∇'s SFun and MST
+
+
+### [kmath-memory](kmath-memory)
+> An API and basic implementation for arranging objects in a continuous memory block.
+>
+> **Maturity**: DEVELOPMENT
+
+### [kmath-multik](kmath-multik)
+> JetBrains Multik connector
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-nd4j](kmath-nd4j)
+> ND4J NDStructure implementation and according NDAlgebra classes
+>
+> **Maturity**: DEPRECATED
+>
+> **Features:**
+> - [nd4jarraystructure](kmath-nd4j/#) : NDStructure wrapper for INDArray
+> - [nd4jarrayrings](kmath-nd4j/#) : Rings over Nd4jArrayStructure of Int and Long
+> - [nd4jarrayfields](kmath-nd4j/#) : Fields over Nd4jArrayStructure of Float and Double
+
+
+### [kmath-optimization](kmath-optimization)
+>
+> **Maturity**: EXPERIMENTAL
+
+### [kmath-stat](kmath-stat)
+>
+> **Maturity**: EXPERIMENTAL
+
+### [kmath-symja](kmath-symja)
+> Symja integration module
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-tensorflow](kmath-tensorflow)
+> Google tensorflow connector
+>
+> **Maturity**: PROTOTYPE
+
+### [kmath-tensors](kmath-tensors)
+>
+> **Maturity**: PROTOTYPE
+>
+> **Features:**
+> - [tensor algebra](kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorAlgebra.kt) : Basic linear algebra operations on tensors (plus, dot, etc.)
+> - [tensor algebra with broadcasting](kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt) : Basic linear algebra operations implemented with broadcasting.
+> - [linear algebra operations](kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt) : Advanced linear algebra operations like LU decomposition, SVD, etc.
+
+
+### [kmath-viktor](kmath-viktor)
+> Binding for https://github.com/JetBrains-Research/viktor
+>
+> **Maturity**: DEPRECATED
+
+### [test-utils](test-utils)
+>
+> **Maturity**: EXPERIMENTAL
 
-* **Fitting** Non-linear curve fitting facilities
 
 ## Multi-platform support
 
-KMath is developed as a multi-platform library, which means that most of interfaces are declared in the [common module](kmath-core/src/commonMain). Implementation is also done in the common module wherever possible. In some cases, features are delegated to platform-specific implementations even if they could be done in the common module for performance reasons. Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome.
+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, 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 of performance is one of the 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 a 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):
+## Requirements
+
+KMath currently relies on JDK 11 for compilation and execution of Kotlin-JVM part. We recommend using GraalVM-CE or
+Oracle GraalVM for execution to get better performance.
+
+### Repositories
+
+Release and development artifacts are accessible from mipt-npm [Space](https://www.jetbrains.com/space/)
+repository `https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven` (see documentation of
+[Kotlin Multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) for more details). The repository could
+be reached through [repo.kotlin.link](https://repo.kotlin.link) proxy:
 
 ```kotlin
-repositories{
-    maven("https://dl.bintray.com/mipt-npm/scientifik")
+repositories {
+    maven("https://repo.kotlin.link")
 }
 
-dependencies{
-    api("scientifik:kmath-core:${kmathVersion}")
-    //api("scientifik:kmath-core-jvm:${kmathVersion}") for jvm-specific version
+dependencies {
+    api("space.kscience:kmath-core:$version")
+    // api("space.kscience:kmath-core-jvm:$version") for jvm-specific version
 }
 ```
 
-Gradle `6.0+` is required for multiplatform artifacts.
-
-### Development
-
-Development builds are accessible from the reposirtory 
-```kotlin
-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 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 [good first issue](hhttps://github.com/SciProgCentre/kmath/issues?q=is%3Aissue+is%3Aopen+label%3A%22good+first+issue%22)
+label.

attributes-kt/README.md

@@ -0,0 +1,21 @@
+# Module attributes-kt
+
+
+
+## Usage
+
+## Artifact:
+
+The Maven coordinates of this project are `space.kscience:attributes-kt:0.1.0`.
+
+**Gradle Kotlin DSL:**
+```kotlin
+repositories {
+    maven("https://repo.kotlin.link")
+    mavenCentral()
+}
+
+dependencies {
+    implementation("space.kscience:attributes-kt:0.1.0")
+}
+```

attributes-kt/api/attributes-kt.api

@@ -0,0 +1,104 @@
+public abstract interface class space/kscience/attributes/Attribute {
+}
+
+public abstract interface class space/kscience/attributes/AttributeContainer {
+	public abstract fun getAttributes ()Lspace/kscience/attributes/Attributes;
+}
+
+public abstract interface class space/kscience/attributes/AttributeScope {
+}
+
+public abstract interface class space/kscience/attributes/AttributeWithDefault : space/kscience/attributes/Attribute {
+	public abstract fun getDefault ()Ljava/lang/Object;
+}
+
+public abstract interface class space/kscience/attributes/Attributes {
+	public static final field Companion Lspace/kscience/attributes/Attributes$Companion;
+	public abstract fun equals (Ljava/lang/Object;)Z
+	public fun get (Lspace/kscience/attributes/Attribute;)Ljava/lang/Object;
+	public abstract fun getContent ()Ljava/util/Map;
+	public fun getKeys ()Ljava/util/Set;
+	public abstract fun hashCode ()I
+	public abstract fun toString ()Ljava/lang/String;
+}
+
+public final class space/kscience/attributes/Attributes$Companion {
+	public final fun equals (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attributes;)Z
+	public final fun getEMPTY ()Lspace/kscience/attributes/Attributes;
+}
+
+public final class space/kscience/attributes/AttributesBuilder : space/kscience/attributes/Attributes {
+	public final fun add (Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)V
+	public final fun build ()Lspace/kscience/attributes/Attributes;
+	public fun equals (Ljava/lang/Object;)Z
+	public fun getContent ()Ljava/util/Map;
+	public fun hashCode ()I
+	public final fun invoke (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)V
+	public final fun put (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)V
+	public final fun putAll (Lspace/kscience/attributes/Attributes;)V
+	public final fun remove (Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)V
+	public final fun set (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)V
+	public fun toString ()Ljava/lang/String;
+}
+
+public final class space/kscience/attributes/AttributesBuilderKt {
+	public static final fun Attributes (Lkotlin/jvm/functions/Function1;)Lspace/kscience/attributes/Attributes;
+}
+
+public final class space/kscience/attributes/AttributesKt {
+	public static final fun Attributes (Lspace/kscience/attributes/Attribute;)Lspace/kscience/attributes/Attributes;
+	public static final fun Attributes (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+	public static final fun getOrDefault (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/AttributeWithDefault;)Ljava/lang/Object;
+	public static final fun isEmpty (Lspace/kscience/attributes/Attributes;)Z
+	public static final fun modified (Lspace/kscience/attributes/Attributes;Lkotlin/jvm/functions/Function1;)Lspace/kscience/attributes/Attributes;
+	public static final fun plus (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attributes;)Lspace/kscience/attributes/Attributes;
+	public static final fun withAttribute (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attribute;)Lspace/kscience/attributes/Attributes;
+	public static final fun withAttribute (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+	public static final fun withAttributeElement (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+	public static final fun withoutAttribute (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attribute;)Lspace/kscience/attributes/Attributes;
+	public static final fun withoutAttributeElement (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+}
+
+public abstract interface class space/kscience/attributes/FlagAttribute : space/kscience/attributes/Attribute {
+}
+
+public abstract class space/kscience/attributes/PolymorphicAttribute : space/kscience/attributes/Attribute {
+	public synthetic fun <init> (Lkotlin/reflect/KType;Lkotlin/jvm/internal/DefaultConstructorMarker;)V
+	public fun equals (Ljava/lang/Object;)Z
+	public final fun getType-V0oMfBY ()Lkotlin/reflect/KType;
+	public fun hashCode ()I
+}
+
+public final class space/kscience/attributes/PolymorphicAttributeKt {
+	public static final fun get (Lspace/kscience/attributes/Attributes;Lkotlin/jvm/functions/Function0;)Ljava/lang/Object;
+	public static final fun set (Lspace/kscience/attributes/AttributesBuilder;Lkotlin/jvm/functions/Function0;Ljava/lang/Object;)V
+}
+
+public final class space/kscience/attributes/SafeType {
+	public static final synthetic fun box-impl (Lkotlin/reflect/KType;)Lspace/kscience/attributes/SafeType;
+	public static fun constructor-impl (Lkotlin/reflect/KType;)Lkotlin/reflect/KType;
+	public fun equals (Ljava/lang/Object;)Z
+	public static fun equals-impl (Lkotlin/reflect/KType;Ljava/lang/Object;)Z
+	public static final fun equals-impl0 (Lkotlin/reflect/KType;Lkotlin/reflect/KType;)Z
+	public final fun getKType ()Lkotlin/reflect/KType;
+	public fun hashCode ()I
+	public static fun hashCode-impl (Lkotlin/reflect/KType;)I
+	public fun toString ()Ljava/lang/String;
+	public static fun toString-impl (Lkotlin/reflect/KType;)Ljava/lang/String;
+	public final synthetic fun unbox-impl ()Lkotlin/reflect/KType;
+}
+
+public final class space/kscience/attributes/SafeTypeKt {
+	public static final fun getKClass-X0YbwmU (Lkotlin/reflect/KType;)Lkotlin/reflect/KClass;
+}
+
+public abstract interface class space/kscience/attributes/SetAttribute : space/kscience/attributes/Attribute {
+}
+
+public abstract interface annotation class space/kscience/attributes/UnstableAttributesAPI : java/lang/annotation/Annotation {
+}
+
+public abstract interface class space/kscience/attributes/WithType {
+	public abstract fun getType-V0oMfBY ()Lkotlin/reflect/KType;
+}
+

attributes-kt/build.gradle.kts

@@ -0,0 +1,20 @@
+plugins {
+    id("space.kscience.gradle.mpp")
+    `maven-publish`
+}
+
+version = rootProject.extra.get("attributesVersion").toString()
+
+kscience {
+    jvm()
+    js()
+    native()
+    wasm()
+}
+
+readme {
+    maturity = space.kscience.gradle.Maturity.DEVELOPMENT
+    description = """
+        An API and basic implementation for arranging objects in a continuous memory block.
+    """.trimIndent()
+}

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt

@@ -0,0 +1,29 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+/**
+ * A marker interface for an attribute. Attributes are used as keys to access contents of type [T] in the container.
+ */
+public interface Attribute<T>
+
+/**
+ * An attribute that could be either present or absent
+ */
+public interface FlagAttribute : Attribute<Unit>
+
+/**
+ * An attribute with a default value
+ */
+public interface AttributeWithDefault<T> : Attribute<T> {
+    public val default: T
+}
+
+/**
+ * Attribute containing a set of values
+ */
+public interface SetAttribute<V> : Attribute<Set<V>>
+

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/AttributeContainer.kt

@@ -0,0 +1,20 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+/**
+ * A container for [Attributes]
+ */
+public interface AttributeContainer {
+    public val attributes: Attributes
+}
+
+/**
+ * A scope, where attribute keys could be resolved.
+ * [O] is used only to resolve types in compile-time.
+ */
+public interface AttributeScope<O>
+

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt

@@ -0,0 +1,157 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+/**
+ * A set of attributes. The implementation must guarantee that [content] keys correspond to their value types.
+ */
+public interface Attributes {
+    /**
+     * Raw content for this [Attributes]
+     */
+    public val content: Map<out Attribute<*>, Any?>
+
+    /**
+     * Attribute keys contained in this [Attributes]
+     */
+    public val keys: Set<Attribute<*>> get() = content.keys
+
+    /**
+     * Provide an attribute value. Return null if attribute is not present or if its value is null.
+     */
+    @Suppress("UNCHECKED_CAST")
+    public operator fun <T> get(attribute: Attribute<T>): T? = content[attribute] as? T
+
+    override fun toString(): String
+    override fun equals(other: Any?): Boolean
+    override fun hashCode(): Int
+
+    public companion object {
+        public val EMPTY: Attributes = object : Attributes {
+            override val content: Map<out Attribute<*>, Any?> get() = emptyMap()
+
+            override fun toString(): String = "Attributes.EMPTY"
+
+            override fun equals(other: Any?): Boolean = (other as? Attributes)?.isEmpty() ?: false
+
+            override fun hashCode(): Int = Unit.hashCode()
+        }
+
+        public fun equals(a1: Attributes, a2: Attributes): Boolean =
+            a1.keys == a2.keys && a1.keys.all { a1[it] == a2[it] }
+    }
+}
+
+internal class MapAttributes(override val content: Map<out Attribute<*>, Any?>) : Attributes {
+    override fun toString(): String = "Attributes(value=${content.entries})"
+    override fun equals(other: Any?): Boolean = other is Attributes && Attributes.equals(this, other)
+    override fun hashCode(): Int = content.hashCode()
+}
+
+public fun Attributes.isEmpty(): Boolean = keys.isEmpty()
+
+/**
+ * Get attribute value or default
+ */
+public fun <T> Attributes.getOrDefault(attribute: AttributeWithDefault<T>): T = get(attribute) ?: attribute.default
+
+/**
+ * Check if there is an attribute that matches given key by type and adheres to [predicate].
+ */
+@Suppress("UNCHECKED_CAST")
+public inline fun <T, reified A : Attribute<T>> Attributes.hasAny(predicate: (value: T) -> Boolean): Boolean =
+    content.any { (mapKey, mapValue) -> mapKey is A && predicate(mapValue as T) }
+
+/**
+ * Check if there is an attribute of given type (subtypes included)
+ */
+public inline fun <reified A : Attribute<*>> Attributes.hasAny(): Boolean =
+    content.any { (mapKey, _) -> mapKey is A }
+
+/**
+ * Check if [Attributes] contains a flag. Multiple keys that are instances of a flag could be present
+ */
+public inline fun <reified A : FlagAttribute> Attributes.hasFlag(): Boolean =
+    content.keys.any { it is A }
+
+/**
+ * Create [Attributes] with an added or replaced attribute key.
+ */
+public fun <T, A : Attribute<T>> Attributes.withAttribute(
+    attribute: A,
+    attrValue: T,
+): Attributes = MapAttributes(content + (attribute to attrValue))
+
+public fun <A : Attribute<Unit>> Attributes.withAttribute(attribute: A): Attributes =
+    withAttribute(attribute, Unit)
+
+/**
+ * Create a new [Attributes] by modifying the current one
+ */
+public fun <O> Attributes.modified(block: AttributesBuilder<O>.() -> Unit): Attributes = Attributes<O> {
+    putAll(this@modified)
+    block()
+}
+
+/**
+ * Create new [Attributes] by removing [attribute] key
+ */
+public fun Attributes.withoutAttribute(attribute: Attribute<*>): Attributes = MapAttributes(content.minus(attribute))
+
+/**
+ * Add an element to a [SetAttribute]
+ */
+public fun <T, A : SetAttribute<T>> Attributes.withAttributeElement(
+    attribute: A,
+    attrValue: T,
+): Attributes {
+    val currentSet: Set<T> = get(attribute) ?: emptySet()
+    return MapAttributes(
+        content + (attribute to (currentSet + attrValue))
+    )
+}
+
+/**
+ * Remove an element from [SetAttribute]
+ */
+public fun <T, A : SetAttribute<T>> Attributes.withoutAttributeElement(
+    attribute: A,
+    attrValue: T,
+): Attributes {
+    val currentSet: Set<T> = get(attribute) ?: emptySet()
+    return MapAttributes(content + (attribute to (currentSet - attrValue)))
+}
+
+/**
+ * Create [Attributes] with a single key
+ */
+public fun <T, A : Attribute<T>> Attributes(
+    attribute: A,
+    attrValue: T,
+): Attributes = MapAttributes(mapOf(attribute to attrValue))
+
+/**
+ * Create Attributes with a single [Unit] valued attribute
+ */
+public fun <A : Attribute<Unit>> Attributes(
+    attribute: A,
+): Attributes = MapAttributes(mapOf(attribute to Unit))
+
+/**
+ * Create a new [Attributes] that overlays [other] on top of this set of attributes. New attributes are added.
+ * Existing attribute keys are replaced.
+ */
+public operator fun Attributes.plus(other: Attributes): Attributes = when {
+    isEmpty() -> other
+    other.isEmpty() -> this
+    else -> MapAttributes(content + other.content)
+}
+
+/**
+ * Create a new [Attributes] with removed [key] (if it is present).
+ */
+public operator fun Attributes.minus(key: Attribute<*>): Attributes =
+    if (content.contains(key)) MapAttributes(content.minus(key)) else this

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/AttributesBuilder.kt

@@ -0,0 +1,68 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+/**
+ * A builder for [Attributes].
+ * The builder is not thread safe
+ *
+ * @param O type marker of an owner object, for which these attributes are made
+ */
+public class AttributesBuilder<out O> internal constructor() : Attributes {
+
+    private val map = mutableMapOf<Attribute<*>, Any?>()
+
+    override fun toString(): String = "Attributes(value=${map.entries})"
+    override fun equals(other: Any?): Boolean = other is Attributes && Attributes.equals(this, other)
+    override fun hashCode(): Int = map.hashCode()
+
+    override val content: Map<out Attribute<*>, Any?> get() = map
+
+    public operator fun <T> set(attribute: Attribute<T>, value: T?) {
+        if (value == null) {
+            map.remove(attribute)
+        } else {
+            map[attribute] = value
+        }
+    }
+
+    public operator fun <V> Attribute<V>.invoke(value: V?) {
+        set(this, value)
+    }
+
+    public infix fun <V> Attribute<V>.put(value: V?) {
+        set(this, value)
+    }
+
+    /**
+     * Put all attributes for given [attributes]
+     */
+    public fun putAll(attributes: Attributes) {
+        map.putAll(attributes.content)
+    }
+
+    public infix fun <V> SetAttribute<V>.add(attrValue: V) {
+        val currentSet: Set<V> = get(this) ?: emptySet()
+        map[this] = currentSet + attrValue
+    }
+
+    /**
+     * Remove an element from [SetAttribute]
+     */
+    public infix fun <V> SetAttribute<V>.remove(attrValue: V) {
+        val currentSet: Set<V> = get(this) ?: emptySet()
+        map[this] = currentSet - attrValue
+    }
+
+    public fun build(): Attributes = MapAttributes(map)
+}
+
+/**
+ * Create [Attributes] with a given [builder]
+ * @param O the type for which attributes are built. The type is used only during compilation phase for static extension dispatch
+ */
+public fun <O> Attributes(builder: AttributesBuilder<O>.() -> Unit): Attributes =
+    AttributesBuilder<O>().apply(builder).build()

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/PolymorphicAttribute.kt

@@ -0,0 +1,34 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+/**
+ * An attribute that has a type parameter for value
+ * @param type parameter-type
+ */
+public abstract class PolymorphicAttribute<T>(public val type: SafeType<T>) : Attribute<T> {
+    override fun equals(other: Any?): Boolean = other != null &&
+            (this::class == other::class) &&
+            (other as? PolymorphicAttribute<*>)?.type == this.type
+
+    override fun hashCode(): Int = this::class.hashCode() + type.hashCode()
+}
+
+
+/**
+ * Get a polymorphic attribute using attribute factory
+ */
+@UnstableAttributesAPI
+public operator fun <T> Attributes.get(attributeKeyBuilder: () -> PolymorphicAttribute<T>): T? =
+    get(attributeKeyBuilder())
+
+/**
+ * Set a polymorphic attribute using its factory
+ */
+@UnstableAttributesAPI
+public operator fun <O, T> AttributesBuilder<O>.set(attributeKeyBuilder: () -> PolymorphicAttribute<T>, value: T) {
+    set(attributeKeyBuilder(), value)
+}

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/SafeType.kt

@@ -0,0 +1,35 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+import kotlin.jvm.JvmInline
+import kotlin.reflect.KClass
+import kotlin.reflect.KType
+import kotlin.reflect.typeOf
+
+/**
+ * Safe variant ok Kotlin [KType] that ensures that the type parameter is of the same type as [kType]
+ *
+ * @param kType raw [KType]
+ */
+@JvmInline
+public value class SafeType<out T> @PublishedApi internal constructor(public val kType: KType)
+
+public inline fun <reified T> safeTypeOf(): SafeType<T> = SafeType(typeOf<T>())
+
+/**
+ * Derive Kotlin [KClass] from this type and fail if the type is not a class (should not happen)
+ */
+@Suppress("UNCHECKED_CAST")
+@UnstableAttributesAPI
+public val <T> SafeType<T>.kClass: KClass<T & Any> get() = kType.classifier as KClass<T & Any>
+
+/**
+ * An interface containing [type] for dynamic type checking.
+ */
+public interface WithType<out T> {
+    public val type: SafeType<T>
+}

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/annotations.kt

@@ -0,0 +1,17 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.attributes
+
+/**
+ * Marks declarations that are still experimental in the Attributes-kt APIs, which means that the design of the corresponding
+ * declarations has open issues that may (or may not) lead to their changes in the future. Roughly speaking, there is
+ * a chance of those declarations will be deprecated in the future or the semantics of their behavior may change
+ * in some way that may break some code.
+ */
+@MustBeDocumented
+@Retention(value = AnnotationRetention.BINARY)
+@RequiresOptIn("This API is unstable and could change in future", RequiresOptIn.Level.WARNING)
+public annotation class UnstableAttributesAPI

benchmarks/README.md

@@ -0,0 +1,119 @@
+# BenchmarksResult
+
+## Report for benchmark configuration <code>main</code>
+
+* Run on Java HotSpot(TM) 64-Bit Server VM (build 21.0.4+8-LTS-jvmci-23.1-b41) with Java process:
+
+```
+C:\Users\altavir\scoop\apps\graalvm-oracle-21jdk\current\bin\java.exe -XX:ThreadPriorityPolicy=1 -XX:+UnlockExperimentalVMOptions -XX:+EnableJVMCIProduct -XX:-UnlockExperimentalVMOptions -Dfile.encoding=UTF-8 -Duser.country=US -Duser.language=en -Duser.variant
+```
+* JMH 1.21 was used in `thrpt` mode with 5 warmup iterations by 10 s and 5 measurement iterations by 10 s.
+### [ArrayBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ArrayBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`benchmarkArrayRead`|1.9E+07 &plusmn; 2.3E+05 ops/s|
+|`benchmarkBufferRead`|1.4E+07 &plusmn; 8.7E+05 ops/s|
+|`nativeBufferRead`|1.4E+07 &plusmn; 1.3E+06 ops/s|
+### [BigIntBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`jvmAdd`|5.1E+07 &plusmn; 1.3E+06 ops/s|
+|`jvmAddLarge`|5.1E+04 &plusmn; 8.2E+02 ops/s|
+|`jvmMultiply`|8.5E+07 &plusmn; 9.7E+06 ops/s|
+|`jvmMultiplyLarge`|2.5E+02 &plusmn; 15 ops/s|
+|`jvmParsing10`|8.7E+06 &plusmn; 5.1E+05 ops/s|
+|`jvmParsing16`|6.4E+06 &plusmn; 1.8E+05 ops/s|
+|`jvmPower`|28 &plusmn; 0.79 ops/s|
+|`jvmSmallAdd`|7.0E+07 &plusmn; 4.3E+06 ops/s|
+|`kmAdd`|4.8E+07 &plusmn; 2.2E+06 ops/s|
+|`kmAddLarge`|3.5E+04 &plusmn; 3.7E+03 ops/s|
+|`kmMultiply`|6.7E+07 &plusmn; 1.5E+07 ops/s|
+|`kmMultiplyLarge`|54 &plusmn; 4.2 ops/s|
+|`kmParsing10`|4.5E+06 &plusmn; 8.3E+04 ops/s|
+|`kmParsing16`|4.9E+06 &plusmn; 1.1E+05 ops/s|
+|`kmPower`|10 &plusmn; 0.96 ops/s|
+|`kmSmallAdd`|4.1E+07 &plusmn; 5.9E+05 ops/s|
+### [BufferBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BufferBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`bufferViewReadWrite`|5.8E+06 &plusmn; 1.6E+05 ops/s|
+|`bufferViewReadWriteSpecialized`|5.6E+06 &plusmn; 2.6E+05 ops/s|
+|`complexBufferReadWrite`|6.6E+06 &plusmn; 2.7E+05 ops/s|
+|`doubleArrayReadWrite`|7.5E+06 &plusmn; 1.0E+06 ops/s|
+|`doubleBufferReadWrite`|8.0E+06 &plusmn; 6.7E+05 ops/s|
+### [DotBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`bufferedDot`|1.3 &plusmn; 0.020 ops/s|
+|`cmDot`|0.47 &plusmn; 0.42 ops/s|
+|`cmDotWithConversion`|0.76 &plusmn; 0.13 ops/s|
+|`ejmlDot`|6.7 &plusmn; 0.091 ops/s|
+|`ejmlDotWithConversion`|6.4 &plusmn; 0.82 ops/s|
+|`multikDot`|40 &plusmn; 6.7 ops/s|
+|`parallelDot`|12 &plusmn; 1.8 ops/s|
+|`tensorDot`|1.2 &plusmn; 0.041 ops/s|
+|`tfDot`|5.9 &plusmn; 0.49 ops/s|
+### [ExpressionsInterpretersBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`asmGenericExpression`|29 &plusmn; 1.2 ops/s|
+|`asmPrimitiveExpression`|43 &plusmn; 1.3 ops/s|
+|`asmPrimitiveExpressionArray`|71 &plusmn; 0.38 ops/s|
+|`functionalExpression`|5.6 &plusmn; 0.11 ops/s|
+|`justCalculate`|69 &plusmn; 9.0 ops/s|
+|`mstExpression`|7.1 &plusmn; 0.020 ops/s|
+|`rawExpression`|41 &plusmn; 1.5 ops/s|
+### [IntegrationBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/IntegrationBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`complexIntegration`|3.6E+03 &plusmn; 1.9E+02 ops/s|
+|`doubleIntegration`|3.7E+03 &plusmn; 12 ops/s|
+### [JafamaBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/JafamaBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`core`|38 &plusmn; 0.64 ops/s|
+|`jafama`|52 &plusmn; 0.36 ops/s|
+|`strictJafama`|52 &plusmn; 4.0 ops/s|
+### [MatrixInverseBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`cmLUPInversion`|2.2E+03 &plusmn; 76 ops/s|
+|`ejmlInverse`|1.3E+03 &plusmn; 5.7 ops/s|
+|`kmathLupInversion`|9.5E+02 &plusmn; 1.8E+02 ops/s|
+|`kmathParallelLupInversion`|9.1E+02 &plusmn; 1.4E+02 ops/s|
+### [NDFieldBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`boxingFieldAdd`|7.7 &plusmn; 0.79 ops/s|
+|`multikAdd`|6.5 &plusmn; 0.33 ops/s|
+|`multikInPlaceAdd`|64 &plusmn; 0.79 ops/s|
+|`specializedFieldAdd`|8.0 &plusmn; 0.090 ops/s|
+|`tensorAdd`|9.2 &plusmn; 0.053 ops/s|
+|`tensorInPlaceAdd`|17 &plusmn; 10 ops/s|
+|`viktorAdd`|7.6 &plusmn; 1.2 ops/s|
+### [ViktorBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`doubleFieldAddition`|7.7 &plusmn; 0.34 ops/s|
+|`rawViktor`|5.9 &plusmn; 1.1 ops/s|
+|`viktorFieldAddition`|7.3 &plusmn; 1.1 ops/s|
+### [ViktorLogBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`rawViktorLog`|1.4 &plusmn; 0.076 ops/s|
+|`realFieldLog`|1.3 &plusmn; 0.069 ops/s|
+|`viktorFieldLog`|1.3 &plusmn; 0.032 ops/s|
+
+
+

benchmarks/build.gradle.kts

@@ -0,0 +1,281 @@
+import com.fasterxml.jackson.module.kotlin.jacksonObjectMapper
+import com.fasterxml.jackson.module.kotlin.readValue
+import kotlinx.benchmark.gradle.BenchmarksExtension
+import java.time.LocalDateTime
+import java.time.ZoneId
+import java.util.*
+
+plugins {
+    kotlin("multiplatform")
+    alias(spclibs.plugins.kotlin.plugin.allopen)
+    alias(spclibs.plugins.kotlinx.benchmark)
+}
+
+allOpen.annotation("org.openjdk.jmh.annotations.State")
+sourceSets.register("benchmarks")
+
+repositories {
+    mavenCentral()
+}
+
+kotlin {
+    jvm()
+
+    js(IR) {
+        nodejs()
+    }
+
+    sourceSets {
+        all {
+            languageSettings {
+                progressiveMode = true
+                optIn("kotlin.contracts.ExperimentalContracts")
+                optIn("kotlin.ExperimentalUnsignedTypes")
+                optIn("space.kscience.kmath.UnstableKMathAPI")
+            }
+        }
+
+        val commonMain by getting {
+            dependencies {
+                implementation(project(":kmath-ast"))
+                implementation(project(":kmath-core"))
+                implementation(project(":kmath-coroutines"))
+                implementation(project(":kmath-complex"))
+                implementation(project(":kmath-stat"))
+                implementation(project(":kmath-dimensions"))
+                implementation(project(":kmath-for-real"))
+                implementation(project(":kmath-tensors"))
+                implementation(project(":kmath-multik"))
+                implementation(libs.multik.default)
+                implementation(spclibs.kotlinx.benchmark.runtime)
+            }
+        }
+
+        val jvmMain by getting {
+            dependencies {
+                implementation(project(":kmath-commons"))
+                implementation(project(":kmath-ejml"))
+                implementation(project(":kmath-nd4j"))
+                implementation(project(":kmath-kotlingrad"))
+                implementation(project(":kmath-viktor"))
+                implementation(project(":kmath-jafama"))
+                implementation(projects.kmath.kmathTensorflow)
+                implementation("org.tensorflow:tensorflow-core-platform:0.4.0")
+                implementation("org.nd4j:nd4j-native:1.0.0-M1")
+                //    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")
+            }
+        }
+    }
+}
+
+// Configure benchmark
+benchmark {
+    // Setup configurations
+    targets {
+        register("jvm")
+        register("js")
+    }
+
+    fun kotlinx.benchmark.gradle.BenchmarkConfiguration.commonConfiguration() {
+        warmups = 2
+        iterations = 5
+        iterationTime = 2000
+        iterationTimeUnit = "ms"
+    }
+
+    configurations.register("buffer") {
+        commonConfiguration()
+        include("BufferBenchmark")
+    }
+
+    configurations.register("nd") {
+        commonConfiguration()
+        include("NDFieldBenchmark")
+    }
+
+    configurations.register("dot") {
+        commonConfiguration()
+        include("DotBenchmark")
+    }
+
+    configurations.register("expressions") {
+        // Some extra precision
+        warmups = 2
+        iterations = 10
+        iterationTime = 10
+        iterationTimeUnit = "s"
+        outputTimeUnit = "s"
+        include("ExpressionsInterpretersBenchmark")
+    }
+
+    configurations.register("matrixInverse") {
+        commonConfiguration()
+        include("MatrixInverseBenchmark")
+    }
+
+    configurations.register("bigInt") {
+        commonConfiguration()
+        include("BigIntBenchmark")
+    }
+
+    configurations.register("jafamaDouble") {
+        commonConfiguration()
+        include("JafamaBenchmark")
+    }
+
+    configurations.register("tensorAlgebra") {
+        commonConfiguration()
+        include("TensorAlgebraBenchmark")
+    }
+
+    configurations.register("viktor") {
+        commonConfiguration()
+        include("ViktorBenchmark")
+    }
+
+    configurations.register("viktorLog") {
+        commonConfiguration()
+        include("ViktorLogBenchmark")
+    }
+
+    configurations.register("integration") {
+        commonConfiguration()
+        include("IntegrationBenchmark")
+    }
+}
+
+kotlin {
+    jvmToolchain(11)
+    compilerOptions {
+        optIn.addAll(
+            "space.kscience.kmath.UnstableKMathAPI"
+        )
+    }
+}
+
+
+private data class JmhReport(
+    val jmhVersion: String,
+    val benchmark: String,
+    val mode: String,
+    val threads: Int,
+    val forks: Int,
+    val jvm: String,
+    val jvmArgs: List<String>,
+    val jdkVersion: String,
+    val vmName: String,
+    val vmVersion: String,
+    val warmupIterations: Int,
+    val warmupTime: String,
+    val warmupBatchSize: Int,
+    val measurementIterations: Int,
+    val measurementTime: String,
+    val measurementBatchSize: Int,
+    val params: Map<String, String> = emptyMap(),
+    val primaryMetric: PrimaryMetric,
+    val secondaryMetrics: Map<String, SecondaryMetric>,
+) {
+    interface Metric {
+        val score: Double
+        val scoreError: Double
+        val scoreConfidence: List<Double>
+        val scorePercentiles: Map<Double, Double>
+        val scoreUnit: String
+    }
+
+    data class PrimaryMetric(
+        override val score: Double,
+        override val scoreError: Double,
+        override val scoreConfidence: List<Double>,
+        override val scorePercentiles: Map<Double, Double>,
+        override val scoreUnit: String,
+        val rawDataHistogram: List<List<List<List<Double>>>>? = null,
+        val rawData: List<List<Double>>? = null,
+    ) : Metric
+
+    data class SecondaryMetric(
+        override val score: Double,
+        override val scoreError: Double,
+        override val scoreConfidence: List<Double>,
+        override val scorePercentiles: Map<Double, Double>,
+        override val scoreUnit: String,
+        val rawData: List<List<Double>>,
+    ) : Metric
+}
+
+readme {
+    maturity = space.kscience.gradle.Maturity.EXPERIMENTAL
+
+    val jsonMapper = jacksonObjectMapper()
+
+    fun noun(number: Number, singular: String, plural: String) = if (number.toLong() == 1L) singular else plural
+
+    extensions.findByType(BenchmarksExtension::class.java)?.configurations?.forEach { cfg ->
+        val propertyName =
+            "benchmark${cfg.name.replaceFirstChar { if (it.isLowerCase()) it.titlecase(Locale.getDefault()) else it.toString() }}"
+
+        logger.info("Processing benchmark data from benchmark ${cfg.name} into readme property $propertyName")
+
+        val launches = layout.buildDirectory.dir("reports/benchmarks/${cfg.name}").get().asFile
+        if (!launches.exists()) return@forEach
+
+        property(propertyName) {
+            val resDirectory = launches.listFiles()?.maxByOrNull {
+                LocalDateTime.parse(it.name).atZone(ZoneId.systemDefault()).toInstant()
+            }
+
+            if (resDirectory == null || !(resDirectory.resolve("jvm.json")).exists()) {
+                "> **Can't find appropriate benchmark data. Try generating readme files after running benchmarks**."
+            } else {
+                val reports: List<JmhReport> =
+                    jsonMapper.readValue<List<JmhReport>>(resDirectory.resolve("jvm.json"))
+
+                buildString {
+                    appendLine("## Report for benchmark configuration <code>${cfg.name}</code>")
+                    appendLine()
+                    val first = reports.first()
+
+                    appendLine("* Run on ${first.vmName} (build ${first.vmVersion}) with Java process:")
+                    appendLine()
+                    appendLine("```")
+                    appendLine(
+                        "${first.jvm} ${
+                            first.jvmArgs.joinToString(" ")
+                        }"
+                    )
+                    appendLine("```")
+
+                    appendLine(
+                        "* JMH ${first.jmhVersion} was used in `${first.mode}` mode with ${first.warmupIterations} warmup ${
+                            noun(first.warmupIterations, "iteration", "iterations")
+                        } by ${first.warmupTime} and ${first.measurementIterations} measurement ${
+                            noun(first.measurementIterations, "iteration", "iterations")
+                        } by ${first.measurementTime}."
+                    )
+
+                    reports.groupBy { it.benchmark.substringBeforeLast(".") }.forEach { (cl, compare) ->
+                        appendLine("### [${cl.substringAfterLast(".")}](src/jvmMain/kotlin/${cl.replace(".","/")}.kt)")
+                        appendLine()
+                        appendLine("| Benchmark | Score |")
+                        appendLine("|:---------:|:-----:|")
+                        compare.forEach { report ->
+                            val benchmarkName = report.benchmark.substringAfterLast(".")
+                            val score = String.format("%.2G", report.primaryMetric.score)
+                            val error = String.format("%.2G", report.primaryMetric.scoreError)
+
+                            appendLine("|`$benchmarkName`|$score &plusmn; $error ${report.primaryMetric.scoreUnit}|")
+                        }
+                    }
+                }
+            }
+        }
+    }
+}

benchmarks/docs/README-TEMPLATE.md

@@ -0,0 +1,5 @@
+# BenchmarksResult
+
+${benchmarkMain}
+
+

benchmarks/src/jsMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt

@@ -0,0 +1,108 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.expressions.*
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.bindSymbol
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Float64
+import kotlin.math.sin
+import kotlin.random.Random
+import space.kscience.kmath.estree.compileToExpression as estreeCompileToExpression
+import space.kscience.kmath.wasm.compileToExpression as wasmCompileToExpression
+
+@State(Scope.Benchmark)
+class ExpressionsInterpretersBenchmark {
+    /**
+     * Benchmark case for [Expression] created with [expressionInExtendedField].
+     */
+    @Benchmark
+    fun functionalExpression(blackhole: Blackhole) = invokeAndSum(functional, blackhole)
+
+    /**
+     * Benchmark case for [Expression] created with [toExpression].
+     */
+    @Benchmark
+    fun mstExpression(blackhole: Blackhole) = invokeAndSum(mst, blackhole)
+
+    /**
+     * Benchmark case for [Expression] created with [compileToExpression].
+     */
+    @Benchmark
+    fun wasmExpression(blackhole: Blackhole) = invokeAndSum(wasm, blackhole)
+
+    /**
+     * Benchmark case for [Expression] created with [compileToExpression].
+     */
+    @Benchmark
+    fun estreeExpression(blackhole: Blackhole) = invokeAndSum(estree, blackhole)
+
+    /**
+     * Benchmark case for [Expression] implemented manually with `kotlin.math` functions.
+     */
+    @Benchmark
+    fun rawExpression(blackhole: Blackhole) = invokeAndSum(raw, blackhole)
+
+    /**
+     * Benchmark case for direct computation w/o [Expression].
+     */
+    @Benchmark
+    fun justCalculate(blackhole: Blackhole) {
+        val random = Random(0)
+        var sum = 0.0
+
+        repeat(times) {
+            val x = random.nextDouble()
+            sum += x * 2.0 + 2.0 / x - 16.0 / sin(x)
+        }
+
+        blackhole.consume(sum)
+    }
+
+    private fun invokeAndSum(expr: Expression<Float64>, blackhole: Blackhole) {
+        val random = Random(0)
+        var sum = 0.0
+        val m = HashMap<Symbol, Double>()
+
+        repeat(times) {
+            m[x] = random.nextDouble()
+            sum += expr(m)
+        }
+
+        blackhole.consume(sum)
+    }
+
+    private companion object {
+        private val x by symbol
+        private const val times = 1_000_000
+
+        private val functional = Float64Field.expression {
+            val x = bindSymbol(Symbol.x)
+            x * number(2.0) + 2.0 / x - 16.0 / sin(x)
+        }
+
+        private val node = MstExtendedField {
+            x * 2.0 + number(2.0) / x - number(16.0) / sin(x)
+        }
+
+        private val mst = node.toExpression(Float64Field)
+
+        @OptIn(UnstableKMathAPI::class)
+        private val wasm = node.wasmCompileToExpression(Float64Field)
+        private val estree = node.estreeCompileToExpression(Float64Field)
+
+        private val raw = Expression<Float64> { args ->
+            val x = args.getValue(x)
+            x * 2.0 + 2.0 / x - 16.0 / sin(x)
+        }
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ArrayBenchmark.kt

@@ -0,0 +1,43 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import java.nio.IntBuffer
+
+@State(Scope.Benchmark)
+internal class ArrayBenchmark {
+    @Benchmark
+    fun benchmarkArrayRead(blackhole: Blackhole) {
+        var res = 0
+        for (i in 1..size) res += array[size - i]
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun benchmarkBufferRead(blackhole: Blackhole) {
+        var res = 0
+        for (i in 1..size) res += arrayBuffer[size - i]
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun nativeBufferRead(blackhole: Blackhole) {
+        var res = 0
+        for (i in 1..size) res += nativeBuffer[size - i]
+        blackhole.consume(res)
+    }
+
+    private companion object {
+        private const val size = 1000
+        private val array = IntArray(size) { it }
+        private val arrayBuffer = IntBuffer.wrap(array)
+        private val nativeBuffer = IntBuffer.allocate(size).also { for (i in 0 until size) it.put(i, i) }
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt

@@ -0,0 +1,112 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+
+import kotlinx.benchmark.Blackhole
+import org.openjdk.jmh.annotations.Benchmark
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.operations.BigIntField
+import space.kscience.kmath.operations.JBigIntegerField
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.operations.parseBigInteger
+import java.math.BigInteger
+
+
+@UnstableKMathAPI
+@State(Scope.Benchmark)
+internal class BigIntBenchmark {
+
+    val kmSmallNumber = BigIntField.number(100)
+    val jvmSmallNumber = JBigIntegerField.number(100)
+    val kmNumber = BigIntField.number(Int.MAX_VALUE)
+    val jvmNumber = JBigIntegerField.number(Int.MAX_VALUE)
+    val kmLargeNumber = BigIntField { number(11).pow(100_000U) }
+    val jvmLargeNumber: BigInteger = JBigIntegerField { number(11).pow(100_000) }
+    val bigExponent = 50_000
+
+    @Benchmark
+    fun kmSmallAdd(blackhole: Blackhole) = BigIntField {
+        blackhole.consume(kmSmallNumber + kmSmallNumber + kmSmallNumber)
+    }
+
+    @Benchmark
+    fun jvmSmallAdd(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume(jvmSmallNumber + jvmSmallNumber + jvmSmallNumber)
+    }
+
+    @Benchmark
+    fun kmAdd(blackhole: Blackhole) = BigIntField {
+        blackhole.consume(kmNumber + kmNumber + kmNumber)
+    }
+
+    @Benchmark
+    fun jvmAdd(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume(jvmNumber + jvmNumber + jvmNumber)
+    }
+
+    @Benchmark
+    fun kmAddLarge(blackhole: Blackhole) = BigIntField {
+        blackhole.consume(kmLargeNumber + kmLargeNumber + kmLargeNumber)
+    }
+
+    @Benchmark
+    fun jvmAddLarge(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume(jvmLargeNumber + jvmLargeNumber + jvmLargeNumber)
+    }
+
+    @Benchmark
+    fun kmMultiply(blackhole: Blackhole) = BigIntField {
+        blackhole.consume(kmNumber * kmNumber * kmNumber)
+    }
+
+    @Benchmark
+    fun kmMultiplyLarge(blackhole: Blackhole) = BigIntField {
+        blackhole.consume(kmLargeNumber * kmLargeNumber)
+    }
+
+    @Benchmark
+    fun jvmMultiply(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume(jvmNumber * jvmNumber * jvmNumber)
+    }
+
+    @Benchmark
+    fun jvmMultiplyLarge(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume(jvmLargeNumber * jvmLargeNumber)
+    }
+
+    @Benchmark
+    fun kmPower(blackhole: Blackhole) = BigIntField {
+        blackhole.consume(kmNumber.pow(bigExponent.toUInt()))
+    }
+
+    @Benchmark
+    fun jvmPower(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume(jvmNumber.pow(bigExponent))
+    }
+
+    @Benchmark
+    fun kmParsing16(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume("0x7f57ed8b89c29a3b9a85c7a5b84ca3929c7b7488593".parseBigInteger())
+    }
+
+    @Benchmark
+    fun kmParsing10(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume("236656783929183747565738292847574838922010".parseBigInteger())
+    }
+
+    @Benchmark
+    fun jvmParsing10(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume("236656783929183747565738292847574838922010".toBigInteger(10))
+    }
+
+    @Benchmark
+    fun jvmParsing16(blackhole: Blackhole) = JBigIntegerField {
+        blackhole.consume("7f57ed8b89c29a3b9a85c7a5b84ca3929c7b7488593".toBigInteger(16))
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BufferBenchmark.kt

@@ -0,0 +1,80 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.complex.ComplexField
+import space.kscience.kmath.complex.complex
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64Buffer
+import space.kscience.kmath.structures.getDouble
+import space.kscience.kmath.structures.permute
+
+@State(Scope.Benchmark)
+internal class BufferBenchmark {
+
+    @Benchmark
+    fun doubleArrayReadWrite(blackhole: Blackhole) {
+        val buffer = DoubleArray(size) { it.toDouble() }
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer[it]
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun doubleBufferReadWrite(blackhole: Blackhole) {
+        val buffer = Float64Buffer(size) { it.toDouble() }
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer[it]
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun bufferViewReadWrite(blackhole: Blackhole) {
+        val buffer = Float64Buffer(size) { it.toDouble() }.permute(reversedIndices)
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer[it]
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun bufferViewReadWriteSpecialized(blackhole: Blackhole) {
+        val buffer = Float64Buffer(size) { it.toDouble() }.permute(reversedIndices)
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer.getDouble(it)
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun complexBufferReadWrite(blackhole: Blackhole) = ComplexField {
+        val buffer = Buffer.complex(size / 2) { Complex(it.toDouble(), -it.toDouble()) }
+
+        var res = zero
+        (0 until size / 2).forEach {
+            res += buffer[it]
+        }
+
+        blackhole.consume(res)
+    }
+
+    private companion object {
+        private const val size = 100
+        private val reversedIndices = IntArray(size) { it }.apply { reverse() }
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt

@@ -0,0 +1,93 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.commons.linear.CMLinearSpace
+import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
+import space.kscience.kmath.linear.Float64ParallelLinearSpace
+import space.kscience.kmath.linear.invoke
+import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.tensorflow.produceWithTF
+import space.kscience.kmath.tensors.core.tensorAlgebra
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+internal class DotBenchmark {
+    companion object {
+        val random = Random(12224)
+        const val dim = 1000
+
+        //creating invertible matrix
+        val matrix1 = Float64Field.linearSpace.buildMatrix(dim, dim) { _, _ ->
+            random.nextDouble()
+        }
+        val matrix2 = Float64Field.linearSpace.buildMatrix(dim, dim) { _, _ ->
+            random.nextDouble()
+        }
+
+        val cmMatrix1 = CMLinearSpace { matrix1.toCM() }
+        val cmMatrix2 = CMLinearSpace { matrix2.toCM() }
+
+        val ejmlMatrix1 = EjmlLinearSpaceDDRM { matrix1.toEjml() }
+        val ejmlMatrix2 = EjmlLinearSpaceDDRM { matrix2.toEjml() }
+    }
+
+
+    @Benchmark
+    fun tfDot(blackhole: Blackhole) {
+        blackhole.consume(
+            Float64Field.produceWithTF {
+                matrix1 dot matrix1
+            }
+        )
+    }
+
+    @Benchmark
+    fun cmDotWithConversion(blackhole: Blackhole) = CMLinearSpace {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
+    @Benchmark
+    fun cmDot(blackhole: Blackhole) = CMLinearSpace {
+        blackhole.consume(cmMatrix1 dot cmMatrix2)
+    }
+
+    @Benchmark
+    fun ejmlDot(blackhole: Blackhole) = EjmlLinearSpaceDDRM {
+        blackhole.consume(ejmlMatrix1 dot ejmlMatrix2)
+    }
+
+    @Benchmark
+    fun ejmlDotWithConversion(blackhole: Blackhole) = EjmlLinearSpaceDDRM {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
+    @Benchmark
+    fun multikDot(blackhole: Blackhole) = with(multikAlgebra) {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
+    @Benchmark
+    fun tensorDot(blackhole: Blackhole) = with(Float64Field.tensorAlgebra) {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
+    @Benchmark
+    fun bufferedDot(blackhole: Blackhole) = with(Float64Field.linearSpace) {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
+    @Benchmark
+    fun parallelDot(blackhole: Blackhole) = with(Float64ParallelLinearSpace) {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt

@@ -0,0 +1,125 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.asm.compileToExpression
+import space.kscience.kmath.expressions.*
+import space.kscience.kmath.operations.Algebra
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.bindSymbol
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Float64
+import kotlin.math.sin
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+internal class ExpressionsInterpretersBenchmark {
+    /**
+     * Benchmark case for [Expression] created with [expressionInExtendedField].
+     */
+    @Benchmark
+    fun functionalExpression(blackhole: Blackhole) = invokeAndSum(functional, blackhole)
+
+    /**
+     * Benchmark case for [Expression] created with [toExpression].
+     */
+    @Benchmark
+    fun mstExpression(blackhole: Blackhole) = invokeAndSum(mst, blackhole)
+
+    /**
+     * Benchmark case for [Expression] created with [compileToExpression].
+     */
+    @Benchmark
+    fun asmGenericExpression(blackhole: Blackhole) = invokeAndSum(asmGeneric, blackhole)
+
+    /**
+     * Benchmark case for [Expression] created with [compileToExpression].
+     */
+    @Benchmark
+    fun asmPrimitiveExpressionArray(blackhole: Blackhole) {
+        val random = Random(0)
+        var sum = 0.0
+        val m = DoubleArray(1)
+
+        repeat(times) {
+            m[xIdx] = random.nextDouble()
+            sum += asmPrimitive(m)
+        }
+
+        blackhole.consume(sum)
+    }
+
+    /**
+     * Benchmark case for [Expression] created with [compileToExpression].
+     */
+    @Benchmark
+    fun asmPrimitiveExpression(blackhole: Blackhole) = invokeAndSum(asmPrimitive, blackhole)
+
+    /**
+     * Benchmark case for [Expression] implemented manually with `kotlin.math` functions.
+     */
+    @Benchmark
+    fun rawExpression(blackhole: Blackhole) = invokeAndSum(raw, blackhole)
+
+    /**
+     * Benchmark case for direct computation w/o [Expression].
+     */
+    @Benchmark
+    fun justCalculate(blackhole: Blackhole) {
+        val random = Random(0)
+        var sum = 0.0
+
+        repeat(times) {
+            val x = random.nextDouble()
+            sum += x * 2.0 + 2.0 / x - 16.0 / sin(x)
+        }
+
+        blackhole.consume(sum)
+    }
+
+    private fun invokeAndSum(expr: Expression<Float64>, blackhole: Blackhole) {
+        val random = Random(0)
+        var sum = 0.0
+        val m = HashMap<Symbol, Double>()
+
+        repeat(times) {
+            m[x] = random.nextDouble()
+            sum += expr(m)
+        }
+
+        blackhole.consume(sum)
+    }
+
+    private companion object {
+        private val x by symbol
+        private const val times = 1_000_000
+
+        private val functional = Float64Field.expression {
+            val x = bindSymbol(Symbol.x)
+            x * number(2.0) + 2.0 / x - 16.0 / sin(x)
+        }
+
+        private val node = MstExtendedField {
+            x * 2.0 + number(2.0) / x - number(16.0) / sin(x)
+        }
+
+        private val mst = node.toExpression(Float64Field)
+
+        private val asmPrimitive = node.compileToExpression(Float64Field)
+        private val xIdx = asmPrimitive.indexer.indexOf(x)
+
+        private val asmGeneric = node.compileToExpression(Float64Field as Algebra<Float64>)
+
+        private val raw = Expression<Float64> { args ->
+            val x = args[x]!!
+            x * 2.0 + 2.0 / x - 16.0 / sin(x)
+        }
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/IntegrationBenchmark.kt

@@ -0,0 +1,40 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import org.openjdk.jmh.annotations.Benchmark
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import org.openjdk.jmh.infra.Blackhole
+import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.complex.algebra
+import space.kscience.kmath.integration.gaussIntegrator
+import space.kscience.kmath.integration.integrate
+import space.kscience.kmath.integration.value
+import space.kscience.kmath.operations.algebra
+
+
+@State(Scope.Benchmark)
+internal class IntegrationBenchmark {
+
+    @Benchmark
+    fun doubleIntegration(blackhole: Blackhole) {
+        val res = Double.algebra.gaussIntegrator.integrate(0.0..1.0, intervals = 1000) { x: Double ->
+            //sin(1 / x)
+            1 / x
+        }.value
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun complexIntegration(blackhole: Blackhole) = with(Complex.algebra) {
+        val res = gaussIntegrator.integrate(0.0..1.0, intervals = 1000) { x: Double ->
+//            sin(1 / x) + i * cos(1 / x)
+            1 / x - i / x
+        }.value
+        blackhole.consume(res)
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/JafamaBenchmark.kt

@@ -0,0 +1,39 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Blackhole
+import org.openjdk.jmh.annotations.Benchmark
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import space.kscience.kmath.jafama.JafamaDoubleField
+import space.kscience.kmath.jafama.StrictJafamaDoubleField
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.invoke
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+internal class JafamaBenchmark {
+    @Benchmark
+    fun jafama(blackhole: Blackhole) = invokeBenchmarks(blackhole) { x ->
+        JafamaDoubleField { x * power(x, 4) * exp(x) / cos(x) + sin(x) }
+    }
+
+    @Benchmark
+    fun core(blackhole: Blackhole) = invokeBenchmarks(blackhole) { x ->
+        Float64Field { x * power(x, 4) * exp(x) / cos(x) + sin(x) }
+    }
+
+    @Benchmark
+    fun strictJafama(blackhole: Blackhole) = invokeBenchmarks(blackhole) { x ->
+        StrictJafamaDoubleField { x * power(x, 4) * exp(x) / cos(x) + sin(x) }
+    }
+}
+
+private inline fun invokeBenchmarks(blackhole: Blackhole, expr: (Double) -> Double) {
+    val rng = Random(0)
+    repeat(1000000) { blackhole.consume(expr(rng.nextDouble())) }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt

@@ -0,0 +1,57 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.commons.linear.CMLinearSpace
+import space.kscience.kmath.commons.linear.lupSolver
+import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
+import space.kscience.kmath.linear.invoke
+import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.linear.lupSolver
+import space.kscience.kmath.linear.parallel
+import space.kscience.kmath.operations.algebra
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+internal class MatrixInverseBenchmark {
+    private companion object {
+        private val random = Random(1224)
+        private const val dim = 100
+
+        private val space = Double.algebra.linearSpace
+
+        //creating invertible matrix
+        private val u = space.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        private val l = space.buildMatrix(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
+        private val matrix = space { l dot u }
+    }
+
+    @Benchmark
+    fun kmathLupInversion(blackhole: Blackhole) {
+        blackhole.consume(Double.algebra.linearSpace.lupSolver().inverse(matrix))
+    }
+
+    @Benchmark
+    fun kmathParallelLupInversion(blackhole: Blackhole) {
+        blackhole.consume(Double.algebra.linearSpace.parallel.lupSolver().inverse(matrix))
+    }
+
+    @Benchmark
+    fun cmLUPInversion(blackhole: Blackhole) = CMLinearSpace {
+        blackhole.consume(lupSolver().inverse(matrix))
+    }
+
+
+    @Benchmark
+    fun ejmlInverse(blackhole: Blackhole) = EjmlLinearSpaceDDRM {
+        blackhole.consume(matrix.toEjml().inverted())
+    }
+
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt

@@ -0,0 +1,97 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import org.jetbrains.kotlinx.multik.api.Multik
+import org.jetbrains.kotlinx.multik.api.ones
+import org.jetbrains.kotlinx.multik.ndarray.data.DN
+import org.jetbrains.kotlinx.multik.ndarray.data.DataType
+import space.kscience.kmath.UnsafeKMathAPI
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.nd4j.nd4j
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.tensors.core.DoubleTensor
+import space.kscience.kmath.tensors.core.one
+import space.kscience.kmath.tensors.core.tensorAlgebra
+import space.kscience.kmath.viktor.viktorAlgebra
+
+@State(Scope.Benchmark)
+internal class NDFieldBenchmark {
+
+    private companion object {
+        private const val dim = 1000
+        private const val n = 100
+        private val shape = ShapeND(dim, dim)
+        private val specializedField = Float64Field.ndAlgebra
+        private val genericField = BufferedFieldOpsND(Float64Field)
+        private val nd4jField = Float64Field.nd4j
+        private val viktorField = Float64Field.viktorAlgebra
+    }
+
+    @Benchmark
+    fun specializedFieldAdd(blackhole: Blackhole) = with(specializedField) {
+        var res: StructureND<Float64> = one(shape)
+        repeat(n) { res += 1.0 }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun boxingFieldAdd(blackhole: Blackhole) = with(genericField) {
+        var res: StructureND<Float64> = one(shape)
+        repeat(n) { res += 1.0 }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun multikAdd(blackhole: Blackhole) = with(multikAlgebra) {
+        var res: StructureND<Float64> = one(shape)
+        repeat(n) { res += 1.0 }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun viktorAdd(blackhole: Blackhole) = with(viktorField) {
+        var res: StructureND<Float64> = one(shape)
+        repeat(n) { res += 1.0 }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun tensorAdd(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+        var res: DoubleTensor = one(shape)
+        repeat(n) { res = res + 1.0 }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun tensorInPlaceAdd(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+        val res: DoubleTensor = one(shape)
+        repeat(n) { res += 1.0 }
+        blackhole.consume(res)
+    }
+
+    @OptIn(UnsafeKMathAPI::class)
+    @Benchmark
+    fun multikInPlaceAdd(blackhole: Blackhole) = with(multikAlgebra) {
+        val res = Multik.ones<Double, DN>(shape.asArray(), DataType.DoubleDataType).wrap()
+        repeat(n) { res += 1.0 }
+        blackhole.consume(res)
+    }
+
+//    @Benchmark
+//    fun nd4jAdd(blackhole: Blackhole) = with(nd4jField) {
+//        var res: StructureND<Float64> = one(dim, dim)
+//        repeat(n) { res += 1.0 }
+//        blackhole.consume(res)
+//    }
+
+
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt

@@ -0,0 +1,39 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.linear.matrix
+import space.kscience.kmath.linear.symmetric
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.tensors.core.symEigJacobi
+import space.kscience.kmath.tensors.core.symEigSvd
+import space.kscience.kmath.tensors.core.tensorAlgebra
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+internal class TensorAlgebraBenchmark {
+    companion object {
+        private val random = Random(12224)
+        private const val dim = 30
+
+        private val matrix = Float64Field.linearSpace.matrix(dim, dim).symmetric { _, _ -> random.nextDouble() }
+    }
+
+    @Benchmark
+    fun tensorSymEigSvd(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+        blackhole.consume(symEigSvd(matrix, 1e-10))
+    }
+
+    @Benchmark
+    fun tensorSymEigJacobi(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+        blackhole.consume(symEigJacobi(matrix, 50, 1e-10))
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt

@@ -0,0 +1,59 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import org.jetbrains.bio.viktor.F64Array
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.ndAlgebra
+import space.kscience.kmath.nd.one
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.viktor.ViktorFieldND
+
+@State(Scope.Benchmark)
+internal class ViktorBenchmark {
+
+    @Benchmark
+    fun doubleFieldAddition(blackhole: Blackhole) {
+        with(doubleField) {
+            var res: StructureND<Float64> = one(shape)
+            repeat(n) { res += 1.0 }
+            blackhole.consume(res)
+        }
+    }
+
+    @Benchmark
+    fun viktorFieldAddition(blackhole: Blackhole) {
+        with(viktorField) {
+            var res = one(shape)
+            repeat(n) { res += 1.0 }
+            blackhole.consume(res)
+        }
+    }
+
+    @Benchmark
+    fun rawViktor(blackhole: Blackhole) {
+        val one = F64Array.full(init = 1.0, shape = intArrayOf(dim, dim))
+        var res = one
+        repeat(n) { res = res + one }
+        blackhole.consume(res)
+    }
+
+    private companion object {
+        private const val dim = 1000
+        private const val n = 100
+        private val shape = ShapeND(dim, dim)
+
+        // automatically build context most suited for given type.
+        private val doubleField = Float64Field.ndAlgebra
+        private val viktorField = ViktorFieldND(dim, dim)
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt

@@ -0,0 +1,58 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
+import kotlinx.benchmark.Scope
+import kotlinx.benchmark.State
+import org.jetbrains.bio.viktor.F64Array
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.ndAlgebra
+import space.kscience.kmath.nd.one
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.viktor.ViktorFieldND
+
+@State(Scope.Benchmark)
+internal class ViktorLogBenchmark {
+    @Benchmark
+    fun realFieldLog(blackhole: Blackhole) {
+        with(doubleField) {
+            val fortyTwo = structureND(shape) { 42.0 }
+            var res = one(shape)
+            repeat(n) { res = ln(fortyTwo) }
+            blackhole.consume(res)
+        }
+    }
+
+    @Benchmark
+    fun viktorFieldLog(blackhole: Blackhole) {
+        with(viktorField) {
+            val fortyTwo = structureND(shape) { 42.0 }
+            var res = one
+            repeat(n) { res = ln(fortyTwo) }
+            blackhole.consume(res)
+        }
+    }
+
+    @Benchmark
+    fun rawViktorLog(blackhole: Blackhole) {
+        val fortyTwo = F64Array.full(dim, dim, init = 42.0)
+        lateinit var res: F64Array
+        repeat(n) { res = fortyTwo.log() }
+        blackhole.consume(res)
+    }
+
+    private companion object {
+        private const val dim = 1000
+        private const val n = 100
+        private val shape = ShapeND(dim, dim)
+
+        // automatically build context most suited for given type.
+        private val doubleField = Float64Field.ndAlgebra
+        private val viktorField = ViktorFieldND(dim, dim)
+    }
+}

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/globals.kt

@@ -0,0 +1,11 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.benchmarks
+
+import org.jetbrains.kotlinx.multik.default.DefaultEngine
+import space.kscience.kmath.multik.MultikDoubleAlgebra
+
+val multikAlgebra = MultikDoubleAlgebra(DefaultEngine())

build.gradle.kts

@@ -1,25 +1,72 @@
+import space.kscience.gradle.useApache2Licence
+import space.kscience.gradle.useSPCTeam
+
 plugins {
-    id("scientifik.publish") apply false
+    alias(spclibs.plugins.kscience.project)
+    alias(spclibs.plugins.kotlinx.kover)
 }
 
-val kmathVersion by extra("0.1.4-dev-8")
-
-val bintrayRepo by extra("scientifik")
-val githubProject by extra("kmath")
+val attributesVersion by extra("0.2.0")
 
 allprojects {
     repositories {
-        jcenter()
-        maven("https://dl.bintray.com/kotlin/kotlinx")
-        maven("https://dl.bintray.com/hotkeytlt/maven")
+        maven("https://repo.kotlin.link")
+        maven("https://oss.sonatype.org/content/repositories/snapshots")
+        mavenCentral()
     }
 
-    group = "scientifik"
-    version = kmathVersion
+    group = "space.kscience"
+    version = "0.4.1-dev"
 }
 
 subprojects {
-    if (name.startsWith("kmath")) {
-        apply(plugin = "scientifik.publish")
+    if (name.startsWith("kmath")) apply<MavenPublishPlugin>()
+
+    plugins.withId("org.jetbrains.dokka") {
+        tasks.withType<org.jetbrains.dokka.gradle.DokkaTaskPartial> {
+            dependsOn(tasks["assemble"])
+
+            dokkaSourceSets.all {
+                val readmeFile = this@subprojects.projectDir.resolve("README.md")
+                if (readmeFile.exists()) includes.from(readmeFile)
+                val kotlinDirPath = "src/$name/kotlin"
+                val kotlinDir = file(kotlinDirPath)
+
+                if (kotlinDir.exists()) sourceLink {
+                    localDirectory.set(kotlinDir)
+
+                    remoteUrl.set(
+                        uri("https://github.com/SciProgCentre/kmath/tree/master/${this@subprojects.name}/$kotlinDirPath").toURL()
+                    )
+                }
+
+                externalDocumentationLink("https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/")
+                externalDocumentationLink("https://deeplearning4j.org/api/latest/")
+                externalDocumentationLink("https://axelclk.bitbucket.io/symja/javadoc/")
+
+                externalDocumentationLink(
+                    "https://kotlin.github.io/kotlinx.coroutines/kotlinx-coroutines-core/",
+                    "https://kotlin.github.io/kotlinx.coroutines/package-list",
+                )
+
+                externalDocumentationLink(
+                    "https://breandan.net/kotlingrad/kotlingrad/",
+                    "https://breandan.net/kotlingrad/kotlingrad/kotlingrad/package-list",
+                )
+            }
+        }
     }
-}
+}
+
+readme.readmeTemplate = file("docs/templates/README-TEMPLATE.md")
+
+ksciencePublish {
+    pom("https://github.com/SciProgCentre/kmath") {
+        useApache2Licence()
+        useSPCTeam()
+    }
+    repository("spc", "https://maven.sciprog.center/kscience")
+    sonatype("https://oss.sonatype.org")
+}
+
+apiValidation.nonPublicMarkers.add("space.kscience.kmath.UnstableKMathAPI")

doc/algebra.md

@@ -1,125 +0,0 @@
-# Algebraic Structures and Algebraic Elements
-
-The mathematical operations in KMath are generally separated from mathematical objects. This means that to perform an 
-operation, say `+`, one needs two objects of a type `T` and an algebra context, which draws appropriate operation up, 
-say `Space<T>`. Next one needs to run the actual operation in the context:
-
-```kotlin
-import scientifik.kmath.operations.*
-
-val a: T = ...
-val b: T = ...
-val space: Space<T> = ...
-
-val c = space { a + b }
-```
-
-At first glance, this distinction seems to be a needless complication, but in fact one needs to remember that in 
-mathematics, one could draw up different operations on same objects. For example, one could use different types of 
-geometry for vectors.
-
-## Algebraic Structures
-
-Mathematical contexts have the following hierarchy:
-
-**Algebra** ← **Space** ← **Ring** ← **Field**
-
-These interfaces follow real algebraic structures:
-
-- [Space](https://mathworld.wolfram.com/VectorSpace.html) defines addition, its neutral element (i.e. 0) and scalar 
-multiplication;
-- [Ring](http://mathworld.wolfram.com/Ring.html) adds multiplication and its neutral element (i.e. 1);
-- [Field](http://mathworld.wolfram.com/Field.html) adds division operation.
-
-A typical implementation of `Field<T>` is the `RealField` which works on doubles, and `VectorSpace` for `Space<T>`.
-
-In some cases algebra context can hold additional operations like `exp` or `sin`, and then it inherits appropriate
-interface. Also, contexts may have operations, which produce elements outside of the context. For example, `Matrix.dot` 
-operation produces a matrix with new dimensions, which can be incompatible with initial matrix in terms of linear 
-operations.
-
-## Algebraic Element
-
-To achieve more familiar behavior (where you apply operations directly to mathematical objects), without involving 
-contexts KMath submits special type objects called `MathElement`. A `MathElement` is basically some object coupled to
-a mathematical context. For example `Complex` is the pair of real numbers representing real and imaginary parts,
-but it also holds reference to the `ComplexField` singleton, which allows performing direct operations on `Complex`
-numbers without explicit involving the context like:
-
-```kotlin
-import scientifik.kmath.operations.*
-
-// Using elements
-val c1 = Complex(1.0, 1.0)
-val c2 = Complex(1.0, -1.0)
-val c3 = c1 + c2 + 3.0.toComplex()
-
-// Using context
-val c4 = ComplexField { c1 + i - 2.0 }
-```
-
-Both notations have their pros and cons.
-
-The hierarchy for algebraic elements follows the hierarchy for the corresponding algebraic structures.
-
-**MathElement** ← **SpaceElement** ← **RingElement** ← **FieldElement**
-
-`MathElement<C>` is the generic common ancestor of the class with context.
-
-One major distinction between algebraic elements and algebraic contexts is that elements have three type 
-parameters:
-
-1. The type of elements, the field operates on.
-2. The self-type of the element returned from operation (which has to be an algebraic element).
-3. The type of the algebra over first type-parameter.
-
-The middle type is needed for of algebra members do not store context. For example, it is impossible to add a context 
-to regular `Double`. The element performs automatic conversions from context types and back. One should use context 
-operations in all performance-critical places. The performance of element operations is not guaranteed.
-
-## Spaces and Fields
-
-KMath submits both contexts and elements for builtin algebraic structures:
-
-```kotlin
-import scientifik.kmath.operations.*
-
-val c1 = Complex(1.0, 2.0)
-val c2 = ComplexField.i
-
-val c3 = c1 + c2
-// or
-val c3 = ComplexField { c1 + c2 }
-```
-
-Also, `ComplexField` features special operations to mix complex and real numbers, for example:
-
-```kotlin
-import scientifik.kmath.operations.*
-
-val c1 = Complex(1.0, 2.0)
-val c2 = ComplexField { c1 - 1.0 } // Returns: Complex(re=0.0, im=2.0)
-val c3 = ComplexField { c1 - i * 2.0 }
-```
-
-**Note**: In theory it is possible to add behaviors directly to the context, but as for now Kotlin does not support 
-that. Watch [KT-10468](https://youtrack.jetbrains.com/issue/KT-10468) and 
-[KEEP-176](https://github.com/Kotlin/KEEP/pull/176) for updates.
-
-## Nested fields
-
-Contexts allow one to build more complex structures. For example, it is possible to create a `Matrix` from complex 
-elements like so:
-
-```kotlin
-val element = NDElement.complex(shape = intArrayOf(2, 2)) { index: IntArray ->
-    Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble())
-}
-```
-
-The `element` in this example is a member of the `Field` of 2D structures, each element of which is a member of its own
-`ComplexField`. It is important one does not need to create a special n-d class to hold complex
-numbers and implement operations on it, one just needs to provide a field for its elements.
-
-**Note**: Fields themselves do not solve the problem of JVM boxing, but it is possible to solve with special contexts like
-`MemorySpec`.

doc/contexts.md

@@ -1,73 +0,0 @@
-# Context-oriented mathematics
-
-## The problem
-
-A known problem for implementing mathematics in statically-typed languages (but not only in them) is that different
-sets of mathematical operators can be defined on the same mathematical objects. Sometimes there is no single way to
-treat some operations, including basic arithmetic operations, on a Java/Kotlin `Number`. Sometimes there are different ways to
-define the same structure, such as Euclidean and elliptic geometry vector spaces over real vectors. Another problem arises when
-one wants to add some kind of behavior to an existing entity. In dynamic languages those problems are usually solved
-by adding dynamic context-specific behaviors at runtime, but this solution has a lot of drawbacks.
-
-## Context-oriented approach
-
-One possible solution to these problems is to divorce numerical representations from behaviors.
-For example in Kotlin one can define a separate class which represents some entity without any operations,
-ex. a complex number:
-
-```kotlin
-data class Complex(val re: Double, val im: Double)
-```
-
-And then to define a separate class or singleton, representing an operation on those complex numbers:
-
-```kotlin
-object ComplexOperations {
-    operator fun Complex.plus(other: Complex) = Complex(re + other.re, im + other.im)
-    operator fun Complex.minus(other: Complex) = Complex(re - other.re, im - other.im)
-}
-```
-
-In Java, applying such external operations could be very cumbersome, but Kotlin has a unique feature which allows us
-implement this naturally: [extensions with receivers](https://kotlinlang.org/docs/reference/extensions.html#extension-functions).
-In Kotlin, an operation on complex number could be implemented as:
-
-```kotlin
-with(ComplexOperations) { c1 + c2 - c3 }
-```
-
-Kotlin also allows the creation of functions with receivers:
-
-```kotlin
-fun ComplexOperations.doSomethingWithComplex(c1: Complex, c2: Complex, c3: Complex) = c1 + c2 - c3
-
-ComplexOperations.doComethingWithComplex(c1, c2, c3)
-```
-
-In fact, whole parts of a program may be run within a mathematical context or even multiple nested contexts.
-
-In KMath, contexts are not only responsible for operations, but also for raw object creation and advanced features.
-
-## Other possibilities
-
-### Type classes
-
-An obvious candidate to get more or less the same functionality is the type class, which allows one to bind a behavior to
-a specific type without modifying the type itself. On the plus side, type classes do not require explicit context
-declaration, so the code looks cleaner. On the minus side, if there are different sets of behaviors for the same types,
-it is impossible to combine them into one module. Also, unlike type classes, context can have parameters or even
-state. For example in KMath, sizes and strides for `NDElement` or `Matrix` could be moved to context to optimize
-performance in case of a large amount of structures.
-
-### Wildcard imports and importing-on-demand
-
-Sometimes, one may wish to use a single context throughout a file. In this case, is possible to import all members
-from a package or file, via `import context.complex.*`. Effectively, this is the same as enclosing an entire file
-with a single context. However when using multiple contexts, this technique can introduce operator ambiguity, due to
-namespace pollution. If there are multiple scoped contexts which define the same operation, it is still possible to
-to import specific operations as needed, without using an explicit context with extension functions, for example:
-
-```
-import context.complex.op1
-import context.quaternion.op2
-```

doc/expressions.md

@@ -1,26 +0,0 @@
-# Expressions
-
-**Experimental: this API is in early stage and could change any time**
-
-Expressions is an experimental feature which allows to construct lazily or immediately calculated parametric mathematical 
-expressions.
-
-The potential use-cases for it (so far) are following:
-
-* Lazy evaluation (in general simple lambda is better, but there are some border cases)
-
-* Automatic differentiation in single-dimension and in multiple dimensions
-
-* Generation of mathematical syntax trees with subsequent code generation for other languages
-
-* Maybe symbolic computations (needs additional research)
-
-The workhorse of this API is `Expression` interface which exposes single `operator fun invoke(arguments: Map<String, T>): T`
-method. `ExpressionContext` is used to generate expressions and introduce variables.
-
-Currently there are two implementations:
-
-* Generic `ExpressionField` in `kmath-core` which allows construction of custom lazy expressions
-
-* Auto-differentiation expression in `kmath-commons` module allows to use full power of `DerivativeStructure` 
-from commons-math. **TODO: add example**

doc/features.md

@@ -1,17 +0,0 @@
-# Features
-
-* [Algebra](./algebra.md) - [Context-based](./contexts.md) operations on different primitives and structures.
-
-* [NDStructures](./nd-structure.md)
-
-* [Linear algebra](./linear.md) - Matrices, operations and linear equations solving. To be moved to separate module. Currently supports basic
-api and multiple library back-ends.
-
-* [Histograms](./histograms.md) - Multidimensional histogram calculation and operations.
-
-* [Expressions](./expressions.md)
-
-* Commons math integration
-
-* Koma integration
-

doc/linear.md

@@ -1,19 +0,0 @@
-## Basic linear algebra layout
-
-KMath support for linear algebra organized in a context-oriented way. Meaning that operations are in most cases declared 
-in context classes, and are not the members of classes that store data. This allows more flexible approach to maintain multiple 
-back-ends. The new operations added as extensions to contexts instead of being member functions of data structures.
-
-Two major contexts used for linear algebra and hyper-geometry:
-
-* `VectorSpace` forms a mathematical space on top of array-like structure (`Buffer` and its typealias `Point` used for geometry). 
-
-* `MatrixContext` forms a space-like context for 2d-structures. It does not store matrix size and therefore does not implement
-`Space` interface (it is not possible to create zero element without knowing the matrix size). 
-
-## Vector spaces
-
-
-## Matrix operations
-
-## Back-end overview

docs/algebra.md

@@ -0,0 +1,86 @@
+# Algebraic Structures and Algebraic Elements
+
+The mathematical operations in KMath are generally separated from mathematical objects. This means that to perform an
+operation, say `+`, one needs two objects of a type `T` and an algebra context, which draws appropriate operation up,
+say `Group<T>`. Next one needs to run the actual operation in the context:
+
+```kotlin
+import space.kscience.kmath.operations.*
+
+val a: T = ...
+val b: T = ...
+val group: Group<T> = ...
+
+val c = group { a + b }
+```
+
+At first glance, this distinction seems to be a needless complication, but in fact one needs to remember that in
+mathematics, one could draw up different operations on same objects. For example, one could use different types of
+geometry for vectors.
+
+## Algebraic Structures
+
+Primary mathematical contexts have the following hierarchy:
+
+`Field <: Ring <: Group <: Algebra`
+
+These interfaces follow real algebraic structures:
+
+- [Group](https://mathworld.wolfram.com/Group.html) defines addition, its identity element (i.e., 0) and additive
+  inverse (-x);
+- [Ring](http://mathworld.wolfram.com/Ring.html) adds multiplication and its identity element (i.e., 1);
+- [Field](http://mathworld.wolfram.com/Field.html) adds division operation.
+
+A typical implementation of `Field<T>` is the `DoubleField` which works on doubles, and `VectorSpace` for `Space<T>`.
+
+In some cases algebra context can hold additional operations like `exp` or `sin`, and then it inherits appropriate
+interface. Also, contexts may have operations, which produce elements outside the context. For example, `Matrix.dot`
+operation produces a matrix with new dimensions, which can be incompatible with initial matrix in linear operations.
+
+## Spaces and Fields
+
+KMath introduces contexts for builtin algebraic structures:
+
+```kotlin
+import space.kscience.kmath.operations.*
+
+val c1 = Complex(1.0, 2.0)
+val c2 = ComplexField.i
+
+val c3 = c1 + c2
+// or
+val c3 = ComplexField { c1 + c2 }
+```
+
+Also, `ComplexField` features special operations to mix complex and real numbers, for example:
+
+```kotlin
+import space.kscience.kmath.operations.*
+
+val c1 = Complex(1.0, 2.0)
+val c2 = ComplexField { c1 - 1.0 } // Returns: Complex(re=0.0, im=2.0)
+val c3 = ComplexField { c1 - i * 2.0 }
+```
+
+**Note**: In theory it is possible to add behaviors directly to the context, but as for now Kotlin does not support
+that. Watch [KT-10468](https://youtrack.jetbrains.com/issue/KT-10468) and
+[KEEP-176](https://github.com/Kotlin/KEEP/pull/176) for updates.
+
+## Nested fields
+
+Contexts allow one to build more complex structures. For example, it is possible to create a `Matrix` from complex
+elements like so:
+
+```kotlin
+val element = NDElement.complex(shape = intArrayOf(2, 2)) { index: IntArray ->
+    Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble())
+}
+```
+
+The `element` in this example is a member of the `Field` of 2D structures, each element of which is a member of its own
+`ComplexField`. It is important one does not need to create a special n-d class to hold complex numbers and implement
+operations on it, one just needs to provide a field for its elements.
+
+**Note**: Fields themselves do not solve the problem of JVM boxing, but it is possible to solve with special contexts
+like
+`MemorySpec`.

docs/buffers.md

@@ -1,17 +1,20 @@
 # Buffers
 
-Buffer is one of main building blocks of kmath. It is a basic interface allowing random-access read and write (with `MutableBuffer`).
-There are different types of buffers:
+Buffer is one of main building blocks of kmath. It is a basic interface allowing random-access read and write (
+with `MutableBuffer`). There are different types of buffers:
 
-* Primitive buffers wrapping like `RealBuffer` which are wrapping primitive arrays.
+* Primitive buffers wrapping like `DoubleBuffer` which are wrapping primitive arrays.
 * Boxing `ListBuffer` wrapping a list
 * Functionally defined `VirtualBuffer` which does not hold a state itself, but provides a function to calculate value
 * `MemoryBuffer` allows direct allocation of objects in continuous memory block.
 
-Some kmath features require a `BufferFactory` class to operate properly. A general convention is to use functions defined in
-`Buffer` and `MutableBuffer` companion classes. For example factory `Buffer.Companion::auto` in most cases creates the most suitable
-buffer for given reified type (for types with custom memory buffer it still better to use their own `MemoryBuffer.create()` factory).
+Some kmath features require a `BufferFactory` class to operate properly. A general convention is to use functions
+defined in
+`Buffer` and `MutableBuffer` companion classes. For example factory `Buffer.Companion::auto` in most cases creates the
+most suitable buffer for given reified type (for types with custom memory buffer it still better to use their
+own `MemoryBuffer.create()` factory).
 
 ## Buffer performance
 
-One should avoid using default boxing buffer wherever it is possible. Try to use primitive buffers or memory buffers instead  
+One should avoid using default boxing buffer wherever it is possible. Try to use primitive buffers or memory buffers
+instead.

docs/codestyle.md

@@ -1,34 +1,35 @@
 # Coding Conventions
 
-KMath code follows general [Kotlin conventions](https://kotlinlang.org/docs/reference/coding-conventions.html), but 
-with a number of small changes and clarifications.
+Generally, KMath code follows
+general [Kotlin coding conventions](https://kotlinlang.org/docs/reference/coding-conventions.html), but with a number of
+small changes and clarifications.
 
 ## Utility Class Naming
 
-Filename should coincide with a name of one of the classes contained in the file or start with small letter and 
-describe its contents.
+Filename should coincide with a name of one of the classes contained in the file or start with small letter and describe
+its contents.
 
-The code convention [here](https://kotlinlang.org/docs/reference/coding-conventions.html#source-file-names) says that 
-file names should start with a capital letter even if file does not contain classes. Yet starting utility classes and 
+The code convention [here](https://kotlinlang.org/docs/reference/coding-conventions.html#source-file-names) says that
+file names should start with a capital letter even if file does not contain classes. Yet starting utility classes and
 aggregators with a small letter seems to be a good way to visually separate those files.
 
 This convention could be changed in future in a non-breaking way.
 
 ## Private Variable Naming
 
-Private variables' names may start with underscore `_` for of the private mutable variable is shadowed by the public 
+Private variables' names may start with underscore `_` for of the private mutable variable is shadowed by the public
 read-only value with the same meaning.
 
-This rule does not permit underscores in names, but it is sometimes useful to "underscore" the fact that public and 
+This rule does not permit underscores in names, but it is sometimes useful to "underscore" the fact that public and
 private versions draw up the same entity. It is allowed only for private variables.
 
 This convention could be changed in future in a non-breaking way.
 
 ## Functions and Properties One-liners
 
-Use one-liners when they occupy single code window line both for functions and properties with getters like 
-`val b: String get() = "fff"`. The same should be performed with multiline expressions when they could be 
+Use one-liners when they occupy single code window line both for functions and properties with getters like
+`val b: String get() = "fff"`. The same should be performed with multiline expressions when they could be
 cleanly separated.
 
-There is no universal consensus whenever use `fun a() = ...` or `fun a() { return ... }`. Yet from reader outlook 
+There is no universal consensus whenever use `fun a() = ...` or `fun a() { return ... }`. Yet from reader outlook
 one-lines seem to better show that the property or function is easily calculated.

docs/contexts.md

@@ -0,0 +1,73 @@
+# Context-oriented mathematics
+
+## The problem
+
+A known problem for implementing mathematics in statically-typed languages (but not only in them) is that different sets
+of mathematical operators can be defined on the same mathematical objects. Sometimes there is no single way to treat
+some operations, including basic arithmetic operations, on a Java/Kotlin `Number`. Sometimes there are different ways to
+define the same structure, such as Euclidean and elliptic geometry vector spaces over real vectors. Another problem
+arises when one wants to add some kind of behavior to an existing entity. In dynamic languages those problems are
+usually solved by adding dynamic context-specific behaviors at runtime, but this solution has a lot of drawbacks.
+
+## Context-oriented approach
+
+One possible solution to these problems is to divorce numerical representations from behaviors. For example in Kotlin
+one can define a separate class representing some entity without any operations, ex. a complex number:
+
+```kotlin
+data class Complex(val re: Double, val im: Double)
+```
+
+And then to define a separate class or singleton, representing an operation on those complex numbers:
+
+```kotlin
+object ComplexOperations {
+    operator fun Complex.plus(other: Complex) = Complex(re + other.re, im + other.im)
+    operator fun Complex.minus(other: Complex) = Complex(re - other.re, im - other.im)
+}
+```
+
+In Java, applying such external operations could be cumbersome, but Kotlin has a unique feature that allows us
+implement this
+naturally: [extensions with receivers](https://kotlinlang.org/docs/reference/extensions.html#extension-functions). In
+Kotlin, an operation on complex number could be implemented as:
+
+```kotlin
+with(ComplexOperations) { c1 + c2 - c3 }
+```
+
+Kotlin also allows the creation of functions with receivers:
+
+```kotlin
+fun ComplexOperations.doSomethingWithComplex(c1: Complex, c2: Complex, c3: Complex) = c1 + c2 - c3
+
+ComplexOperations.doComethingWithComplex(c1, c2, c3)
+```
+
+In fact, whole parts of a program may be run within a mathematical context or even multiple nested contexts.
+
+In KMath, contexts are not only responsible for operations, but also for raw object creation and advanced features.
+
+## Other possibilities
+
+### Type classes
+
+An obvious candidate to get more or less the same functionality is the type class, which allows one to bind a behavior
+to a specific type without modifying the type itself. On the plus side, type classes do not require explicit context
+declaration, so the code looks cleaner. On the minus side, if there are different sets of behaviors for the same types,
+it is impossible to combine them into one module. Also, unlike type classes, context can have parameters or even state.
+For example in KMath, sizes and strides for `NDElement` or `Matrix` could be moved to context to optimize performance in
+case of a large amount of structures.
+
+### Wildcard imports and importing-on-demand
+
+Sometimes, one may wish to use a single context throughout a file. In this case, is possible to import all members from
+a package or file, via `import context.complex.*`. Effectively, this is the same as enclosing an entire file with a
+single context. However, when using multiple contexts, this technique can introduce operator ambiguity, due to namespace
+pollution. If there are multiple scoped contexts that define the same operation, it is still possible to import
+specific operations as needed, without using an explicit context with extension functions, for example:
+
+```
+import context.complex.op1
+import context.quaternion.op2
+```

docs/expressions.md

@@ -0,0 +1,24 @@
+# Expressions
+
+Expressions is a feature, which allows constructing lazily or immediately calculated parametric mathematical
+expressions.
+
+The potential use-cases for it (so far) are following:
+
+* lazy evaluation (in general simple lambda is better, but there are some border cases);
+* automatic differentiation in single-dimension and in multiple dimensions;
+* generation of mathematical syntax trees with subsequent code generation for other languages;
+* symbolic computations, especially differentiation (and some other actions with `kmath-symja` integration with
+  Symja's `IExpr`—integration, simplification, and more);
+* visualization with `kmath-jupyter`.
+
+The workhorse of this API is `Expression` interface, which exposes
+single `operator fun invoke(arguments: Map<Symbol, T>): T`
+method. `ExpressionAlgebra` is used to generate expressions and introduce variables.
+
+Currently there are two implementations:
+
+* Generic `ExpressionField` in `kmath-core` which allows construction of custom lazy expressions
+
+* Auto-differentiation expression in `kmath-commons` module allows using full power of `DerivativeStructure`
+  from commons-math. **TODO: add example**

docs/linear.md

@@ -0,0 +1,36 @@
+## Basic linear algebra layout
+
+KMath support for linear algebra organized in a context-oriented way, which means that operations are in most cases
+declared in context classes, and are not the members of classes that store data. This allows more flexible approach to
+maintain multiple back-ends. The new operations added as extensions to contexts instead of being member functions of
+data structures.
+
+The main context for linear algebra over matrices and vectors is `LinearSpace`, which defines addition and dot products
+of matrices and vectors:
+
+```kotlin
+import space.kscience.kmath.linear.*
+
+LinearSpace.Companion.real {
+    val vec = buildVector(10) { i -> i.toDouble() }
+    val mat = buildMatrix(10, 10) { i, j -> i.toDouble() + j }
+
+    // Addition
+    vec + vec
+    mat + mat
+
+    // Multiplication by scalar
+    vec * 2.0
+    mat * 2.0
+
+    // Dot product
+    mat dot vec
+    mat dot mat
+}
+```
+
+## Backends overview
+
+### EJML
+
+### Commons Math

docs/nd-structure.md

@@ -8,38 +8,45 @@ One of the most sought after features of mathematical libraries is the high-perf
 structures. In `kmath` performance depends on which particular context was used for operation.
 
 Let us consider following contexts:
+
 ```kotlin
     // automatically build context most suited for given type.
-    val autoField = NDField.auto(RealField, dim, dim)
-    // specialized nd-field for Double. It works as generic Double field as well
+    val autoField = NDField.auto(DoubleField, dim, dim)
+    // specialized nd-field for Double. It works as generic Double field as well.
     val specializedField = NDField.real(dim, dim)
     //A generic boxing field. It should be used for objects, not primitives.
-    val genericField = NDField.buffered(RealField, dim, dim)
+    val genericField = NDField.buffered(DoubleField, dim, dim)
 ```
-Now let us perform several tests and see which implementation is best suited for each case:
+
+Now let us perform several tests and see, which implementation is best suited for each case:
 
 ## Test case
 
-In order to test performance we will take 2d-structures with `dim = 1000` and add a structure filled with `1.0`
+To test performance we will take 2d-structures with `dim = 1000` and add a structure filled with `1.0`
 to it `n = 1000` times.
 
 ## Specialized
+
 The code to run this looks like:
+
 ```kotlin
     specializedField.run {
-        var res: NDBuffer<Double> = one
+        var res: NDBuffer<Float64> = one
         repeat(n) {
             res += 1.0
         }
     }
 ```
+
 The performance of this code is the best of all tests since it inlines all operations and is specialized for operation
 with doubles. We will measure everything else relative to this one, so time for this test will be `1x` (real time
-on my computer is about 4.5 seconds). The only problem with this approach is that it requires to specify type
-from the beginning. Everyone do so anyway, so it is the recommended approach.
+on my computer is about 4.5 seconds). The only problem with this approach is that it requires specifying type
+from the beginning. Everyone does so anyway, so it is the recommended approach.
 
 ## Automatic
+
 Let's do the same with automatic field inference:
+
 ```kotlin
     autoField.run {
         var res = one
@@ -48,13 +55,16 @@ Let's do the same with automatic field inference:
         }
     }
 ```
+
 Ths speed of this operation is approximately the same as for specialized case since `NDField.auto` just
-returns the same `RealNDField` in this case. Of course it is usually better to use specialized method to be sure.
+returns the same `RealNDField` in this case. Of course, it is usually better to use specialized method to be sure.
 
 ## Lazy
+
 Lazy field does not produce a structure when asked, instead it generates an empty structure and fills it on-demand
 using coroutines to parallelize computations.
 When one calls
+
 ```kotlin
     lazyField.run {
         var res = one
@@ -63,12 +73,14 @@ When one calls
         }
     }
 ```
-The result will be calculated almost immediately but the result will be empty. In order to get the full result
+
+The result will be calculated almost immediately but the result will be empty. To get the full result
 structure one needs to call all its elements. In this case computation overhead will be huge. So this field never
 should be used if one expects to use the full result structure. Though if one wants only small fraction, it could
 save a lot of time.
 
 This field still could be used with reasonable performance if call code is changed:
+
 ```kotlin
     lazyField.run {
         val res = one.map {
@@ -82,30 +94,37 @@ This field still could be used with reasonable performance if call code is chang
         res.elements().forEach { it.second }
     }
 ```
+
 In this case it completes in about `4x-5x` time due to boxing.
 
 ## Boxing
+
 The boxing field produced by
+
 ```kotlin
     genericField.run {
-        var res: NDBuffer<Double> = one
+        var res: NDBuffer<Float64> = one
         repeat(n) {
             res += 1.0
         }
     }
 ```
-obviously is the slowest one, because it requires to box and unbox the `double` on each operation. It takes about
+
+is the slowest one, because it requires boxing and unboxing the `double` on each operation. It takes about
 `15x` time (**TODO: there seems to be a problem here, it should be slow, but not that slow**). This field should
 never be used for primitives.
 
 ## Element operation
+
 Let us also check the speed for direct operations on elements:
+
 ```kotlin
     var res = genericField.one
     repeat(n) {
         res += 1.0
     }
 ```
+
 One would expect to be at least as slow as field operation, but in fact, this one takes only `2x` time to complete.
 It happens, because in this particular case it does not use actual `NDField` but instead calculated directly
 via extension function.
@@ -114,13 +133,18 @@ via extension function.
 
 Usually it is bad idea to compare the direct numerical operation performance in different languages, but it hard to
 work completely without frame of reference. In this case, simple numpy code:
+
 ```python
+import numpy as np
+
 res = np.ones((1000,1000))
 for i in range(1000):
     res = res + 1.0
 ```
-gives the completion time of about `1.1x`, which means that specialized kotlin code in fact is working faster (I think it is
+
+gives the completion time of about `1.1x`, which means that specialized kotlin code in fact is working faster (I think
+it is
 because better memory management). Of course if one writes `res += 1.0`, the performance will be different,
-but it would be differenc case, because numpy overrides `+=` with in-place operations. In-place operations are
+but it would be different case, because numpy overrides `+=` with in-place operations. In-place operations are
 available in `kmath` with `MutableNDStructure` but there is no field for it (one can still work with mapping
 functions).

docs/polynomials.md

@@ -0,0 +1,223 @@
+# Polynomials and Rational Functions
+
+KMath provides a way to work with uni- and multivariate polynomials and rational functions. It includes full support of
+arithmetic operations of integers, **constants** (elements of ring polynomials are build over), variables (for certain
+multivariate implementations), polynomials and rational functions encapsulated in so-called **polynomial space** and *
+*rational function space** and some other utilities such as algebraic differentiation and substitution.
+
+## Concrete realizations
+
+There are 3 approaches to represent polynomials:
+
+1. For univariate polynomials one can represent and store polynomial as a list of coefficients for each power of the
+   variable. I.e. polynomial $a_0 + \dots + a_n x^n $ can be represented as a finite sequence $(a_0; \dots; a_n)$. (
+   Compare to sequential definition of polynomials.)
+2. For multivariate polynomials one can represent and store polynomial as a matching (in programming it is called "map"
+   or "dictionary", in math it is
+   called [functional relation](https://en.wikipedia.org/wiki/Binary_relation#Special_types_of_binary_relations)) of
+   each "**term signature**" (that describes what variables and in what powers appear in the term) with corresponding
+   coefficient of the term. But there are 2 possible approaches of term signature representation:
+    1. One can number all the variables, so term signature can be represented as a sequence describing powers of the
+       variables. I.e. signature of term $c \\; x_0^{d_0} \dots x_n^{d_n} $ (for natural or zero $d_i $) can be
+       represented as a finite sequence $(d_0; \dots; d_n)$.
+    2. One can represent variables as objects ("**labels**"), so term signature can be also represented as a matching of
+       each appeared variable with its power in the term. I.e. signature of term $c \\; x_0^{d_0} \dots x_n^{d_n} $ (for
+       natural non-zero $d_i $) can be represented as a finite matching $(x_0 \to d_1; \dots; x_n \to d_n)$.
+
+All that three approaches are implemented by "list", "numbered", and "labeled" versions of polynomials and polynomial
+spaces respectively. Whereas all rational functions are represented as fractions with corresponding polynomial numerator
+and denominator, and rational functions' spaces are implemented in the same way as usual field of rational numbers (or
+more precisely, as any field of fractions over integral domain) should be implemented.
+
+So here are a bit of details. Let `C` by type of constants. Then:
+
+1. `ListPolynomial`, `ListPolynomialSpace`, `ListRationalFunction` and `ListRationalFunctionSpace` implement the first
+   scenario. `ListPolynomial` stores polynomial $a_0 + \dots + a_n x^n $ as a coefficients
+   list `listOf(a_0, ..., a_n)` (of type `List<C>`).
+
+   They also have variation `ScalableListPolynomialSpace` that replaces former polynomials and
+   implements `ScaleOperations`.
+2. `NumberedPolynomial`, `NumberedPolynomialSpace`, `NumberedRationalFunction` and `NumberedRationalFunctionSpace`
+   implement second scenario. `NumberedPolynomial` stores polynomials as structures of type `Map<List<UInt>, C>`.
+   Signatures are stored as `List<UInt>`. To prevent ambiguity signatures should not end with zeros.
+3. `LabeledPolynomial`, `LabeledPolynomialSpace`, `LabeledRationalFunction` and `LabeledRationalFunctionSpace` implement
+   third scenario using common `Symbol` as variable type. `LabeledPolynomial` stores polynomials as structures of
+   type `Map<Map<Symbol, UInt>, C>`. Signatures are stored as `Map<Symbol, UInt>`. To prevent ambiguity each signature
+   should not map any variable to zero.
+
+### Example: `ListPolynomial`
+
+For example, polynomial $2 - 3x + x^2 $ (with `Int` coefficients) is represented
+
+```kotlin
+val polynomial: ListPolynomial<Int> = ListPolynomial(listOf(2, -3, 1))
+// or
+val polynomial: ListPolynomial<Int> = ListPolynomial(2, -3, 1)
+```
+
+All algebraic operations can be used in corresponding space:
+
+```kotlin
+val computationResult = Int.algebra.listPolynomialSpace {
+   ListPolynomial(2, -3, 1) + ListPolynomial(0, 6) == ListPolynomial(2, 3, 1)
+}
+
+println(computationResult) // true
+```
+
+For more see [examples](../examples/src/main/kotlin/space/kscience/kmath/functions/polynomials.kt).
+
+### Example: `NumberedPolynomial`
+
+For example, polynomial $3 + 5 x_1 - 7 x_0^2 x_2 $ (with `Int` coefficients) is represented
+
+```kotlin
+val polynomial: NumberedPolynomial<Int> = NumberedPolynomial(
+   mapOf(
+      listOf<UInt>() to 3,
+      listOf(0u, 1u) to 5,
+      listOf(2u, 0u, 1u) to -7,
+   )
+)
+// or
+val polynomial: NumberedPolynomial<Int> = NumberedPolynomial(
+   listOf<UInt>() to 3,
+   listOf(0u, 1u) to 5,
+   listOf(2u, 0u, 1u) to -7,
+)
+```
+
+All algebraic operations can be used in corresponding space:
+
+```kotlin
+val computationResult = Int.algebra.numberedPolynomialSpace {
+   NumberedPolynomial(
+      listOf<UInt>() to 3,
+      listOf(0u, 1u) to 5,
+      listOf(2u, 0u, 1u) to -7,
+   ) + NumberedPolynomial(
+      listOf(0u, 1u) to -5,
+      listOf(0u, 0u, 0u, 4u) to 4,
+   ) == NumberedPolynomial(
+      listOf<UInt>() to 3,
+      listOf(0u, 1u) to 0,
+      listOf(2u, 0u, 1u) to -7,
+      listOf(0u, 0u, 0u, 4u) to 4,
+   )
+}
+
+println(computationResult) // true
+```
+
+For more see [examples](../examples/src/main/kotlin/space/kscience/kmath/functions/polynomials.kt).
+
+### Example: `LabeledPolynomial`
+
+For example, polynomial $3 + 5 y - 7 x^2 z $ (with `Int` coefficients) is represented
+
+```kotlin
+val polynomial: LabeledPolynomial<Int> = LabeledPolynomial(
+   mapOf(
+      mapOf<Symbol, UInt>() to 3,
+      mapOf(y to 1u) to 5,
+      mapOf(x to 2u, z to 1u) to -7,
+   )
+)
+// or
+val polynomial: LabeledPolynomial<Int> = LabeledPolynomial(
+   mapOf<Symbol, UInt>() to 3,
+   mapOf(y to 1u) to 5,
+   mapOf(x to 2u, z to 1u) to -7,
+)
+```
+
+All algebraic operations can be used in corresponding space:
+
+```kotlin
+val computationResult = Int.algebra.labeledPolynomialSpace {
+   LabeledPolynomial(
+      listOf<UInt>() to 3,
+      listOf(0u, 1u) to 5,
+      listOf(2u, 0u, 1u) to -7,
+   ) + LabeledPolynomial(
+      listOf(0u, 1u) to -5,
+      listOf(0u, 0u, 0u, 4u) to 4,
+   ) == LabeledPolynomial(
+      listOf<UInt>() to 3,
+      listOf(0u, 1u) to 0,
+      listOf(2u, 0u, 1u) to -7,
+      listOf(0u, 0u, 0u, 4u) to 4,
+   )
+}
+
+println(computationResult) // true
+```
+
+For more see [examples](../examples/src/main/kotlin/space/kscience/kmath/functions/polynomials.kt).
+
+## Abstract entities (interfaces and abstract classes)
+
+```mermaid
+classDiagram
+   Polynomial <|-- ListPolynomial
+   Polynomial <|-- NumberedPolynomial
+   Polynomial <|-- LabeledPolynomial
+   
+   RationalFunction <|-- ListRationalFunction
+   RationalFunction <|-- NumberedRationalFunction
+   RationalFunction <|-- LabeledRationalFunction
+   
+   Ring <|-- PolynomialSpace
+   PolynomialSpace <|-- MultivariatePolynomialSpace
+   PolynomialSpace <|-- PolynomialSpaceOverRing
+   
+   Ring <|-- RationalFunctionSpace
+   RationalFunctionSpace <|-- MultivariateRationalFunctionSpace
+   RationalFunctionSpace <|-- RationalFunctionSpaceOverRing
+   RationalFunctionSpace <|-- RationalFunctionSpaceOverPolynomialSpace
+   RationalFunctionSpace <|-- PolynomialSpaceOfFractions
+   RationalFunctionSpaceOverPolynomialSpace <|-- MultivariateRationalFunctionSpaceOverMultivariatePolynomialSpace
+   MultivariateRationalFunctionSpace <|-- MultivariateRationalFunctionSpaceOverMultivariatePolynomialSpace
+   MultivariateRationalFunctionSpace <|-- MultivariatePolynomialSpaceOfFractions
+   PolynomialSpaceOfFractions <|-- MultivariatePolynomialSpaceOfFractions
+```
+
+There are implemented `Polynomial` and `RationalFunction` interfaces as abstractions of polynomials and rational
+functions respectively (although, there is not a lot of logic in them) and `PolynomialSpace`
+and `RationalFunctionSpace` (that implement `Ring` interface) as abstractions of polynomials' and rational functions'
+spaces respectively. More precisely, that means they allow to declare common logic of interaction with such objects and
+spaces:
+
+- `Polynomial` does not provide any logic. It is marker interface.
+- `RationalFunction` provides numerator and denominator of rational function and destructuring declaration for them.
+- `PolynomialSpace` provides all possible arithmetic interactions of integers, constants (of type `C`), and
+  polynomials (of type `P`) like addition, subtraction, multiplication, and some others and common properties like
+  degree of polynomial.
+- `RationalFunctionSpace` provides the same as `PolynomialSpace` but also for rational functions: all possible
+  arithmetic interactions of integers, constants (of type `C`), polynomials (of type `P`), and rational functions (of
+  type `R`) like addition, subtraction, multiplication, division (in some cases), and some others and common properties
+  like degree of polynomial.
+
+Then to add abstraction of similar behaviour with variables (in multivariate case) there are
+implemented `MultivariatePolynomialSpace` and `MultivariateRationalFunctionSpace`. They just include variables (of
+type `V`) in the interactions of the entities.
+
+Also, to remove boilerplates there were provided helping subinterfaces and abstract subclasses:
+
+- `PolynomialSpaceOverRing` allows to replace implementation of interactions of integers and constants with
+  implementations from provided ring over constants (of type `A: Ring<C>`).
+- `RationalFunctionSpaceOverRing` &mdash; the same but for `RationalFunctionSpace`.
+- `RationalFunctionSpaceOverPolynomialSpace` &mdash; the same but "the inheritance" includes interactions with
+  polynomials from provided `PolynomialSpace`.
+- `PolynomialSpaceOfFractions` is actually abstract subclass of `RationalFunctionSpace` that implements all fractions
+  boilerplates with provided (`protected`) constructor of rational functions by polynomial numerator and denominator.
+- `MultivariateRationalFunctionSpaceOverMultivariatePolynomialSpace` and `MultivariatePolynomialSpaceOfFractions`
+  &mdash; the same stories of operators inheritance and fractions boilerplates respectively but in multivariate case.
+
+## Utilities
+
+For all kinds of polynomials there are provided (implementation details depend on kind of polynomials) such common
+utilities as:
+
+1. differentiation and anti-differentiation,
+2. substitution, invocation and functional representation.

docs/readme.md

@@ -0,0 +1,14 @@
+# Documentation
+
+* [Algebra](algebra.md): [context-based](contexts.md) operations on different primitives and structures.
+
+* [NDStructures](nd-structure.md)
+
+* [Linear algebra](linear.md): matrices, operations and linear equations solving. To be moved to separate module.
+  Currently, supports basic API and multiple library back-ends.
+
+* [Histograms](histograms.md): multidimensional histogram calculation and operations.
+
+* [Expressions](expressions.md)
+
+* Commons math integration

docs/templates/ARTIFACT-TEMPLATE.md

@@ -0,0 +1,16 @@
+## Artifact:
+
+The Maven coordinates of this project are `${group}:${name}:${version}`.
+
+**Gradle:**
+
+```kotlin
+repositories {
+    maven("https://repo.kotlin.link")
+    mavenCentral()
+}
+
+dependencies {
+    implementation("${group}:${name}:${version}")
+}
+```

docs/templates/README-TEMPLATE.md

@@ -0,0 +1,109 @@
+[![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub)
+[![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
+![Gradle build](https://github.com/SciProgCentre/kmath/workflows/Gradle%20build/badge.svg)
+[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22)
+[![Space](https://img.shields.io/badge/dynamic/xml?color=orange&label=Space&query=//metadata/versioning/latest&url=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fspace%2Fkscience%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/)
+
+# KMath
+
+Could be pronounced as `key-math`. The **K**otlin **Math**ematics 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.
+
+[Documentation site](https://SciProgCentre.github.io/kmath/)
+
+## 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)
+* [A talk at KotlinConf 2019 about using kotlin for science](https://youtu.be/LI_5TZ7tnOE?si=4LknX41gl_YeUbIe)
+* [A talk on architecture at Joker-2021 (in Russian)](https://youtu.be/1bZ2doHiRRM?si=9w953ro9yu98X_KJ)
+* [The same talk in English](https://youtu.be/yP5DIc2fVwQ?si=louZzQ1dcXV6gP10)
+* [A seminar on tensor API](https://youtu.be/0H99wUs0xTM?si=6c__04jrByFQtVpo)
+
+# Goal
+
+* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS, Native and
+  Wasm).
+* 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 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 `kmath-for-real`, which will give better
+  experience for those, who want to work with specific types.
+
+## Features and stability
+
+KMath is a modular library. Different modules provide different features with different API stability guarantees. All
+core modules are released with the same version, but with different API change policy. The features are described in
+module definitions below. The module stability could have the following levels:
+
+* **PROTOTYPE**. On this level there are no compatibility guarantees. All methods and classes form those modules could
+  break any moment. You can still use it, but be sure to fix the specific version.
+* **EXPERIMENTAL**. The general API is decided, but some changes could be made. Volatile API is marked
+  with `@UnstableKMathAPI` or other stability warning annotations.
+* **DEVELOPMENT**. API breaking generally follows semantic versioning ideology. There could be changes in minor
+  versions, but not in patch versions. API is protected
+  with [binary-compatibility-validator](https://github.com/Kotlin/binary-compatibility-validator) tool.
+* **STABLE**. The API stabilized. Breaking changes are allowed only in major releases.
+
+## Modules
+
+${modules}
+
+## Multi-platform support
+
+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, Kotlin/JVM is the primary platform, however, Kotlin/Native and Kotlin/JS contributions and
+feedback are also welcome.
+
+## Performance
+
+Calculation of performance is one of the major goals of KMath in the future, but in some cases it is impossible to
+achieve both
+performance and flexibility.
+
+We expect to focus on creating a 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.
+
+## Requirements
+
+KMath currently relies on JDK 11 for compilation and execution of Kotlin-JVM part. We recommend using GraalVM-CE or
+Oracle GraalVM for execution to get better performance.
+
+### Repositories
+
+Release and development artifacts are accessible from mipt-npm [Space](https://www.jetbrains.com/space/)
+repository `https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven` (see documentation of
+[Kotlin Multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) for more details). The repository could
+be reached through [repo.kotlin.link](https://repo.kotlin.link) proxy:
+
+```kotlin
+repositories {
+    maven("https://repo.kotlin.link")
+}
+
+dependencies {
+    api("${group}:kmath-core:$version")
+    // api("${group}:kmath-core-jvm:$version") for jvm-specific version
+}
+```
+
+## Contributing
+
+The project requires a lot of additional work. The most important thing we need is 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 [good first issue](hhttps://github.com/SciProgCentre/kmath/issues?q=is%3Aissue+is%3Aopen+label%3A%22good+first+issue%22)
+label.

examples/README.md

@@ -0,0 +1,4 @@
+# Module examples
+
+
+

examples/build.gradle.kts

@@ -1,74 +1,75 @@
-import org.jetbrains.kotlin.allopen.gradle.AllOpenExtension
-import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
+import org.jetbrains.kotlin.gradle.tasks.KotlinJvmCompile
 
 plugins {
-    java
     kotlin("jvm")
-    kotlin("plugin.allopen") version "1.3.72"
-    id("kotlinx.benchmark") version "0.2.0-dev-8"
-}
-
-configure<AllOpenExtension> {
-    annotation("org.openjdk.jmh.annotations.State")
 }
 
 repositories {
-    maven("http://dl.bintray.com/kyonifer/maven")
-    maven("https://dl.bintray.com/mipt-npm/scientifik")
-    maven("https://dl.bintray.com/mipt-npm/dev")
     mavenCentral()
-}
-
-sourceSets {
-    register("benchmarks")
+    maven("https://repo.kotlin.link")
+    maven("https://maven.pkg.jetbrains.space/kotlin/p/kotlin/kotlin-js-wrappers")
 }
 
 dependencies {
     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-koma"))
+    implementation(project(":kmath-complex"))
+    implementation(project(":kmath-functions"))
+    implementation(project(":kmath-optimization"))
+    implementation(project(":kmath-stat"))
     implementation(project(":kmath-viktor"))
     implementation(project(":kmath-dimensions"))
-    implementation(project(":kmath-bignum"))
-    implementation(project(":kmath-memory"))
-    implementation("com.kyonifer:koma-core-ejml:0.12")
-    implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:0.2.0-npm-dev-6")
-    implementation("org.slf4j:slf4j-simple:1.7.30")
-    "benchmarksImplementation"("org.jetbrains.kotlinx:kotlinx.benchmark.runtime-jvm:0.2.0-dev-8")
-    "benchmarksImplementation"(sourceSets.main.get().output.plus(sourceSets.main.get().runtimeClasspath))
+    implementation(project(":kmath-ejml"))
+    implementation(project(":kmath-nd4j"))
+    implementation(project(":kmath-tensors"))
+    implementation(project(":kmath-symja"))
+    implementation(project(":kmath-for-real"))
+    //jafama
+    implementation(project(":kmath-jafama"))
+    //multik
+    implementation(project(":kmath-multik"))
+    implementation(libs.multik.default)
+
+    //datetime
+    implementation("org.jetbrains.kotlinx:kotlinx-datetime:0.4.0")
+
+    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.slf4j:slf4j-simple:1.7.32")
+    // plotting
+    implementation("space.kscience:plotlykt-server:0.7.0")
 }
 
-// Configure benchmark
-benchmark {
-    // Setup configurations
-    targets {
-        // This one matches sourceSet name above
-        register("benchmarks")
-    }
-
-    configurations {
-        register("fast") {
-            warmups = 5 // 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
+kotlin {
+    jvmToolchain(11)
+    sourceSets.all {
+        languageSettings {
+            optIn("kotlin.contracts.ExperimentalContracts")
+            optIn("kotlin.ExperimentalUnsignedTypes")
+            optIn("space.kscience.kmath.UnstableKMathAPI")
         }
     }
 }
 
-kotlin.sourceSets.all {
-    with(languageSettings) {
-        useExperimentalAnnotation("kotlin.contracts.ExperimentalContracts")
-        useExperimentalAnnotation("kotlin.ExperimentalUnsignedTypes")
+tasks.withType<KotlinJvmCompile> {
+    compilerOptions {
+        freeCompilerArgs.addAll("-Xjvm-default=all", "-Xopt-in=kotlin.RequiresOptIn", "-Xlambdas=indy")
     }
 }
 
-tasks.withType<KotlinCompile> {
-    kotlinOptions {
-        jvmTarget = Scientifik.JVM_TARGET.toString()
-        freeCompilerArgs = freeCompilerArgs + "-Xopt-in=kotlin.RequiresOptIn"
-    }
+readme {
+    maturity = space.kscience.gradle.Maturity.EXPERIMENTAL
 }

examples/notebooks/Naive classifier.ipynb

@@ -0,0 +1,418 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "source": [
+    "%use kmath(0.3.1-dev-5)\n",
+    "%use plotly(0.5.0)\n",
+    "@file:DependsOn(\"space.kscience:kmath-commons:0.3.1-dev-5\")"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "lQbSB87rNAn9lV6poArVWW",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "//Uncomment to work in Jupyter classic or DataLore\n",
+    "//Plotly.jupyter.notebook()"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "0UP158hfccGgjQtHz0wAi6",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# The model\n",
+    "\n",
+    "Defining the input data format, the statistic abstraction and the statistic implementation based on a weighted sum of elements."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "class XYValues(val xValues: DoubleArray, val yValues: DoubleArray) {\n",
+    "    init {\n",
+    "        require(xValues.size == yValues.size)\n",
+    "    }\n",
+    "}\n",
+    "\n",
+    "fun interface XYStatistic {\n",
+    "    operator fun invoke(values: XYValues): Double\n",
+    "}\n",
+    "\n",
+    "class ConvolutionalXYStatistic(val weights: DoubleArray) : XYStatistic {\n",
+    "    override fun invoke(values: XYValues): Double {\n",
+    "        require(weights.size == values.yValues.size)\n",
+    "        val norm = values.yValues.sum()\n",
+    "        return values.yValues.zip(weights) { value, weight -> value * weight }.sum()/norm\n",
+    "    }\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "Zhgz1Ui91PWz0meJiQpHol",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# Generator\n",
+    "Generate sample data for parabolas and hyperbolas"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "outputs": [],
+   "source": [
+    "fun generateParabolas(xValues: DoubleArray, a: Double, b: Double, c: Double): XYValues {\n",
+    "    val yValues = xValues.map { x -> a * x * x + b * x + c }.toDoubleArray()\n",
+    "    return XYValues(xValues, yValues)\n",
+    "}\n",
+    "\n",
+    "fun generateHyperbols(xValues: DoubleArray, gamma: Double, x0: Double, y0: Double): XYValues {\n",
+    "    val yValues = xValues.map { x -> y0 + gamma / (x - x0) }.toDoubleArray()\n",
+    "    return XYValues(xValues, yValues)\n",
+    "}"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "val xValues = (1.0..10.0).step(1.0).toDoubleArray()\n",
+    "\n",
+    "val xy = generateHyperbols(xValues, 1.0, 0.0, 0.0)\n",
+    "\n",
+    "Plotly.plot {\n",
+    "    scatter {\n",
+    "        this.x.doubles = xValues\n",
+    "        this.y.doubles = xy.yValues\n",
+    "    }\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "ZE2atNvFzQsCvpAF8KK4ch",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Create a default statistic with uniform weights"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "val statistic = ConvolutionalXYStatistic(DoubleArray(xValues.size){1.0})\n",
+    "statistic(xy)"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "EA5HaydTddRKYrtAUwd29h",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "import kotlin.random.Random\n",
+    "\n",
+    "val random = Random(1288)\n",
+    "\n",
+    "val parabolas = buildList{\n",
+    "    repeat(500){\n",
+    "        add(\n",
+    "            generateParabolas(\n",
+    "                xValues, \n",
+    "                random.nextDouble(), \n",
+    "                random.nextDouble(), \n",
+    "                random.nextDouble()\n",
+    "            )\n",
+    "        )\n",
+    "    }\n",
+    "}\n",
+    "\n",
+    "val hyperbolas: List<XYValues> =  buildList{\n",
+    "    repeat(500){\n",
+    "        add(\n",
+    "            generateHyperbols(\n",
+    "                xValues, \n",
+    "                random.nextDouble()*10, \n",
+    "                random.nextDouble(), \n",
+    "                random.nextDouble()\n",
+    "            )\n",
+    "        )\n",
+    "    }\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "t5t6IYmD7Q1ykeo9uijFfQ",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "Plotly.plot { \n",
+    "    scatter { \n",
+    "        x.doubles = xValues\n",
+    "        y.doubles = parabolas[257].yValues\n",
+    "    }\n",
+    "    scatter { \n",
+    "        x.doubles = xValues\n",
+    "        y.doubles = hyperbolas[252].yValues\n",
+    "    }\n",
+    " }"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "oXB8lmju7YVYjMRXITKnhO",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "Plotly.plot { \n",
+    "    histogram { \n",
+    "        name = \"parabolae\"\n",
+    "        x.numbers = parabolas.map { statistic(it) }\n",
+    "    }\n",
+    "    histogram { \n",
+    "        name = \"hyperbolae\"\n",
+    "        x.numbers = hyperbolas.map { statistic(it) }\n",
+    "    }\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "8EIIecUZrt2NNrOkhxG5P0",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "val lossFunction: (XYStatistic) -> Double = { statistic ->\n",
+    "    - abs(parabolas.sumOf { statistic(it) } - hyperbolas.sumOf { statistic(it) })\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "h7UmglJW5zXkAfKHK40oIL",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Using commons-math optimizer to optimize weights"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "import org.apache.commons.math3.optim.*\n",
+    "import org.apache.commons.math3.optim.nonlinear.scalar.*\n",
+    "import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.*\n",
+    "\n",
+    "val optimizer = SimplexOptimizer(1e-1, Double.MAX_VALUE)\n",
+    "\n",
+    "val result = optimizer.optimize(\n",
+    "    ObjectiveFunction { point ->\n",
+    "        lossFunction(ConvolutionalXYStatistic(point))\n",
+    "    },\n",
+    "    NelderMeadSimplex(xValues.size),\n",
+    "    InitialGuess(DoubleArray(xValues.size){ 1.0 }),\n",
+    "    GoalType.MINIMIZE,\n",
+    "    MaxEval(100000)\n",
+    ")"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "0EG3K4aCUciMlgGQKPvJ57",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "Print resulting weights of optimization"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "result.point"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "LelUlY0ZSlJEO9yC6SLk5B",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "Plotly.plot { \n",
+    "    scatter { \n",
+    "        y.doubles = result.point\n",
+    "     }\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "AuFOq5t9KpOIkGrOLsVXNf",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "# The resulting statistic distribution"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "source": [
+    "val resultStatistic = ConvolutionalXYStatistic(result.point)\n",
+    "Plotly.plot { \n",
+    "    histogram { \n",
+    "        name = \"parabolae\"\n",
+    "        x.numbers = parabolas.map { resultStatistic(it) }\n",
+    "    }\n",
+    "    histogram { \n",
+    "        name = \"hyperbolae\"\n",
+    "        x.numbers = hyperbolas.map { resultStatistic(it) }\n",
+    "    }\n",
+    "}"
+   ],
+   "execution_count": null,
+   "outputs": [],
+   "metadata": {
+    "datalore": {
+     "node_id": "zvmq42DRdM5mZ3SpzviHwI",
+     "type": "CODE",
+     "hide_input_from_viewers": false,
+     "hide_output_from_viewers": false
+    }
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "outputs": [],
+   "source": [],
+   "metadata": {
+    "collapsed": false
+   }
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Kotlin",
+   "language": "kotlin",
+   "name": "kotlin"
+  },
+  "datalore": {
+   "version": 1,
+   "computation_mode": "JUPYTER",
+   "package_manager": "pip",
+   "base_environment": "default",
+   "packages": []
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

examples/src/benchmarks/kotlin/scientifik/kmath/bignum/BigDecimalBenchmark.kt

@@ -1,41 +0,0 @@
-package scientifik.kmath.bignum
-
-import org.openjdk.jmh.annotations.Benchmark
-import org.openjdk.jmh.annotations.Scope
-import org.openjdk.jmh.annotations.State
-import scientifik.kmath.operations.JBigDecimalField
-import scientifik.kmath.operations.invoke
-import kotlin.random.Random
-
-@State(Scope.Benchmark)
-class BigDecimalBenchmark {
-    final var times: Int = 10000
-
-    @Benchmark
-    fun java() {
-        val random = Random(0)
-        var sum = JBigDecimalField.zero
-
-        repeat(times) {
-            sum += JBigDecimalField {
-                number(random.nextDouble(-1000000000.0, 1000000000.0)) / random.nextDouble(-1000000.0, 1000000.0)
-            }
-        }
-
-        println("java:$sum")
-    }
-
-    @Benchmark
-    fun bignum() {
-        val random = Random(0)
-        var sum = BigDecimalField.zero
-
-        repeat(times) {
-            sum += BigDecimalField {
-                number(random.nextDouble(-1000000000.0, 1000000000.0)) / random.nextDouble(-1000000.0, 1000000.0)
-            }
-        }
-
-        println("bignum:$sum")
-    }
-}

examples/src/benchmarks/kotlin/scientifik/kmath/bignum/BigIntegerBenchmark.kt

@@ -1,38 +0,0 @@
-package scientifik.kmath.bignum
-
-import org.openjdk.jmh.annotations.Benchmark
-import org.openjdk.jmh.annotations.Scope
-import org.openjdk.jmh.annotations.State
-import scientifik.kmath.operations.BigIntField
-import scientifik.kmath.operations.JBigIntegerField
-import scientifik.kmath.operations.invoke
-import kotlin.random.Random
-
-@State(Scope.Benchmark)
-class BigIntegerBenchmark {
-    final var times: Int = 1000000
-
-    @Benchmark
-    fun java() {
-        val random = Random(0)
-        var sum = JBigIntegerField.zero
-        repeat(times) { sum += JBigIntegerField { number(random.nextInt()) * random.nextInt() } }
-        println("java:$sum")
-    }
-
-    @Benchmark
-    fun bignum() {
-        val random = Random(0)
-        var sum = BigIntegerRing.zero
-        repeat(times) { sum += BigIntegerRing { number(random.nextInt()) * random.nextInt() } }
-        println("bignum:$sum")
-    }
-
-    @Benchmark
-    fun bigint() {
-        val random = Random(0)
-        var sum = BigIntField.zero
-        repeat(times) { sum += BigIntField { number(random.nextInt()) * random.nextInt() } }
-        println("bigint:$sum")
-    }
-}

examples/src/benchmarks/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt

@@ -1,48 +0,0 @@
-package scientifik.kmath.structures
-
-import org.openjdk.jmh.annotations.Benchmark
-import org.openjdk.jmh.annotations.Scope
-import org.openjdk.jmh.annotations.State
-import java.nio.IntBuffer
-
-
-@State(Scope.Benchmark)
-class ArrayBenchmark {
-
-    @Benchmark
-    fun benchmarkArrayRead() {
-        var res = 0
-        for (i in 1..size) {
-            res += array[size - i]
-        }
-    }
-
-    @Benchmark
-    fun benchmarkBufferRead() {
-        var res = 0
-        for (i in 1..size) {
-            res += arrayBuffer.get(size - i)
-        }
-    }
-
-    @Benchmark
-    fun nativeBufferRead() {
-        var res = 0
-        for (i in 1..size) {
-            res += nativeBuffer.get(size - i)
-        }
-    }
-
-    companion object {
-        val size = 1000
-
-        val array = IntArray(size) { it }
-        val arrayBuffer = IntBuffer.wrap(array)
-        val nativeBuffer = IntBuffer.allocate(size).also {
-            for (i in 0 until size) {
-                it.put(i, i)
-            }
-
-        }
-    }
-}

examples/src/benchmarks/kotlin/scientifik/kmath/structures/BufferBenchmark.kt

@@ -1,33 +0,0 @@
-package scientifik.kmath.structures
-
-import org.openjdk.jmh.annotations.Benchmark
-import org.openjdk.jmh.annotations.Scope
-import org.openjdk.jmh.annotations.State
-import scientifik.kmath.operations.Complex
-import scientifik.kmath.operations.complex
-
-@State(Scope.Benchmark)
-class BufferBenchmark {
-
-    @Benchmark
-    fun genericRealBufferReadWrite() {
-        val buffer = RealBuffer(size){it.toDouble()}
-
-        (0 until size).forEach {
-            buffer[it]
-        }
-    }
-
-    @Benchmark
-    fun complexBufferReadWrite() {
-        val buffer = MutableBuffer.complex(size / 2){Complex(it.toDouble(), -it.toDouble())}
-
-        (0 until size / 2).forEach {
-            buffer[it]
-        }
-    }
-
-    companion object {
-        const val size = 100
-    }
-}

examples/src/benchmarks/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt

@@ -1,50 +0,0 @@
-package scientifik.kmath.structures
-
-import org.openjdk.jmh.annotations.Benchmark
-import org.openjdk.jmh.annotations.Scope
-import org.openjdk.jmh.annotations.State
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.invoke
-
-@State(Scope.Benchmark)
-class NDFieldBenchmark {
-    @Benchmark
-    fun autoFieldAdd() {
-        bufferedField {
-            var res: NDBuffer<Double> = one
-            repeat(n) { res += one }
-        }
-    }
-
-    @Benchmark
-    fun autoElementAdd() {
-        var res = genericField.one
-        repeat(n) { res += 1.0 }
-    }
-
-    @Benchmark
-    fun specializedFieldAdd() {
-        specializedField {
-            var res: NDBuffer<Double> = one
-            repeat(n) { res += 1.0 }
-        }
-    }
-
-
-    @Benchmark
-    fun boxingFieldAdd() {
-        genericField {
-            var res: NDBuffer<Double> = one
-            repeat(n) { res += one }
-        }
-    }
-
-    companion object {
-        val dim = 1000
-        val n = 100
-
-        val bufferedField = NDField.auto(RealField, dim, dim)
-        val specializedField = NDField.real(dim, dim)
-        val genericField = NDField.boxing(RealField, dim, dim)
-    }
-}

examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt

@@ -1,61 +0,0 @@
-package scientifik.kmath.structures
-
-import org.jetbrains.bio.viktor.F64Array
-import org.openjdk.jmh.annotations.Benchmark
-import org.openjdk.jmh.annotations.Scope
-import org.openjdk.jmh.annotations.State
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.invoke
-import scientifik.kmath.viktor.ViktorNDField
-
-@State(Scope.Benchmark)
-class ViktorBenchmark {
-    final val dim = 1000
-    final val n = 100
-
-    // automatically build context most suited for given type.
-    final val autoField: BufferedNDField<Double, RealField> = NDField.auto(RealField, dim, dim)
-    final val realField: RealNDField = NDField.real(dim, dim)
-    final val viktorField: ViktorNDField = ViktorNDField(intArrayOf(dim, dim))
-
-    @Benchmark
-    fun automaticFieldAddition() {
-        autoField {
-            var res = one
-            repeat(n) { res += one }
-        }
-    }
-
-    @Benchmark
-    fun viktorFieldAddition() {
-        viktorField {
-            var res = one
-            repeat(n) { res += one }
-        }
-    }
-
-    @Benchmark
-    fun rawViktor() {
-        val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim))
-        var res = one
-        repeat(n) { res = res + one }
-    }
-
-    @Benchmark
-    fun realdFieldLog() {
-        realField {
-            val fortyTwo = produce { 42.0 }
-            var res = one
-            repeat(n) { res = ln(fortyTwo) }
-        }
-    }
-
-    @Benchmark
-    fun rawViktorLog() {
-        val fortyTwo = F64Array.full(dim, dim, init = 42.0)
-        var res: F64Array
-        repeat(n) {
-            res = fortyTwo.log()
-        }
-    }
-}

examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt

@@ -1,11 +0,0 @@
-package scientifik.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")
-}

examples/src/main/kotlin/scientifik/kmath/ast/ExpressionsInterpretersBenchmark.kt

@@ -1,71 +0,0 @@
-package scientifik.kmath.ast
-
-import scientifik.kmath.asm.compile
-import scientifik.kmath.expressions.Expression
-import scientifik.kmath.expressions.expressionInField
-import scientifik.kmath.expressions.invoke
-import scientifik.kmath.operations.Field
-import scientifik.kmath.operations.RealField
-import kotlin.random.Random
-import kotlin.system.measureTimeMillis
-
-private 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")
-}

examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionBenchmark.kt

@@ -1,71 +0,0 @@
-package scientifik.kmath.commons.prob
-
-import kotlinx.coroutines.Dispatchers
-import kotlinx.coroutines.async
-import kotlinx.coroutines.runBlocking
-import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler
-import org.apache.commons.rng.simple.RandomSource
-import scientifik.kmath.chains.BlockingRealChain
-import scientifik.kmath.prob.*
-import java.time.Duration
-import java.time.Instant
-
-
-private suspend fun runChain(): Duration {
-    val generator = RandomGenerator.fromSource(RandomSource.MT, 123L)
-
-    val normal = Distribution.normal(NormalSamplerMethod.Ziggurat)
-    val chain = normal.sample(generator) as BlockingRealChain
-
-    val startTime = Instant.now()
-    var sum = 0.0
-    repeat(10000001) { counter ->
-
-        sum += chain.nextDouble()
-
-        if (counter % 100000 == 0) {
-            val duration = Duration.between(startTime, Instant.now())
-            val meanValue = sum / counter
-            println("Chain sampler completed $counter elements in $duration: $meanValue")
-        }
-    }
-    return Duration.between(startTime, Instant.now())
-}
-
-private fun runDirect(): Duration {
-    val provider = RandomSource.create(RandomSource.MT, 123L)
-    val sampler = ZigguratNormalizedGaussianSampler(provider)
-    val startTime = Instant.now()
-
-    var sum = 0.0
-    repeat(10000001) { counter ->
-
-        sum += sampler.sample()
-
-        if (counter % 100000 == 0) {
-            val duration = Duration.between(startTime, Instant.now())
-            val meanValue = sum / counter
-            println("Direct sampler completed $counter elements in $duration: $meanValue")
-        }
-    }
-    return Duration.between(startTime, Instant.now())
-}
-
-/**
- * Comparing chain sampling performance with direct sampling performance
- */
-fun main() {
-    runBlocking(Dispatchers.Default) {
-        val chainJob = async {
-            runChain()
-        }
-
-        val directJob = async {
-            runDirect()
-        }
-
-        println("Chain: ${chainJob.await()}")
-        println("Direct: ${directJob.await()}")
-    }
-
-}

examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt

@@ -1,32 +0,0 @@
-package scientifik.kmath.commons.prob
-
-import kotlinx.coroutines.runBlocking
-import scientifik.kmath.chains.Chain
-import scientifik.kmath.chains.collectWithState
-import scientifik.kmath.prob.Distribution
-import scientifik.kmath.prob.RandomGenerator
-import scientifik.kmath.prob.normal
-
-data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
-
-fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChainState(), { it.copy() }) { chain ->
-    val next = chain.next()
-    num++
-    value += next
-    return@collectWithState value / num
-}
-
-
-fun main() {
-    val normal = Distribution.normal()
-    val chain = normal.sample(RandomGenerator.default).mean()
-
-    runBlocking {
-        repeat(10001) { counter ->
-            val mean = chain.next()
-            if (counter % 1000 == 0) {
-                println("[$counter] Average value is $mean")
-            }
-        }
-    }
-}

examples/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt

@@ -1,55 +0,0 @@
-package scientifik.kmath.linear
-
-import koma.matrix.ejml.EJMLMatrixFactory
-import scientifik.kmath.commons.linear.CMMatrixContext
-import scientifik.kmath.commons.linear.inverse
-import scientifik.kmath.commons.linear.toCM
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.invoke
-import scientifik.kmath.structures.Matrix
-import kotlin.contracts.ExperimentalContracts
-import kotlin.random.Random
-import kotlin.system.measureTimeMillis
-
-@ExperimentalContracts
-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) { val res = inverse(matrix) }
-        val inverseTime = measureTimeMillis { repeat(n) { val res = 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) { val res = inverse(cm) }
-        }
-    }
-
-
-    println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis")
-
-    //koma-ejml
-
-    val komaTime = measureTimeMillis {
-        (KomaMatrixContext(EJMLMatrixFactory(), RealField)) {
-            val km = matrix.toKoma()      //avoid overhead on conversion
-            repeat(n) {
-                val res = inverse(km)
-            }
-        }
-    }
-
-    println("[koma-ejml] Inversion of $n matrices $dim x $dim finished in $komaTime millis")
-}

examples/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt

@@ -1,49 +0,0 @@
-package scientifik.kmath.linear
-
-import koma.matrix.ejml.EJMLMatrixFactory
-import scientifik.kmath.commons.linear.CMMatrixContext
-import scientifik.kmath.commons.linear.toCM
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.invoke
-import scientifik.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")
-    }
-
-    (KomaMatrixContext(EJMLMatrixFactory(), RealField)) {
-        val komaMatrix1 = matrix1.toKoma()
-        val komaMatrix2 = matrix2.toKoma()
-
-        val komaTime = measureTimeMillis {
-            komaMatrix1 dot komaMatrix2
-        }
-
-        println("Koma-ejml implementation time: $komaTime")
-    }
-
-    val genericTime = measureTimeMillis {
-        val res = matrix1 dot matrix2
-    }
-
-    println("Generic implementation time: $genericTime")
-}

examples/src/main/kotlin/scientifik/kmath/operations/ComplexDemo.kt

@@ -1,19 +0,0 @@
-package scientifik.kmath.operations
-
-import scientifik.kmath.structures.NDElement
-import scientifik.kmath.structures.NDField
-import scientifik.kmath.structures.complex
-
-fun main() {
-    val element = NDElement.complex(2, 2) { index: IntArray ->
-        Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble())
-    }
-
-    val compute = (NDField.complex(8)) {
-        val a = produce { (it) -> i * it - it.toDouble() }
-        val b = 3
-        val c = Complex(1.0, 1.0)
-
-        (a pow b) + c
-    }
-}

examples/src/main/kotlin/scientifik/kmath/structures/ComplexND.kt

@@ -1,53 +0,0 @@
-package scientifik.kmath.structures
-
-import scientifik.kmath.linear.transpose
-import scientifik.kmath.operations.Complex
-import scientifik.kmath.operations.ComplexField
-import scientifik.kmath.operations.invoke
-import kotlin.system.measureTimeMillis
-
-fun main() {
-    val dim = 1000
-    val n = 1000
-
-    val realField = NDField.real(dim, dim)
-    val complexField = NDField.complex(dim, dim)
-
-    val realTime = measureTimeMillis {
-        realField {
-            var res: NDBuffer<Double> = one
-            repeat(n) {
-                res += 1.0
-            }
-        }
-    }
-
-    println("Real addition completed in $realTime millis")
-
-    val complexTime = measureTimeMillis {
-        complexField {
-            var res: NDBuffer<Complex> = one
-            repeat(n) { res += 1.0 }
-        }
-    }
-
-    println("Complex addition completed in $complexTime millis")
-}
-
-fun complexExample() {
-    //Create a context for 2-d structure with complex values
-    ComplexField {
-        nd(4, 8) {
-            //a constant real-valued structure
-            val x = one * 2.5
-            operator fun Number.plus(other: Complex) = Complex(this.toDouble() + other.re, other.im)
-            //a structure generator specific to this context
-            val matrix = produce { (k, l) -> k + l * i }
-            //Perform sum
-            val sum = matrix + x + 1.0
-
-            //Represent the sum as 2d-structure and transpose
-            sum.as2D().transpose()
-        }
-    }
-}

examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt

@@ -1,68 +0,0 @@
-package scientifik.kmath.structures
-
-import kotlinx.coroutines.GlobalScope
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.invoke
-import kotlin.contracts.ExperimentalContracts
-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")
-}
-
-fun main() {
-    val dim = 1000
-    val n = 1000
-
-    // automatically build context most suited for given type.
-    val autoField = NDField.auto(RealField, dim, dim)
-    // specialized nd-field for Double. It works as generic Double field as well
-    val specializedField = NDField.real(dim, dim)
-    //A generic boxing field. It should be used for objects, not primitives.
-    val genericField = NDField.boxing(RealField, dim, dim)
-
-    measureAndPrint("Automatic field addition") {
-        autoField {
-            var res: NDBuffer<Double> = one
-            repeat(n) { res += number(1.0) }
-        }
-    }
-
-    measureAndPrint("Element addition") {
-        var res = genericField.one
-        repeat(n) { res += 1.0 }
-    }
-
-    measureAndPrint("Specialized addition") {
-        specializedField {
-            var res: NDBuffer<Double> = one
-            repeat(n) { res += 1.0 }
-        }
-    }
-
-    measureAndPrint("Lazy addition") {
-        val res = specializedField.one.mapAsync(GlobalScope) {
-            var c = 0.0
-            repeat(n) {
-                c += 1.0
-            }
-            c
-        }
-
-        res.elements().forEach { it.second }
-    }
-
-    measureAndPrint("Generic addition") {
-        //genericField.run(action)
-        genericField {
-            var res: NDBuffer<Double> = one
-            repeat(n) {
-                res += one // couldn't avoid using `one` due to resolution ambiguity }
-            }
-        }
-    }
-}

examples/src/main/kotlin/scientifik/kmath/structures/StructureReadBenchmark.kt

@@ -1,36 +0,0 @@
-package scientifik.kmath.structures
-
-import kotlin.system.measureTimeMillis
-
-fun main(args: Array<String>) {
-    val n = 6000
-
-    val array = DoubleArray(n * n) { 1.0 }
-    val buffer = RealBuffer(array)
-    val strides = DefaultStrides(intArrayOf(n, n))
-
-    val structure = BufferNDStructure(strides, buffer)
-
-    measureTimeMillis {
-        var res: Double = 0.0
-        strides.indices().forEach { res = structure[it] }
-    } // warmup
-
-    val time1 = measureTimeMillis {
-        var res: Double = 0.0
-        strides.indices().forEach { res = structure[it] }
-    }
-    println("Structure reading finished in $time1 millis")
-
-    val time2 = measureTimeMillis {
-        var res: Double = 0.0
-        strides.indices().forEach { res = buffer[strides.offset(it)] }
-    }
-    println("Buffer reading finished in $time2 millis")
-
-    val time3 = measureTimeMillis {
-        var res: Double = 0.0
-        strides.indices().forEach { res = array[strides.offset(it)] }
-    }
-    println("Array reading finished in $time3 millis")
-}

examples/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt

@@ -1,38 +0,0 @@
-package scientifik.kmath.structures
-
-import kotlin.system.measureTimeMillis
-
-
-fun main(args: Array<String>) {
-
-    val n = 6000
-
-    val structure = NDStructure.build(intArrayOf(n, n), Buffer.Companion::auto) { 1.0 }
-
-    structure.mapToBuffer { it + 1 } // warm-up
-
-    val time1 = measureTimeMillis {
-        val res = structure.mapToBuffer { it + 1 }
-    }
-    println("Structure mapping finished in $time1 millis")
-
-    val array = DoubleArray(n * n) { 1.0 }
-
-    val time2 = measureTimeMillis {
-        val target = DoubleArray(n * n)
-        val res = array.forEachIndexed { index, value ->
-            target[index] = value + 1
-        }
-    }
-    println("Array mapping finished in $time2 millis")
-
-    val buffer = RealBuffer(DoubleArray(n * n) { 1.0 })
-
-    val time3 = measureTimeMillis {
-        val target = RealBuffer(DoubleArray(n * n))
-        val res = array.forEachIndexed { index, value ->
-            target[index] = value + 1
-        }
-    }
-    println("Buffer mapping finished in $time3 millis")
-}

examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt

@@ -1,32 +0,0 @@
-package scientifik.kmath.structures
-
-import scientifik.kmath.dimensions.D2
-import scientifik.kmath.dimensions.D3
-import scientifik.kmath.dimensions.DMatrixContext
-import scientifik.kmath.dimensions.Dimension
-import scientifik.kmath.operations.RealField
-
-fun DMatrixContext<Double, RealField>.simple() {
-    val m1 = produce<D2, D3> { i, j -> (i + j).toDouble() }
-    val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
-
-    //Dimension-safe addition
-    m1.transpose() + m2
-}
-
-
-object D5 : Dimension {
-    override val dim: UInt = 5u
-}
-
-fun DMatrixContext<Double, RealField>.custom() {
-    val m1 = produce<D2, D5> { i, j -> (i + j).toDouble() }
-    val m2 = produce<D5, D2> { i, j -> (i - j).toDouble() }
-    val m3 = produce<D2, D2> { i, j -> (i - j).toDouble() }
-    (m1 dot m2) + m3
-}
-
-fun main(): Unit = with(DMatrixContext.real) {
-    simple()
-    custom()
-}

examples/src/main/kotlin/space/kscience/kmath/ast/astRendering.kt

@@ -0,0 +1,26 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ast
+
+import space.kscience.kmath.ast.rendering.FeaturedMathRendererWithPostProcess
+import space.kscience.kmath.ast.rendering.LatexSyntaxRenderer
+import space.kscience.kmath.ast.rendering.MathMLSyntaxRenderer
+import space.kscience.kmath.ast.rendering.renderWithStringBuilder
+
+fun main() {
+    val mst = "exp(sqrt(x))-asin(2*x)/(2e10+x^3)/(-12)".parseMath()
+    val syntax = FeaturedMathRendererWithPostProcess.Default.render(mst)
+    println("MathSyntax:")
+    println(syntax)
+    println()
+    val latex = LatexSyntaxRenderer.renderWithStringBuilder(syntax)
+    println("LaTeX:")
+    println(latex)
+    println()
+    val mathML = MathMLSyntaxRenderer.renderWithStringBuilder(syntax)
+    println("MathML:")
+    println(mathML)
+}

examples/src/main/kotlin/space/kscience/kmath/ast/expressions.kt

@@ -0,0 +1,26 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ast
+
+import space.kscience.kmath.asm.compileToExpression
+import space.kscience.kmath.expressions.MstExtendedField
+import space.kscience.kmath.expressions.Symbol.Companion.x
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.invoke
+
+fun main() {
+    val expr = MstExtendedField {
+        x * 2.0 + number(2.0) / x - number(16.0) + asinh(x) / sin(x)
+    }.compileToExpression(Float64Field)
+
+    val m = DoubleArray(expr.indexer.symbols.size)
+    val xIdx = expr.indexer.indexOf(x)
+
+    repeat(10000000) {
+        m[xIdx] = 1.0
+        expr(m)
+    }
+}

examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt

@@ -0,0 +1,27 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ast
+
+import space.kscience.kmath.expressions.Symbol.Companion.x
+import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.expressions.invoke
+import space.kscience.kmath.expressions.toExpression
+import space.kscience.kmath.kotlingrad.toKotlingradExpression
+import space.kscience.kmath.operations.Float64Field
+
+/**
+ * In this example, *x<sup>2</sup> &minus; 4 x &minus; 44* function is differentiated with Kotlin∇, and the
+ * derivation result is compared with valid derivative in a certain point.
+ */
+fun main() {
+    val actualDerivative = "x^2-4*x-44"
+        .parseMath()
+        .toKotlingradExpression(Float64Field)
+        .derivative(x)
+
+    val expectedDerivative = "2*x-4".parseMath().toExpression(Float64Field)
+    check(actualDerivative(x to 123.0) == expectedDerivative(x to 123.0))
+}

examples/src/main/kotlin/space/kscience/kmath/ast/symjaSupport.kt

@@ -0,0 +1,27 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.ast
+
+import space.kscience.kmath.expressions.Symbol.Companion.x
+import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.expressions.invoke
+import space.kscience.kmath.expressions.toExpression
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.symja.toSymjaExpression
+
+/**
+ * In this example, *x<sup>2</sup> &minus; 4 x &minus; 44* function is differentiated with Symja, and the
+ * derivation result is compared with valid derivative in a certain point.
+ */
+fun main() {
+    val actualDerivative = "x^2-4*x-44"
+        .parseMath()
+        .toSymjaExpression(Float64Field)
+        .derivative(x)
+
+    val expectedDerivative = "2*x-4".parseMath().toExpression(Float64Field)
+    check(actualDerivative(x to 123.0) == expectedDerivative(x to 123.0))
+}

examples/src/main/kotlin/space/kscience/kmath/expressions/autodiff.kt

@@ -0,0 +1,92 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.expressions
+
+import space.kscience.kmath.UnstableKMathAPI
+
+// Only kmath-core is needed.
+
+// Let's declare some variables
+val x by symbol
+val y by symbol
+val z by symbol
+
+@OptIn(UnstableKMathAPI::class)
+fun main() {
+    // Let's define some random expression.
+    val someExpression = Double.autodiff.differentiate {
+        // We bind variables `x` and `y` to the builder scope,
+        val x = bindSymbol(x)
+        val y = bindSymbol(y)
+
+        // Then we use the bindings to define expression `xy + x + y - 1`
+        x * y + x + y - 1
+    }
+
+    // Then we can evaluate it at any point ((-1, -1) in the case):
+    println(someExpression(x to -1.0, y to -1.0))
+    // >>> -2.0
+
+    // We can also construct its partial derivatives:
+    val dxExpression = someExpression.derivative(x) // ∂/∂x. Must be `y+1`
+    val dyExpression = someExpression.derivative(y) // ∂/∂y. Must be `x+1`
+    val dxdxExpression = someExpression.derivative(x, x) // ∂^2/∂x^2. Must be `0`
+
+    // We can evaluate them as well
+    println(dxExpression(x to 57.0, y to 6.0))
+    // >>> 7.0
+    println(dyExpression(x to -1.0, y to 179.0))
+    // >>> 0.0
+    println(dxdxExpression(x to 239.0, y to 30.0))
+    // >>> 0.0
+
+    // You can also provide extra arguments that obviously won't affect the result:
+    println(dxExpression(x to 57.0, y to 6.0, z to 42.0))
+    // >>> 7.0
+    println(dyExpression(x to -1.0, y to 179.0, z to 0.0))
+    // >>> 0.0
+    println(dxdxExpression(x to 239.0, y to 30.0, z to 100_000.0))
+    // >>> 0.0
+
+    // But in case you forgot to specify bound symbol's value, exception is thrown:
+    println(runCatching { someExpression(z to 4.0) })
+    // >>> Failure(java.lang.IllegalStateException: Symbol 'x' is not supported in ...)
+
+    // The reason is that the expression is evaluated lazily,
+    // and each `bindSymbol` operation actually substitutes the provided symbol with the corresponding value.
+
+    // For example, let there be an expression
+    val simpleExpression = Double.autodiff.differentiate {
+        val x = bindSymbol(x)
+        x pow 2
+    }
+    // When you evaluate it via
+    simpleExpression(x to 1.0, y to 57.0, z to 179.0)
+    // lambda above has the context of map `{x: 1.0, y: 57.0, z: 179.0}`.
+    // When x is bound, you can think of it as substitution `x -> 1.0`.
+    // Other values are unused which does not make any problem to us.
+    // But in the case the corresponding value is not provided,
+    // we cannot bind the variable. Thus, exception is thrown.
+
+    // There is also a function `bindSymbolOrNull` that fixes the problem:
+    val fixedExpression = Double.autodiff.differentiate {
+        val x = bindSymbolOrNull(x) ?: const(8.0)
+        x pow -2
+    }
+    println(fixedExpression())
+    // >>> 0.015625
+    // It works!
+
+    // The expression provides a bunch of operations:
+    // 1. Constant bindings (via `const` and `number`).
+    // 2. Variable bindings (via `bindVariable`, `bindVariableOrNull`).
+    // 3. Arithmetic operations (via `+`, `-`, `*`, and `-`).
+    // 4. Exponentiation (via `pow` or `power`).
+    // 5. `exp` and `ln`.
+    // 6. Trigonometrical functions (`sin`, `cos`, `tan`, `cot`).
+    // 7. Inverse trigonometrical functions (`asin`, `acos`, `atan`, `acot`).
+    // 8. Hyperbolic functions and inverse hyperbolic functions.
+}

examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt

@@ -0,0 +1,113 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.fit
+
+import kotlinx.html.br
+import kotlinx.html.h3
+import space.kscience.kmath.commons.optimization.CMOptimizer
+import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.expressions.autodiff
+import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.operations.asIterable
+import space.kscience.kmath.operations.toList
+import space.kscience.kmath.optimization.*
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.real.DoubleVector
+import space.kscience.kmath.real.map
+import space.kscience.kmath.real.step
+import space.kscience.kmath.stat.chiSquaredExpression
+import space.kscience.plotly.*
+import space.kscience.plotly.models.ScatterMode
+import space.kscience.plotly.models.TraceValues
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+// Forward declaration of symbols that will be used in expressions.
+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: DoubleVector) {
+    numbers = vector.asIterable()
+}
+
+/**
+ * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules.
+ */
+suspend fun main() {
+    //A generator for a normally distributed values
+    val generator = NormalDistribution(0.0, 1.0)
+
+    //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 { it ->
+        val value = it.pow(2) + it + 1
+        value + chain.next() * 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 = Double.autodiff.chiSquaredExpression(x, y, yErr) { arg ->
+        //bind variables to autodiff context
+        val a = bindSymbol(a)
+        val b = bindSymbol(b)
+        //Include default value for c if it is not provided as a parameter
+        val c = bindSymbolOrNull(c) ?: one
+        a * arg.pow(2) + b * arg + c
+    }
+
+    //minimize the chi^2 in given starting point. Derivatives are not required, they are already included.
+    val result = chi2.optimizeWith(
+        CMOptimizer,
+        mapOf(a to 1.5, b to 0.9, c to 1.0),
+    ) {
+        FunctionOptimizationTarget(OptimizationDirection.MINIMIZE)
+    }
+
+    //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.result[a]!! * it.pow(2) + result.result[b]!! * it + 1 })
+                name = "fit"
+            }
+        }
+        br()
+        h3 {
+            +"Fit result: $result"
+        }
+        h3 {
+            +"Chi2/dof = ${result.resultValue / (x.size - 3)}"
+        }
+    }
+
+    page.makeFile()
+}

examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt

@@ -0,0 +1,112 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.fit
+
+import kotlinx.html.br
+import kotlinx.html.h3
+import space.kscience.attributes.Attributes
+import space.kscience.kmath.data.XYErrorColumnarData
+import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.autodiff
+import space.kscience.kmath.expressions.binding
+import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.operations.asIterable
+import space.kscience.kmath.operations.toList
+import space.kscience.kmath.optimization.*
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.real.map
+import space.kscience.kmath.real.step
+import space.kscience.plotly.*
+import space.kscience.plotly.models.ScatterMode
+import kotlin.math.abs
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+// Forward declaration of symbols that will be used in expressions.
+private val a by symbol
+private val b by symbol
+private val c by symbol
+private val d by symbol
+private val e by symbol
+
+
+/**
+ * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules.
+ */
+suspend fun main() {
+    //A generator for a normally distributed values
+    val generator = NormalDistribution(0.0, 1.0)
+
+    //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 { it ->
+        val value = it.pow(2) + it + 1
+        value + chain.next() * 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(abs(it)) }
+    require(yErr.asIterable().all { it > 0 }) { "All errors must be strictly positive" }
+
+    val result = XYErrorColumnarData.of(x, y, yErr).fitWith(
+        QowOptimizer,
+        Double.autodiff,
+        mapOf(a to 0.9, b to 1.2, c to 2.0, e to 1.0, d to 1.0, e to 0.0),
+        attributes = Attributes(OptimizationParameters, listOf(a, b, c, d))
+    ) { arg ->
+        //bind variables to autodiff context
+        val a by binding
+        val b by binding
+        //Include default value for c if it is not provided as a parameter
+        val c = bindSymbolOrNull(c) ?: one
+        val d by binding
+        val e by binding
+
+        a * arg.pow(2) + b * arg + c + d * arg.pow(3) + e / arg
+    }
+
+    println("Resulting chi2/dof: ${result.chiSquaredOrNull}/${result.dof}")
+
+    //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.model(result.startPoint + result.result + (Symbol.x to it)) })
+                name = "fit"
+            }
+        }
+        br()
+        h3 {
+            +"Fit result: ${result.result}"
+        }
+        h3 {
+            +"Chi2/dof = ${result.chiSquaredOrNull!! / result.dof}"
+        }
+    }
+
+    page.makeFile()
+}

examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt

@@ -0,0 +1,37 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.functions
+
+import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.complex.ComplexField
+import space.kscience.kmath.complex.ComplexField.div
+import space.kscience.kmath.complex.ComplexField.minus
+import space.kscience.kmath.complex.algebra
+import space.kscience.kmath.integration.gaussIntegrator
+import space.kscience.kmath.integration.integrate
+import space.kscience.kmath.integration.value
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
+import kotlin.math.pow
+
+fun main() {
+    //Define a function
+    val function: Function1D<Float64> = { x -> 3 * x.pow(2) + 2 * x + 1 }
+
+    //get the result of the integration
+    val result = Float64Field.gaussIntegrator.integrate(0.0..10.0, function = function)
+
+    //the value is nullable because in some cases the integration could not succeed
+    println(result.value)
+
+
+    repeat(100000) {
+        Complex.algebra.gaussIntegrator.integrate(0.0..1.0, intervals = 1000) { x: Double ->
+//            sin(1 / x) + i * cos(1 / x)
+            1 / x - ComplexField.i / x
+        }.value
+    }
+}

examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt

@@ -0,0 +1,52 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.functions
+
+import space.kscience.kmath.interpolation.SplineInterpolator
+import space.kscience.kmath.interpolation.interpolatePolynomials
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
+import space.kscience.plotly.Plotly
+import space.kscience.plotly.UnstablePlotlyAPI
+import space.kscience.plotly.makeFile
+import space.kscience.plotly.models.functionXY
+import space.kscience.plotly.scatter
+import kotlin.math.PI
+import kotlin.math.sin
+
+@OptIn(UnstablePlotlyAPI::class)
+fun main() {
+    val data = (0..10).map {
+        val x = it.toDouble() / 5 * PI
+        x to sin(x)
+    }
+
+    val polynomial: PiecewisePolynomial<Float64> = SplineInterpolator(Float64Field).interpolatePolynomials(data)
+
+    val function = polynomial.asFunction(Float64Field, 0.0)
+
+    val cmInterpolate = org.apache.commons.math3.analysis.interpolation.SplineInterpolator().interpolate(
+        data.map { it.first }.toDoubleArray(),
+        data.map { it.second }.toDoubleArray()
+    )
+
+    Plotly.plot {
+        scatter {
+            name = "interpolated"
+            x.numbers = data.map { it.first }
+            y.numbers = x.doubles.map { function(it) }
+        }
+        scatter {
+            name = "original"
+            functionXY(0.0..(2 * PI), 0.1) { sin(it) }
+        }
+        scatter {
+            name = "cm"
+            x.numbers = data.map { it.first }
+            y.numbers = x.doubles.map { cmInterpolate.value(it) }
+        }
+    }.makeFile()
+}

examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt

@@ -0,0 +1,46 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.functions
+
+import space.kscience.kmath.interpolation.interpolatePolynomials
+import space.kscience.kmath.interpolation.splineInterpolator
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.real.map
+import space.kscience.kmath.real.step
+import space.kscience.kmath.structures.Float64
+import space.kscience.plotly.Plotly
+import space.kscience.plotly.UnstablePlotlyAPI
+import space.kscience.plotly.makeFile
+import space.kscience.plotly.models.functionXY
+import space.kscience.plotly.scatter
+
+@OptIn(UnstablePlotlyAPI::class)
+fun main() {
+    val function: Function1D<Float64> = { x ->
+        if (x in 30.0..50.0) {
+            1.0
+        } else {
+            0.0
+        }
+    }
+    val xs = 0.0..100.0 step 0.5
+    val ys = xs.map(function)
+
+    val polynomial: PiecewisePolynomial<Float64> = Float64Field.splineInterpolator.interpolatePolynomials(xs, ys)
+
+    val polyFunction = polynomial.asFunction(Float64Field, 0.0)
+
+    Plotly.plot {
+        scatter {
+            name = "interpolated"
+            functionXY(25.0..55.0, 0.1) { polyFunction(it) }
+        }
+        scatter {
+            name = "original"
+            functionXY(25.0..55.0, 0.1) { function(it) }
+        }
+    }.makeFile()
+}

examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt

@@ -0,0 +1,33 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.functions
+
+import space.kscience.kmath.integration.gaussIntegrator
+import space.kscience.kmath.integration.integrate
+import space.kscience.kmath.integration.value
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.structureND
+import space.kscience.kmath.nd.withNdAlgebra
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
+import kotlin.math.pow
+
+fun main(): Unit = Double.algebra.withNdAlgebra(2, 2) {
+
+    //Produce a diagonal StructureND
+    fun diagonal(v: Double) = structureND { (i, j) ->
+        if (i == j) v else 0.0
+    }
+
+    //Define a function in a nd space
+    val function: (Double) -> StructureND<Float64> = { x: Double -> 3 * x.pow(2) + 2 * diagonal(x) + 1 }
+
+    //get the result of the integration
+    val result = gaussIntegrator.integrate(0.0..10.0, function = function)
+
+    //the value is nullable because in some cases the integration could not succeed
+    println(result.value)
+}

examples/src/main/kotlin/space/kscience/kmath/jafama/JafamaDemo.kt

@@ -0,0 +1,15 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.jafama
+
+import space.kscience.kmath.operations.invoke
+
+fun main() {
+    val a = 2.0
+    val b = StrictJafamaDoubleField { exp(a) }
+    println(JafamaDoubleField { b + a })
+    println(StrictJafamaDoubleField { ln(b) })
+}

examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt

@@ -0,0 +1,31 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import kotlin.random.Random
+import kotlin.time.measureTime
+
+fun main() = with(Float64ParallelLinearSpace) {
+    val random = Random(12224)
+    val dim = 1000
+
+    //creating invertible matrix
+    val matrix1 = buildMatrix(dim, dim) { i, j ->
+        if (i <= j) random.nextDouble() else 0.0
+    }
+    val matrix2 = buildMatrix(dim, dim) { i, j ->
+        if (i <= j) random.nextDouble() else 0.0
+    }
+
+    val time = measureTime {
+        repeat(30) {
+            matrix1 dot matrix2
+        }
+    }
+
+    println(time)
+
+}

examples/src/main/kotlin/space/kscience/kmath/linear/eigenValueDecomposition.kt

@@ -0,0 +1,40 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.kmath.commons.linear.CMLinearSpace
+import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
+import kotlin.random.Random
+
+fun main() {
+    val dim = 46
+
+    val random = Random(123)
+
+    val u = Float64.algebra.linearSpace.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+
+    listOf(CMLinearSpace, EjmlLinearSpaceDDRM).forEach { algebra ->
+        with(algebra) {
+            //create a simmetric matrix
+            val matrix = buildMatrix(dim, dim) { row, col ->
+                if (row >= col) u[row, col] else u[col, row]
+            }
+            val eigen = matrix.getOrComputeAttribute(EIG) ?: error("Failed to compute eigenvalue decomposition")
+            check(
+                StructureND.contentEquals(
+                    matrix,
+                    eigen.v dot eigen.d dot eigen.v.transposed(),
+                    1e-4
+                )
+            ) { "$algebra decomposition failed" }
+            println("$algebra eigenvalue decomposition complete and checked" )
+        }
+    }
+
+}

examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt

@@ -0,0 +1,34 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.kmath.real.*
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.structures.Float64Buffer
+
+fun main() {
+    val x0 = DoubleVector(0.0, 0.0, 0.0)
+    val sigma = DoubleVector(1.0, 1.0, 1.0)
+
+    val gaussian: (Point<Float64>) -> Double = { x ->
+        require(x.size == x0.size)
+        kotlin.math.exp(-((x - x0) / sigma).square().sum())
+    }
+
+    fun ((Point<Float64>) -> Double).grad(x: Point<Float64>): Point<Float64> {
+        require(x.size == x0.size)
+        return Float64Buffer(x.size) { i ->
+            val h = sigma[i] / 5
+            val dVector = Float64Buffer(x.size) { if (it == i) h else 0.0 }
+            val f1 = this(x + dVector / 2)
+            val f0 = this(x - dVector / 2)
+            (f1 - f0) / h
+        }
+    }
+
+    println(gaussian.grad(x0))
+
+}

examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt

@@ -0,0 +1,27 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import kotlin.random.Random
+import kotlin.time.measureTime
+
+fun main(): Unit = with(Float64LinearSpace) {
+    val random = Random(1224)
+    val dim = 500
+
+    //creating invertible matrix
+    val u = buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+    val l = buildMatrix(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
+    val matrix = l dot u
+
+    val time = measureTime {
+        repeat(20) {
+            lupSolver().inverse(matrix)
+        }
+    }
+
+    println(time)
+}