4622 changed files with 72038 additions and 587018 deletions
--- a/.github/CODEOWNERS
+++ b/.github/CODEOWNERS
@ -0,0 +1,3 @@
+@altavir
+
+/kmath-trajectory @ESchouten
--- a/.github/workflows/build.yml
+++ b/.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
--- a/.github/workflows/pages.yml
+++ b/.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
--- a/.github/workflows/publish.yml
+++ b/.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 }}
--- a/.gitignore
+++ b/.gitignore
@ -0,0 +1,24 @@
+.gradle
+build/
+out/
+
+.idea/
+.vscode/
+.fleet/
+.kotlin/
+
+
+# Avoid ignoring Gradle wrapper jar file (.jar files are usually ignored)
+!gradle-wrapper.jar
+
+# Cache of project
+.gradletasknamecache
+
+# Generated by javac -h and runtime
+*.class
+*.log
+
+!/.idea/copyright/
+!/.idea/scopes/
+/gradle/yarn.lock
+
--- a/.idea/copyright/kmath.xml
+++ b/.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>
--- a/.idea/copyright/profiles_settings.xml
+++ b/.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>
--- a/.idea/scopes/Apply_copyright.xml
+++ b/.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>
--- a/.space.kts
+++ b/.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\"",
+            )
+        }
+    }
+}
--- a/.space/CODEOWNERS
+++ b/.space/CODEOWNERS
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@ -0,0 +1,272 @@
+# KMath
+
+## 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)
+- 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)
+- New space extensions: `average` and `averageWith`
+- Local coding conventions
+- Geometric Domains API in `kmath-core`
+- Blocking chains in `kmath-coroutines`
+- Full hyperbolic functions support and default implementations within `ExtendedField`
+- 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
+- `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)
+- Gradle version: 6.3 -> 6.6
+- 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)
+- Commons RNG compatibility (https://github.com/mipt-npm/kmath/issues/93)
+- Multiplication of BigInt by scalar
--- a/README.md
+++ b/README.md
@ -0,0 +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/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
+
+
+### [attributes-kt](attributes-kt)
+> An API and basic implementation for arranging objects in a continuous memory block.
+>
+> **Maturity**: DEVELOPMENT
+
+### [benchmarks](benchmarks)
+>
+> **Maturity**: EXPERIMENTAL
+
+### [examples](examples)
+>
+> **Maturity**: EXPERIMENTAL
+
+### [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
+
+
+### [kmath-commons](kmath-commons)
+> Commons math binding for kmath
+>
+> **Maturity**: EXPERIMENTAL
+
+### [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
+
+
+### [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
+
+
+## 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("space.kscience:kmath-core:$version")
+    // api("space.kscience: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.
--- a/attributes-kt/README.md
+++ b/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")
+}
+```
--- a/attributes-kt/api/attributes-kt.api
+++ b/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;
+}
+
--- a/attributes-kt/build.gradle.kts
+++ b/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()
+}
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt
+++ b/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>>
+
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/AttributeContainer.kt
+++ b/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>
+
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt
+++ b/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
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/AttributesBuilder.kt
+++ b/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()
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/PolymorphicAttribute.kt
+++ b/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)
+}
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/SafeType.kt
+++ b/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>
+}
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/annotations.kt
+++ b/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
--- a/benchmarks/README.md
+++ b/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|
+
+
+
--- a/benchmarks/build.gradle.kts
+++ b/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}|")
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
--- a/benchmarks/docs/README-TEMPLATE.md
+++ b/benchmarks/docs/README-TEMPLATE.md
@ -0,0 +1,5 @@
+# BenchmarksResult
+
+${benchmarkMain}
+
+
--- a/benchmarks/src/jsMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt
+++ b/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)
+        }
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ArrayBenchmark.kt
+++ b/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) }
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt
+++ b/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))
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BufferBenchmark.kt
+++ b/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() }
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt
+++ b/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)
+    }
+
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt
+++ b/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)
+        }
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/IntegrationBenchmark.kt
+++ b/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)
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/JafamaBenchmark.kt
+++ b/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())) }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt
+++ b/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())
+    }
+
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt
+++ b/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)
+//    }
+
+
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/TensorAlgebraBenchmark.kt
+++ b/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))
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt
+++ b/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)
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt
+++ b/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)
+    }
+}
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/globals.kt
+++ b/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())
--- a/build.gradle.kts
+++ b/build.gradle.kts
@ -0,0 +1,72 @@
+import space.kscience.gradle.useApache2Licence
+import space.kscience.gradle.useSPCTeam
+
+plugins {
+    alias(spclibs.plugins.kscience.project)
+    alias(spclibs.plugins.kotlinx.kover)
+}
+
+val attributesVersion by extra("0.2.0")
+
+allprojects {
+    repositories {
+        maven("https://repo.kotlin.link")
+        maven("https://oss.sonatype.org/content/repositories/snapshots")
+        mavenCentral()
+    }
+
+    group = "space.kscience"
+    version = "0.4.1-dev"
+}
+
+subprojects {
+    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")
--- a/docs/algebra.md
+++ b/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`.
--- a/docs/buffers.md
+++ b/docs/buffers.md
@ -0,0 +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:
+
+* 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).
+
+## Buffer performance
+
+One should avoid using default boxing buffer wherever it is possible. Try to use primitive buffers or memory buffers
+instead.
--- a/docs/codestyle.md
+++ b/docs/codestyle.md
@ -0,0 +1,35 @@
+# Coding Conventions
+
+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.
+
+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
+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
+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
+cleanly separated.
+
+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.
--- a/docs/contexts.md
+++ b/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
+```
--- a/docs/diagrams/core.puml
+++ b/docs/diagrams/core.puml
--- a/docs/expressions.md
+++ b/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`&mdash;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**
--- a/docs/histograms.md
+++ b/docs/histograms.md
@ -0,0 +1 @@
+**TODO**
--- a/docs/images/KM.svg
+++ b/docs/images/KM.svg
--- a/docs/images/KM_mono.svg
+++ b/docs/images/KM_mono.svg
--- a/docs/images/KMath.svg
+++ b/docs/images/KMath.svg
--- a/docs/images/KMath_mono.svg
+++ b/docs/images/KMath_mono.svg
--- a/docs/linear.md
+++ b/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
--- a/docs/nd-structure.md
+++ b/docs/nd-structure.md
@ -0,0 +1,150 @@
+# ND-structure generation and operations
+
+**TODO**
+
+# Performance for n-dimensional structures operations
+
+One of the most sought after features of mathematical libraries is the high-performance operations on n-dimensional
+structures. In `kmath` performance depends on which particular context was used for operation.
+
+Let us consider following contexts:
+
+```kotlin
+    // automatically build context most suited for given type.
+    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(DoubleField, dim, dim)
+```
+
+Now let us perform several tests and see, which implementation is best suited for each case:
+
+## Test case
+
+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<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 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
+        repeat(n) {
+            res += 1.0
+        }
+    }
+```
+
+Ths speed of this operation is approximately the same as for specialized case since `NDField.auto` just
+returns the same `RealNDField` in this case. Of course, it is usually better to use specialized method to be sure.
+
+## Lazy
+
+Lazy field does not produce a structure when asked, instead it generates an empty structure and fills it on-demand
+using coroutines to parallelize computations.
+When one calls
+
+```kotlin
+    lazyField.run {
+        var res = one
+        repeat(n) {
+            res += 1.0
+        }
+    }
+```
+
+The result will be calculated almost immediately but the result will be empty. To get the full result
+structure one needs to call all its elements. In this case computation overhead will be huge. So this field never
+should be used if one expects to use the full result structure. Though if one wants only small fraction, it could
+save a lot of time.
+
+This field still could be used with reasonable performance if call code is changed:
+
+```kotlin
+    lazyField.run {
+        val res = one.map {
+            var c = 0.0
+            repeat(n) {
+                c += 1.0
+            }
+            c
+        }
+
+        res.elements().forEach { it.second }
+    }
+```
+
+In this case it completes in about `4x-5x` time due to boxing.
+
+## Boxing
+
+The boxing field produced by
+
+```kotlin
+    genericField.run {
+        var res: NDBuffer<Float64> = one
+        repeat(n) {
+            res += 1.0
+        }
+    }
+```
+
+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.
+
+## What about python?
+
+Usually it is bad idea to compare the direct numerical operation performance in different languages, but it hard to
+work completely without frame of reference. In this case, simple numpy code:
+
+```python
+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
+because better memory management). Of course if one writes `res += 1.0`, the performance will be different,
+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).
--- a/docs/polynomials.md
+++ b/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.
--- a/docs/readme.md
+++ b/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
--- a/docs/templates/ARTIFACT-TEMPLATE.md
+++ b/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}")
+}
+```
--- a/docs/templates/README-TEMPLATE.md
+++ b/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.
--- a/examples/README.md
+++ b/examples/README.md
@ -0,0 +1,4 @@
+# Module examples
+
+
+
--- a/examples/build.gradle.kts
+++ b/examples/build.gradle.kts
@ -0,0 +1,75 @@
+import org.jetbrains.kotlin.gradle.tasks.KotlinJvmCompile
+
+plugins {
+    kotlin("jvm")
+}
+
+repositories {
+    mavenCentral()
+    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-complex"))
+    implementation(project(":kmath-functions"))
+    implementation(project(":kmath-optimization"))
+    implementation(project(":kmath-stat"))
+    implementation(project(":kmath-viktor"))
+    implementation(project(":kmath-dimensions"))
+    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")
+}
+
+kotlin {
+    jvmToolchain(11)
+    sourceSets.all {
+        languageSettings {
+            optIn("kotlin.contracts.ExperimentalContracts")
+            optIn("kotlin.ExperimentalUnsignedTypes")
+            optIn("space.kscience.kmath.UnstableKMathAPI")
+        }
+    }
+}
+
+tasks.withType<KotlinJvmCompile> {
+    compilerOptions {
+        freeCompilerArgs.addAll("-Xjvm-default=all", "-Xopt-in=kotlin.RequiresOptIn", "-Xlambdas=indy")
+    }
+}
+
+readme {
+    maturity = space.kscience.gradle.Maturity.EXPERIMENTAL
+}
--- a/examples/notebooks/Naive
+++ b/examples/notebooks/Naive
@ -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
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/ast/astRendering.kt
+++ b/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)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/ast/expressions.kt
+++ b/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)
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt
+++ b/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))
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/ast/symjaSupport.kt
+++ b/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))
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/expressions/autodiff.kt
+++ b/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.
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt
+++ b/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()
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt
+++ b/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()
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt
+++ b/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
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt
+++ b/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()
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt
+++ b/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()
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt
+++ b/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)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/jafama/JafamaDemo.kt
+++ b/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) })
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt
+++ 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)
+
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/eigenValueDecomposition.kt
+++ b/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" )
+        }
+    }
+
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt
+++ b/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))
+
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt
+++ b/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)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/operations/BigIntDemo.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/operations/BigIntDemo.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.operations
+
+fun main() {
+    val res = BigIntField { number(1) * 2 }
+    println("bigint:$res")
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.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.operations
+
+import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.complex.algebra
+import space.kscience.kmath.complex.ndAlgebra
+import space.kscience.kmath.nd.BufferND
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.structureND
+
+fun main() = Complex.algebra {
+    val complex = 2 + 2 * i
+    println(complex * 8 - 5 * i)
+
+    //flat buffer
+    val buffer = with(bufferAlgebra) {
+        buffer(8) { Complex(it, -it) }.map { Complex(it.im, it.re) }
+    }
+    println(buffer)
+
+    // 2d element
+    val element: BufferND<Complex> = ndAlgebra.structureND(2, 2) { (i, j) ->
+        Complex(i - j, i + j)
+    }
+    println(element)
+
+    // 1d element operation
+    val result: StructureND<Complex> = ndAlgebra {
+        val a = structureND(8) { (it) -> i * it - it.toDouble() }
+        val b = 3
+        val c = Complex(1.0, 1.0)
+
+        (a pow b) + c
+    }
+    println(result)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.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.operations
+
+import space.kscience.kmath.commons.linear.CMLinearSpace
+import space.kscience.kmath.linear.matrix
+import space.kscience.kmath.nd.Float64BufferND
+import space.kscience.kmath.nd.Structure2D
+import space.kscience.kmath.nd.mutableStructureND
+import space.kscience.kmath.nd.ndAlgebra
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.viktor.viktorAlgebra
+import kotlin.collections.component1
+import kotlin.collections.component2
+
+fun main() {
+    val viktorStructure = Float64Field.viktorAlgebra.mutableStructureND(2, 2) { (i, j) ->
+        if (i == j) 2.0 else 0.0
+    }
+
+    val cmMatrix: Structure2D<Float64> = CMLinearSpace.matrix(2, 2)(0.0, 1.0, 0.0, 3.0)
+
+    val res: Float64BufferND = Float64Field.ndAlgebra {
+        exp(viktorStructure) + 2.0 * cmMatrix
+    }
+
+    println(res)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/series/DateTimeSeries.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/series/DateTimeSeries.kt
@ -0,0 +1,17 @@
+/*
+ * 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.series
+
+import kotlinx.datetime.Instant
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import kotlin.time.Duration
+
+fun SeriesAlgebra.Companion.time(zero: Instant, step: Duration) = MonotonicSeriesAlgebra(
+    bufferAlgebra = Double.algebra.bufferAlgebra,
+    offsetToLabel = { zero + step * it },
+    labelToOffset = { (it - zero) / step }
+)
--- a/examples/src/main/kotlin/space/kscience/kmath/series/analyzeDif.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/series/analyzeDif.kt
@ -0,0 +1,65 @@
+/*
+ * 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.series
+
+
+import kotlinx.html.h1
+import kotlinx.html.p
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.toList
+import space.kscience.kmath.stat.KMComparisonResult
+import space.kscience.kmath.stat.ksComparisonStatistic
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.structures.slice
+import space.kscience.plotly.*
+import kotlin.math.PI
+
+fun Double.Companion.seriesAlgebra() = Double.algebra.bufferAlgebra.seriesAlgebra()
+
+
+fun main() = with(Double.seriesAlgebra()) {
+
+
+    fun Plot.plotSeries(name: String, buffer: Buffer<Float64>) {
+        scatter {
+            this.name = name
+            x.numbers = buffer.labels
+            y.numbers = buffer.toList()
+        }
+    }
+
+
+    val s1 = series(100) { sin(2 * PI * it / 100) + 1.0 }
+
+    val s2 = s1.slice(20..50).moveTo(40)
+
+    val s3: Buffer<Float64> = s1.zip(s2) { l, r -> l + r } //s1 + s2
+    val s4 = s3.map { ln(it) }
+
+    val kmTest: KMComparisonResult<Float64> = ksComparisonStatistic(s1, s2)
+
+    Plotly.page {
+        h1 { +"This is my plot" }
+        p {
+            +"Kolmogorov-smirnov test for s1 and s2: ${kmTest.value}"
+        }
+        plot {
+            plotSeries("s1", s1)
+            plotSeries("s2", s2)
+            plotSeries("s3", s3)
+            plotSeries("s4", s4)
+            layout {
+                xaxis {
+                    range(0.0..100.0)
+                }
+            }
+        }
+
+    }.makeFile()
+
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt
@ -0,0 +1,47 @@
+/*
+ * 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.kmath.series
+
+
+import space.kscience.kmath.structures.*
+import space.kscience.plotly.*
+import space.kscience.plotly.models.Scatter
+import space.kscience.plotly.models.ScatterMode
+import kotlin.random.Random
+
+fun main(): Unit = with(Double.seriesAlgebra()) {
+
+    val random = Random(1234)
+
+    val arrayOfRandoms = DoubleArray(20) { random.nextDouble() }
+
+    val series1: Float64Buffer = arrayOfRandoms.asBuffer()
+    val series2: Series<Float64> = series1.moveBy(3)
+
+    val res = series2 - series1
+
+    println(res.size)
+
+    println(res)
+
+    fun Plot.series(name: String, buffer: Buffer<Float64>, block: Scatter.() -> Unit = {}) {
+        scatter {
+            this.name = name
+            x.numbers = buffer.offsetIndices
+            y.doubles = buffer.toDoubleArray()
+            block()
+        }
+    }
+
+    Plotly.plot {
+        series("series1", series1)
+        series("series2", series2)
+        series("dif", res) {
+            mode = ScatterMode.lines
+            line.color("magenta")
+        }
+    }.makeFile(resourceLocation = ResourceLocation.REMOTE)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt
@ -0,0 +1,94 @@
+/*
+ * 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.stat
+
+import kotlinx.coroutines.Dispatchers
+import kotlinx.coroutines.async
+import kotlinx.coroutines.runBlocking
+import org.apache.commons.rng.sampling.distribution.BoxMullerNormalizedGaussianSampler
+import org.apache.commons.rng.simple.RandomSource
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.samplers.GaussianSampler
+import java.time.Duration
+import java.time.Instant
+import org.apache.commons.rng.sampling.distribution.GaussianSampler as CMGaussianSampler
+
+private suspend fun runKMathChained(): Duration {
+    val generator = RandomGenerator.fromSource(RandomSource.MT, 123L)
+    val normal = GaussianSampler(7.0, 2.0)
+    val chain = normal.sample(generator)
+    val startTime = Instant.now()
+    var sum = 0.0
+
+    repeat(10000001) { counter ->
+        sum += chain.next()
+
+        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 runKMathBlocking(): Duration {
+    val generator = RandomGenerator.fromSource(RandomSource.MT, 123L)
+    val normal = GaussianSampler(7.0, 2.0)
+    val chain = normal.sample(generator)
+    val startTime = Instant.now()
+    var sum = 0.0
+
+    repeat(10000001) { counter ->
+        sum += chain.nextBlocking()
+
+        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 runCMDirect(): Duration {
+    val rng = RandomSource.MT.create(123L)
+
+    val sampler = CMGaussianSampler.of(
+        BoxMullerNormalizedGaussianSampler.of(rng),
+        7.0,
+        2.0
+    )
+
+    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(): Unit = runBlocking(Dispatchers.Default) {
+    val directJob = async { runCMDirect() }
+    val chainJob = async { runKMathChained() }
+    val blockingJob = async { runKMathBlocking() }
+    println("KMath Chained: ${chainJob.await()}")
+    println("KMath Blocking: ${blockingJob.await()}")
+    println("Apache Direct: ${directJob.await()}")
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.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.stat
+
+import kotlinx.coroutines.runBlocking
+import space.kscience.kmath.chains.Chain
+import space.kscience.kmath.chains.combineWithState
+import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.structures.Float64
+
+private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
+
+/**
+ * Averaging.
+ */
+private fun Chain<Float64>.mean(): Chain<Float64> = combineWithState(AveragingChainState(), { it.copy() }) { chain ->
+    val next = chain.next()
+    num++
+    value += next
+    return@combineWithState value / num
+}
+
+
+fun main() {
+    val normal = NormalDistribution(0.0, 2.0)
+    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")
+            }
+        }
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt
@ -0,0 +1,66 @@
+/*
+ * 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.
+ */
+
+@file:Suppress("unused")
+
+package space.kscience.kmath.structures
+
+import space.kscience.kmath.complex.*
+import space.kscience.kmath.linear.transposed
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.ndAlgebra
+import space.kscience.kmath.nd.structureND
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.invoke
+import kotlin.system.measureTimeMillis
+
+fun main() {
+    val dim = 1000
+    val n = 1000
+
+    val realField = Float64Field.ndAlgebra(dim, dim)
+    val complexField: ComplexFieldND = ComplexField.ndAlgebra(dim, dim)
+
+    val realTime = measureTimeMillis {
+        realField {
+            var res: StructureND<Float64> = one
+            repeat(n) {
+                res += 1.0
+            }
+        }
+    }
+
+    println("Real addition completed in $realTime millis")
+
+    val complexTime = measureTimeMillis {
+        complexField {
+            var res: StructureND<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 {
+        withNdAlgebra(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 = structureND { (k, l) -> k + l * i }
+            //Perform sum
+            val sum = matrix + x + 1.0
+
+            //Represent the sum as 2d-structure and transpose
+            sum.as2D().transposed()
+        }
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.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.structures
+
+import kotlinx.coroutines.DelicateCoroutinesApi
+import kotlinx.coroutines.GlobalScope
+import org.nd4j.linalg.factory.Nd4j
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.nd4j.nd4j
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.viktor.ViktorFieldND
+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")
+}
+
+@OptIn(DelicateCoroutinesApi::class)
+fun main() {
+    // initializing Nd4j
+    Nd4j.zeros(0)
+    val dim = 1000
+    val n = 1000
+    val shape = ShapeND(dim, dim)
+
+
+    // specialized nd-field for Double. It works as generic Double field as well.
+    val doubleField = Float64Field.ndAlgebra
+    //A generic field. It should be used for objects, not primitives.
+    val genericField = BufferedFieldOpsND(Float64Field)
+    // Nd4j specialized field.
+    val nd4jField = Float64Field.nd4j
+    //viktor field
+    val viktorField = ViktorFieldND(dim, dim)
+    //parallel processing based on Java Streams
+    val parallelField = Float64Field.ndStreaming(dim, dim)
+
+    measureAndPrint("Boxing addition") {
+        genericField {
+            var res: StructureND<Float64> = one(shape)
+            repeat(n) { res += 1.0 }
+        }
+    }
+
+    measureAndPrint("Specialized addition") {
+        doubleField {
+            var res: StructureND<Float64> = one(shape)
+            repeat(n) { res += 1.0 }
+        }
+    }
+
+    measureAndPrint("Nd4j specialized addition") {
+        nd4jField {
+            var res: StructureND<Float64> = one(shape)
+            repeat(n) { res += 1.0 }
+        }
+    }
+
+    measureAndPrint("Viktor addition") {
+        viktorField {
+            var res: StructureND<Float64> = one
+            repeat(n) { res += 1.0 }
+        }
+    }
+
+    measureAndPrint("Parallel stream addition") {
+        parallelField {
+            var res: StructureND<Float64> = one
+            repeat(n) { res += 1.0 }
+        }
+    }
+
+    measureAndPrint("Lazy addition") {
+        val res = doubleField.one(shape).mapAsync(GlobalScope) {
+            var c = 0.0
+            repeat(n) {
+                c += 1.0
+            }
+            c
+        }
+
+        res.elements().forEach { it.second }
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/StreamDoubleFieldND.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/StreamDoubleFieldND.kt
@ -0,0 +1,127 @@
+/*
+ * 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.structures
+
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.NumbersAddOps
+import java.util.*
+import java.util.stream.IntStream
+
+/**
+ * A demonstration implementation of NDField over Real using Java [java.util.stream.DoubleStream] for parallel
+ * execution.
+ */
+class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, Float64Field>,
+    NumbersAddOps<StructureND<Float64>>,
+    ExtendedField<StructureND<Float64>> {
+
+    private val strides = ColumnStrides(shape)
+    override val elementAlgebra: Float64Field get() = Float64Field
+    override val zero: BufferND<Float64> by lazy { structureND(shape) { zero } }
+    override val one: BufferND<Float64> by lazy { structureND(shape) { one } }
+
+    override fun number(value: Number): BufferND<Float64> {
+        val d = value.toDouble() // minimize conversions
+        return structureND(shape) { d }
+    }
+
+    @OptIn(PerformancePitfall::class)
+    private val StructureND<Float64>.buffer: Float64Buffer
+        get() = when {
+            shape != this@StreamDoubleFieldND.shape -> throw ShapeMismatchException(
+                this@StreamDoubleFieldND.shape,
+                shape
+            )
+
+            this is BufferND && indices == this@StreamDoubleFieldND.strides -> this.buffer as Float64Buffer
+            else -> Float64Buffer(strides.linearSize) { offset -> get(strides.index(offset)) }
+        }
+
+    override fun structureND(shape: ShapeND, initializer: Float64Field.(IntArray) -> Double): BufferND<Float64> {
+        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
+            val index = strides.index(offset)
+            Float64Field.initializer(index)
+        }.toArray()
+
+        return BufferND(strides, array.asBuffer())
+    }
+
+    override fun mutableStructureND(
+        shape: ShapeND,
+        initializer: DoubleField.(IntArray) -> Double,
+    ): MutableBufferND<Float64> {
+        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
+            val index = strides.index(offset)
+            DoubleField.initializer(index)
+        }.toArray()
+
+        return MutableBufferND(strides, array.asBuffer())
+    }
+
+    @OptIn(PerformancePitfall::class)
+    override fun StructureND<Float64>.map(
+        transform: Float64Field.(Double) -> Double,
+    ): BufferND<Float64> {
+        val array = Arrays.stream(buffer.array).parallel().map { Float64Field.transform(it) }.toArray()
+        return BufferND(strides, array.asBuffer())
+    }
+
+    @OptIn(PerformancePitfall::class)
+    override fun StructureND<Float64>.mapIndexed(
+        transform: Float64Field.(index: IntArray, Double) -> Double,
+    ): BufferND<Float64> {
+        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
+            Float64Field.transform(
+                strides.index(offset),
+                buffer.array[offset]
+            )
+        }.toArray()
+
+        return BufferND(strides, array.asBuffer())
+    }
+
+    @OptIn(PerformancePitfall::class)
+    override fun zip(
+        left: StructureND<Float64>,
+        right: StructureND<Float64>,
+        transform: Float64Field.(Double, Double) -> Double,
+    ): BufferND<Float64> {
+        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
+            Float64Field.transform(left.buffer.array[offset], right.buffer.array[offset])
+        }.toArray()
+        return BufferND(strides, array.asBuffer())
+    }
+
+    override fun StructureND<Float64>.unaryMinus(): StructureND<Float64> = map { -it }
+
+    override fun scale(a: StructureND<Float64>, value: Double): StructureND<Float64> = a.map { it * value }
+
+    override fun power(arg: StructureND<Float64>, pow: Number): BufferND<Float64> = arg.map { power(it, pow) }
+
+    override fun exp(arg: StructureND<Float64>): BufferND<Float64> = arg.map { exp(it) }
+
+    override fun ln(arg: StructureND<Float64>): BufferND<Float64> = arg.map { ln(it) }
+
+    override fun sin(arg: StructureND<Float64>): BufferND<Float64> = arg.map { sin(it) }
+    override fun cos(arg: StructureND<Float64>): BufferND<Float64> = arg.map { cos(it) }
+    override fun tan(arg: StructureND<Float64>): BufferND<Float64> = arg.map { tan(it) }
+    override fun asin(arg: StructureND<Float64>): BufferND<Float64> = arg.map { asin(it) }
+    override fun acos(arg: StructureND<Float64>): BufferND<Float64> = arg.map { acos(it) }
+    override fun atan(arg: StructureND<Float64>): BufferND<Float64> = arg.map { atan(it) }
+
+    override fun sinh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { sinh(it) }
+    override fun cosh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { cosh(it) }
+    override fun tanh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { tanh(it) }
+    override fun asinh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { asinh(it) }
+    override fun acosh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { acosh(it) }
+    override fun atanh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { atanh(it) }
+}
+
+fun Float64Field.ndStreaming(vararg shape: Int): StreamDoubleFieldND = StreamDoubleFieldND(ShapeND(shape))
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/StructureReadBenchmark.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/StructureReadBenchmark.kt
@ -0,0 +1,45 @@
+/*
+ * 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.structures
+
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.nd.BufferND
+import space.kscience.kmath.nd.ColumnStrides
+import space.kscience.kmath.nd.ShapeND
+import kotlin.system.measureTimeMillis
+
+@Suppress("ASSIGNED_BUT_NEVER_ACCESSED_VARIABLE")
+@OptIn(PerformancePitfall::class)
+fun main() {
+    val n = 6000
+    val array = DoubleArray(n * n) { 1.0 }
+    val buffer = Float64Buffer(array)
+    val strides = ColumnStrides(ShapeND(n, n))
+    val structure = BufferND(strides, buffer)
+
+    measureTimeMillis {
+        var res = 0.0
+        strides.asSequence().forEach { res = structure[it] }
+    } // warmup
+
+    val time1 = measureTimeMillis {
+        var res = 0.0
+        strides.asSequence().forEach { res = structure[it] }
+    }
+    println("Structure reading finished in $time1 millis")
+
+    val time2 = measureTimeMillis {
+        var res = 0.0
+        strides.asSequence().forEach { res = buffer[strides.offset(it)] }
+    }
+    println("Buffer reading finished in $time2 millis")
+
+    val time3 = measureTimeMillis {
+        var res = 0.0
+        strides.asSequence().forEach { res = array[strides.offset(it)] }
+    }
+    println("Array reading finished in $time3 millis")
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/StructureWriteBenchmark.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/StructureWriteBenchmark.kt
@ -0,0 +1,42 @@
+/*
+ * 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.structures
+
+import space.kscience.kmath.nd.BufferND
+import space.kscience.kmath.operations.mapToBuffer
+import kotlin.system.measureTimeMillis
+
+private inline fun <T, reified R : Any> BufferND<T>.mapToBufferND(
+    bufferFactory: BufferFactory<R> = BufferFactory(),
+    crossinline block: (T) -> R,
+): BufferND<R> = BufferND(indices, buffer.mapToBuffer(bufferFactory, block))
+
+@Suppress("UNUSED_VARIABLE")
+fun main() {
+    val n = 6000
+    val structure = BufferND(n, n) { 1.0 }
+    structure.mapToBufferND { it + 1 } // warm-up
+    val time1 = measureTimeMillis { val res = structure.mapToBufferND { 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 = Float64Buffer(DoubleArray(n * n) { 1.0 })
+
+    val time3 = measureTimeMillis {
+        val target = Float64Buffer(DoubleArray(n * n))
+        val res = array.forEachIndexed { index, value ->
+            target[index] = value + 1
+        }
+    }
+    println("Buffer mapping finished in $time3 millis")
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt
@ -0,0 +1,23 @@
+/*
+ * 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.structures
+
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.buffer
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.withSize
+
+inline fun <reified R : Any> MutableBuffer.Companion.same(
+    n: Int,
+    value: R,
+): MutableBuffer<R> = MutableBuffer(n) { value }
+
+
+fun main() {
+    with(Float64Field.bufferAlgebra.withSize(5)) {
+        println(number(2.0) + buffer(1, 2, 3, 4, 5))
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/mutableNd.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/mutableNd.kt
@ -0,0 +1,26 @@
+/*
+ * 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.kmath.structures
+
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.operations.algebra
+
+@OptIn(PerformancePitfall::class)
+fun main(): Unit = with(Double.algebra.ndAlgebra) {
+    val structure: MutableStructure2D<Float64> = mutableStructureND(ShapeND(2, 2)) { (i, j) ->
+        i.toDouble() + j.toDouble()
+    }.as2D()
+
+    structure[0, 1] = -2.0
+
+    val structure2 = mutableStructureND(2, 2) { (i, j) -> i.toDouble() + j.toDouble() }.as2D()
+
+    structure2[0, 1] = 2.0
+
+
+    println(structure + structure2)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/typeSafeDimensions.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/typeSafeDimensions.kt
@ -0,0 +1,35 @@
+/*
+ * 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.structures
+
+import space.kscience.kmath.dimensions.D2
+import space.kscience.kmath.dimensions.D3
+import space.kscience.kmath.dimensions.DMatrixContext
+import space.kscience.kmath.dimensions.Dimension
+
+private fun DMatrixContext<Double, *>.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.transposed() + m2
+}
+
+private object D5 : Dimension {
+    override val dim: Int = 5
+}
+
+private fun DMatrixContext<Double, *>.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()
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.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.tensors.LevenbergMarquardt.StaticLm
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.math.roundToInt
+
+fun main() {
+    val NData = 200
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1) - 104
+    }
+
+    val Nparams = 15
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat = funcDifficultForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 0.9)
+    }
+
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0 / Nparams * 1.0 - 0.085
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-2, 11.0, 9.0, 1.0)
+//    val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-3, 11.0, 9.0, 1.0)
+
+    val inputData = LMInput(
+        ::funcDifficultForLm,
+        p_init.as2D(),
+        t,
+        y_dat,
+        weight,
+        dp,
+        p_min.as2D(),
+        p_max.as2D(),
+        opts[1].toInt(),
+        doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+        doubleArrayOf(opts[6], opts[7], opts[8]),
+        opts[9].toInt(),
+        10,
+        1
+    )
+
+    val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+    println("Parameters:")
+    for (i in 0 until result.resultParameters.shape.component1()) {
+        val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0
+        print("$x ")
+    }
+    println()
+
+    println("Y true and y received:")
+    var y_hat_after = funcDifficultForLm(t_example, result.resultParameters, exampleNumber)
+    for (i in 0 until y_hat.shape.component1()) {
+        val x = (y_hat[i, 0] * 10000).roundToInt() / 10000.0
+        val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0
+        println("$x $y")
+    }
+
+    println("Сhi_sq:")
+    println(result.resultChiSq)
+    println("Number of iterations:")
+    println(result.iterations)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.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.tensors.LevenbergMarquardt.StaticLm
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcEasyForLm
+import space.kscience.kmath.tensors.LevenbergMarquardt.getStartDataForFuncEasy
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.math.roundToInt
+
+fun main() {
+    val startedData = getStartDataForFuncEasy()
+    val inputData = LMInput(
+        ::funcEasyForLm,
+        DoubleTensorAlgebra.ones(ShapeND(intArrayOf(4, 1))).as2D(),
+        startedData.t,
+        startedData.y_dat,
+        startedData.weight,
+        startedData.dp,
+        startedData.p_min,
+        startedData.p_max,
+        startedData.opts[1].toInt(),
+        doubleArrayOf(startedData.opts[2], startedData.opts[3], startedData.opts[4], startedData.opts[5]),
+        doubleArrayOf(startedData.opts[6], startedData.opts[7], startedData.opts[8]),
+        startedData.opts[9].toInt(),
+        10,
+        startedData.example_number
+    )
+
+    val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+    println("Parameters:")
+    for (i in 0 until result.resultParameters.shape.component1()) {
+        val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0
+        print("$x ")
+    }
+    println()
+
+    println("Y true and y received:")
+    var y_hat_after = funcDifficultForLm(startedData.t, result.resultParameters, startedData.example_number)
+    for (i in 0 until startedData.y_dat.shape.component1()) {
+        val x = (startedData.y_dat[i, 0] * 10000).roundToInt() / 10000.0
+        val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0
+        println("$x $y")
+    }
+
+    println("Сhi_sq:")
+    println(result.resultChiSq)
+    println("Number of iterations:")
+    println(result.iterations)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt
@ -0,0 +1,91 @@
+/*
+ * 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.tensors.LevenbergMarquardt.StaticLm
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcMiddleForLm
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.math.roundToInt
+
+fun main() {
+    val NData = 100
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1)
+    }
+
+    val Nparams = 20
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat = funcMiddleForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 0.9)
+    }
+
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0)
+
+    val inputData = LMInput(
+        ::funcMiddleForLm,
+        p_init.as2D(),
+        t,
+        y_dat,
+        weight,
+        dp,
+        p_min.as2D(),
+        p_max.as2D(),
+        opts[1].toInt(),
+        doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+        doubleArrayOf(opts[6], opts[7], opts[8]),
+        opts[9].toInt(),
+        10,
+        1
+    )
+
+    val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+    println("Parameters:")
+    for (i in 0 until result.resultParameters.shape.component1()) {
+        val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0
+        print("$x ")
+    }
+    println()
+
+
+    var y_hat_after = funcMiddleForLm(t_example, result.resultParameters, exampleNumber)
+    for (i in 0 until y_hat.shape.component1()) {
+        val x = (y_hat[i, 0] * 10000).roundToInt() / 10000.0
+        val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0
+        println("$x $y")
+    }
+
+    println("Сhi_sq:")
+    println(result.resultChiSq)
+    println("Number of iterations:")
+    println(result.iterations)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt
@ -0,0 +1,76 @@
+/*
+ * 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.tensors.LevenbergMarquardt.StreamingLm
+
+import kotlinx.coroutines.delay
+import kotlinx.coroutines.flow.Flow
+import kotlinx.coroutines.flow.flow
+import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.tensors.LevenbergMarquardt.StartDataLm
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.zeros
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.random.Random
+
+fun streamLm(
+    lm_func: (MutableStructure2D<Float64>, MutableStructure2D<Float64>, Int) -> (MutableStructure2D<Float64>),
+    startData: StartDataLm, launchFrequencyInMs: Long, numberOfLaunches: Int,
+): Flow<MutableStructure2D<Float64>> = flow {
+
+    var example_number = startData.example_number
+    var p_init = startData.p_init
+    var t = startData.t
+    var y_dat = startData.y_dat
+    val weight = startData.weight
+    val dp = startData.dp
+    val p_min = startData.p_min
+    val p_max = startData.p_max
+    val opts = startData.opts
+
+    var steps = numberOfLaunches
+    val isEndless = (steps <= 0)
+
+    val inputData = LMInput(
+        lm_func,
+        p_init,
+        t,
+        y_dat,
+        weight,
+        dp,
+        p_min,
+        p_max,
+        opts[1].toInt(),
+        doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+        doubleArrayOf(opts[6], opts[7], opts[8]),
+        opts[9].toInt(),
+        10,
+        example_number
+    )
+
+    while (isEndless || steps > 0) {
+        val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+        emit(result.resultParameters)
+        delay(launchFrequencyInMs)
+        inputData.realValues = generateNewYDat(y_dat, 0.1)
+        inputData.startParameters = result.resultParameters
+        if (!isEndless) steps -= 1
+    }
+}
+
+fun generateNewYDat(y_dat: MutableStructure2D<Float64>, delta: Double): MutableStructure2D<Float64> {
+    val n = y_dat.shape.component1()
+    val y_dat_new = zeros(ShapeND(intArrayOf(n, 1))).as2D()
+    for (i in 0 until n) {
+        val randomEps = Random.nextDouble(delta + delta) - delta
+        y_dat_new[i, 0] = y_dat[i, 0] + randomEps
+    }
+    return y_dat_new
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.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.tensors.LevenbergMarquardt.StreamingLm
+
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm
+import space.kscience.kmath.tensors.LevenbergMarquardt.getStartDataForFuncDifficult
+import kotlin.math.roundToInt
+
+suspend fun main() {
+    val startData = getStartDataForFuncDifficult()
+    // Создание потока:
+    val lmFlow = streamLm(::funcDifficultForLm, startData, 0, 100)
+    var initialTime = System.currentTimeMillis()
+    var lastTime: Long
+    val launches = mutableListOf<Long>()
+    // Запуск потока
+    lmFlow.collect { parameters ->
+        lastTime = System.currentTimeMillis()
+        launches.add(lastTime - initialTime)
+        initialTime = lastTime
+        for (i in 0 until parameters.shape.component1()) {
+            val x = (parameters[i, 0] * 10000).roundToInt() / 10000.0
+            print("$x ")
+            if (i == parameters.shape.component1() - 1) println()
+        }
+    }
+
+    println("Average without first is: ${launches.subList(1, launches.size - 1).average()}")
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt
@ -0,0 +1,233 @@
+/*
+ * 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.tensors.LevenbergMarquardt
+
+import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.max
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.pow
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.times
+import space.kscience.kmath.tensors.core.asDoubleTensor
+
+public data class StartDataLm(
+    var lm_matx_y_dat: MutableStructure2D<Float64>,
+    var example_number: Int,
+    var p_init: MutableStructure2D<Float64>,
+    var t: MutableStructure2D<Float64>,
+    var y_dat: MutableStructure2D<Float64>,
+    var weight: Double,
+    var dp: MutableStructure2D<Float64>,
+    var p_min: MutableStructure2D<Float64>,
+    var p_max: MutableStructure2D<Float64>,
+    var consts: MutableStructure2D<Float64>,
+    var opts: DoubleArray,
+)
+
+fun funcEasyForLm(
+    t: MutableStructure2D<Float64>,
+    p: MutableStructure2D<Float64>,
+    exampleNumber: Int,
+): MutableStructure2D<Float64> {
+    val m = t.shape.component1()
+    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))
+
+    if (exampleNumber == 1) {
+        y_hat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))).times(p[0, 0]) + t.times(p[2, 0]).times(
+            DoubleTensorAlgebra.exp((t.times(-1.0 / p[3, 0])))
+        )
+    } else if (exampleNumber == 2) {
+        val mt = t.max()
+        y_hat = (t.times(1.0 / mt)).times(p[0, 0]) +
+                (t.times(1.0 / mt)).pow(2).times(p[1, 0]) +
+                (t.times(1.0 / mt)).pow(3).times(p[2, 0]) +
+                (t.times(1.0 / mt)).pow(4).times(p[3, 0])
+    } else if (exampleNumber == 3) {
+        y_hat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0])))
+            .times(p[0, 0]) + DoubleTensorAlgebra.sin((t.times(1.0 / p[3, 0]))).times(p[2, 0])
+    }
+
+    return y_hat.as2D()
+}
+
+fun funcMiddleForLm(
+    t: MutableStructure2D<Float64>,
+    p: MutableStructure2D<Float64>,
+    exampleNumber: Int,
+): MutableStructure2D<Float64> {
+    val m = t.shape.component1()
+    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))
+
+    val mt = t.max()
+    for (i in 0 until p.shape.component1()) {
+        y_hat += (t.times(1.0 / mt)).times(p[i, 0])
+    }
+
+    for (i in 0 until 5) {
+        y_hat = funcEasyForLm(y_hat.as2D(), p, exampleNumber).asDoubleTensor()
+    }
+
+    return y_hat.as2D()
+}
+
+fun funcDifficultForLm(
+    t: MutableStructure2D<Float64>,
+    p: MutableStructure2D<Float64>,
+    exampleNumber: Int,
+): MutableStructure2D<Float64> {
+    val m = t.shape.component1()
+    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))
+
+    val mt = t.max()
+    for (i in 0 until p.shape.component1()) {
+        y_hat = y_hat.plus((t.times(1.0 / mt)).times(p[i, 0]))
+    }
+
+    for (i in 0 until 4) {
+        y_hat = funcEasyForLm((y_hat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor()
+    }
+
+    return y_hat.as2D()
+}
+
+
+fun getStartDataForFuncDifficult(): StartDataLm {
+    val NData = 200
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1) - 104
+    }
+
+    val Nparams = 15
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat = funcDifficultForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 0.9)
+    }
+
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0 / Nparams * 1.0 - 0.085
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val consts = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0)
+    ).as2D()
+    val opts = doubleArrayOf(3.0, 10000.0, 1e-2, 1e-3, 1e-2, 1e-2, 1e-2, 11.0, 9.0, 1.0)
+
+    return StartDataLm(y_dat, 1, p_init, t, y_dat, weight, dp, p_min.as2D(), p_max.as2D(), consts, opts)
+}
+
+fun getStartDataForFuncMiddle(): StartDataLm {
+    val NData = 100
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1)
+    }
+
+    val Nparams = 20
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat = funcMiddleForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 10.0)
+    }
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val consts = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0)
+    ).as2D()
+    val opts = doubleArrayOf(3.0, 10000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0)
+
+    var example_number = 1
+
+    return StartDataLm(y_dat, example_number, p_init, t, y_dat, weight, dp, p_min.as2D(), p_max.as2D(), consts, opts)
+}
+
+fun getStartDataForFuncEasy(): StartDataLm {
+    val lm_matx_y_dat = doubleArrayOf(
+        19.6594, 18.6096, 17.6792, 17.2747, 16.3065, 17.1458, 16.0467, 16.7023, 15.7809, 15.9807,
+        14.7620, 15.1128, 16.0973, 15.1934, 15.8636, 15.4763, 15.6860, 15.1895, 15.3495, 16.6054,
+        16.2247, 15.9854, 16.1421, 17.0960, 16.7769, 17.1997, 17.2767, 17.5882, 17.5378, 16.7894,
+        17.7648, 18.2512, 18.1581, 16.7037, 17.8475, 17.9081, 18.3067, 17.9632, 18.2817, 19.1427,
+        18.8130, 18.5658, 18.0056, 18.4607, 18.5918, 18.2544, 18.3731, 18.7511, 19.3181, 17.3066,
+        17.9632, 19.0513, 18.7528, 18.2928, 18.5967, 17.8567, 17.7859, 18.4016, 18.9423, 18.4959,
+        17.8000, 18.4251, 17.7829, 17.4645, 17.5221, 17.3517, 17.4637, 17.7563, 16.8471, 17.4558,
+        17.7447, 17.1487, 17.3183, 16.8312, 17.7551, 17.0942, 15.6093, 16.4163, 15.3755, 16.6725,
+        16.2332, 16.2316, 16.2236, 16.5361, 15.3721, 15.3347, 15.5815, 15.6319, 14.4538, 14.6044,
+        14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157
+    )
+
+    var example_number = 1
+    val p_init = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0)
+    ).as2D()
+
+    var t = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(100, 1))).as2D()
+    for (i in 0 until 100) {
+        t[i, 0] = t[i, 0] * (i + 1)
+    }
+
+    val y_dat = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(100, 1)), lm_matx_y_dat
+    ).as2D()
+
+    val weight = 4.0
+
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+
+    val p_min = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(4, 1)), doubleArrayOf(-50.0, -20.0, -2.0, -100.0)
+    ).as2D()
+
+    val p_max = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0)
+    ).as2D()
+
+    val consts = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0)
+    ).as2D()
+
+    val opts = doubleArrayOf(3.0, 100.0, 1e-3, 1e-3, 1e-1, 1e-1, 1e-2, 11.0, 9.0, 1.0)
+
+    return StartDataLm(y_dat, example_number, p_init, t, y_dat, weight, dp, p_min, p_max, consts, opts)
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/OLSWithSVD.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/OLSWithSVD.kt
@ -0,0 +1,72 @@
+/*
+ * 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.tensors
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.tensors.core.DoubleTensor
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.randomNormal
+import space.kscience.kmath.tensors.core.randomNormalLike
+import kotlin.math.abs
+
+// OLS estimator using SVD
+
+fun main() {
+    //seed for random
+    val randSeed = 100500L
+
+    // work in context with linear operations
+    DoubleTensorAlgebra {
+        // take coefficient vector from normal distribution
+        val alpha = randomNormal(
+            ShapeND(5),
+            randSeed
+        ) + fromArray(
+            ShapeND(5),
+            doubleArrayOf(1.0, 2.5, 3.4, 5.0, 10.1)
+        )
+
+        println("Real alpha:\n$alpha")
+
+        // also take sample of size 20 from normal distribution for x
+        val x = randomNormal(
+            ShapeND(20, 5),
+            randSeed
+        )
+
+        // calculate y and add gaussian noise (N(0, 0.05))
+        val y = x dot alpha
+        y += randomNormalLike(y, randSeed) * 0.05
+
+        // now restore the coefficient vector with OSL estimator with SVD
+        val (u, singValues, v) = svd(x)
+
+        // we have to make sure the singular values of the matrix are not close to zero
+        println("Singular values:\n$singValues")
+
+
+        // inverse Sigma matrix can be restored from singular values with diagonalEmbedding function
+        val sigma = diagonalEmbedding(singValues.map { if (abs(it) < 1e-3) 0.0 else 1.0 / it })
+
+        val alphaOLS = v dot sigma dot u.transposed() dot y
+        println(
+            "Estimated alpha:\n" +
+                    "$alphaOLS"
+        )
+
+        // figure out MSE of approximation
+        fun mse(yTrue: DoubleTensor, yPred: DoubleTensor): Double {
+            require(yTrue.shape.size == 1)
+            require(yTrue.shape == yPred.shape)
+
+            val diff = yTrue - yPred
+            return sqrt(diff.dot(diff)).value()
+        }
+
+        println("MSE: ${mse(alpha, alphaOLS)}")
+    }
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/PCA.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/PCA.kt
@ -0,0 +1,74 @@
+/*
+ * 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.tensors
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.tensors.core.*
+
+
+// simple PCA
+
+fun main(): Unit = Double.tensorAlgebra.withBroadcast {  // work in context with broadcast methods
+    val seed = 100500L
+
+    // assume x is range from 0 until 10
+    val x = fromArray(
+        ShapeND(10),
+        DoubleArray(10) { it.toDouble() }
+    )
+
+    // take y dependent on x with noise
+    val y = 2.0 * x + (3.0 + randomNormalLike(x, seed) * 1.5)
+
+    println("x:\n$x")
+    println("y:\n$y")
+
+    // stack them into single dataset
+    val dataset = stack(listOf(x, y)).transposed()
+
+    // normalize both x and y
+    val xMean = mean(x)
+    val yMean = mean(y)
+
+    val xStd = std(x)
+    val yStd = std(y)
+
+    val xScaled: DoubleTensor = (x - xMean) / xStd
+    val yScaled: DoubleTensor = (y - yMean) / yStd
+
+    // save means ans standard deviations for further recovery
+    val mean = fromArray(
+        ShapeND(2),
+        doubleArrayOf(xMean, yMean)
+    )
+    println("Means:\n$mean")
+
+    val std = fromArray(
+        ShapeND(2),
+        doubleArrayOf(xStd, yStd)
+    )
+    println("Standard deviations:\n$std")
+
+    // calculate the covariance matrix of scaled x and y
+    val covMatrix = covariance(listOf(xScaled.asDoubleTensor1D(), yScaled.asDoubleTensor1D()))
+    println("Covariance matrix:\n$covMatrix")
+
+    // and find out eigenvector of it
+    val (_, evecs) = symEig(covMatrix)
+    val v = evecs.getTensor(0)
+    println("Eigenvector:\n$v")
+
+    // reduce dimension of dataset
+    val datasetReduced = v dot stack(listOf(xScaled, yScaled))
+    println("Reduced data:\n$datasetReduced")
+
+    // we can restore original data from reduced data;
+    // for example, find 7th element of dataset.
+    val n = 7
+    val restored = (datasetReduced.getTensor(n) dot v.view(ShapeND(1, 2))) * std + mean
+    println("Original value:\n${dataset.getTensor(n)}")
+    println("Restored value:\n$restored")
+}
--- a/Show More
+++ b/Show More