From 00f6c027b41fae02e24e55055d821a40ce21e093 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 21 Jun 2019 12:34:04 +0300 Subject: [PATCH 01/22] Experimental type-safe dimensions --- build.gradle.kts | 2 +- buildSrc/build.gradle.kts | 2 +- kmath-dimensions/build.gradle.kts | 14 ++ .../scientifik/kmath/dimensions/Dimensions.kt | 34 +++++ .../scientifik/kmath/dimensions/Wrappers.kt | 130 ++++++++++++++++++ .../kmath/dimensions/DMatrixContextTest.kt | 27 ++++ settings.gradle.kts | 1 + 7 files changed, 208 insertions(+), 2 deletions(-) create mode 100644 kmath-dimensions/build.gradle.kts create mode 100644 kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt create mode 100644 kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt create mode 100644 kmath-dimensions/src/commonTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt diff --git a/build.gradle.kts b/build.gradle.kts index d7b7ac78a..33fe7c33a 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,4 +1,4 @@ -val kmathVersion by extra("0.1.3-dev-1") +val kmathVersion by extra("0.1.3-dev-2") allprojects { repositories { diff --git a/buildSrc/build.gradle.kts b/buildSrc/build.gradle.kts index 318189162..a2760f609 100644 --- a/buildSrc/build.gradle.kts +++ b/buildSrc/build.gradle.kts @@ -7,7 +7,7 @@ repositories { jcenter() } -val kotlinVersion = "1.3.31" +val kotlinVersion = "1.3.40" // Add plugins used in buildSrc as dependencies, also we should specify version only here dependencies { diff --git a/kmath-dimensions/build.gradle.kts b/kmath-dimensions/build.gradle.kts new file mode 100644 index 000000000..59a4c0fc3 --- /dev/null +++ b/kmath-dimensions/build.gradle.kts @@ -0,0 +1,14 @@ +plugins { + `npm-multiplatform` +} + +description = "A proof of concept module for adding typ-safe dimensions to structures" + +kotlin.sourceSets { + commonMain { + dependencies { + implementation(kotlin("reflect")) + api(project(":kmath-core")) + } + } +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt new file mode 100644 index 000000000..421ccb853 --- /dev/null +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt @@ -0,0 +1,34 @@ +package scientifik.kmath.dimensions + +interface Dimension { + val dim: UInt + + companion object { + inline fun of(dim: UInt): Dimension { + return when (dim) { + 1u -> D1 + 2u -> D2 + 3u -> D3 + else -> object : Dimension { + override val dim: UInt = dim + } + } + } + + inline fun dim(): UInt{ + return D::class.objectInstance!!.dim + } + } +} + +object D1 : Dimension { + override val dim: UInt = 1u +} + +object D2 : Dimension { + override val dim: UInt = 2u +} + +object D3 : Dimension { + override val dim: UInt = 3u +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt new file mode 100644 index 000000000..3138558b8 --- /dev/null +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt @@ -0,0 +1,130 @@ +package scientifik.kmath.dimensions + +import scientifik.kmath.linear.GenericMatrixContext +import scientifik.kmath.linear.MatrixContext +import scientifik.kmath.linear.Point +import scientifik.kmath.linear.transpose +import scientifik.kmath.operations.Ring +import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.Structure2D + +interface DMatrix : Structure2D { + companion object { + /** + * Coerces a regular matrix to a matrix with type-safe dimensions and throws a error if coercion failed + */ + inline fun coerce(structure: Structure2D): DMatrix { + if (structure.rowNum != Dimension.dim().toInt()) error("Row number mismatch: expected ${Dimension.dim()} but found ${structure.rowNum}") + if (structure.colNum != Dimension.dim().toInt()) error("Column number mismatch: expected ${Dimension.dim()} but found ${structure.colNum}") + return DMatrixWrapper(structure) + } + + /** + * The same as [coerce] but without dimension checks. Use with caution + */ + inline fun coerceUnsafe(structure: Structure2D): DMatrix { + return DMatrixWrapper(structure) + } + } +} + +inline class DMatrixWrapper( + val structure: Structure2D +) : DMatrix { + override val shape: IntArray get() = structure.shape + override fun get(i: Int, j: Int): T = structure[i, j] +} + +interface DPoint : Point { + companion object { + inline fun coerce(point: Point): DPoint { + if (point.size != Dimension.dim().toInt()) error("Vector dimension mismatch: expected ${Dimension.dim()}, but found ${point.size}") + return DPointWrapper(point) + } + + inline fun coerceUnsafe(point: Point): DPoint { + return DPointWrapper(point) + } + } +} + +inline class DPointWrapper(val point: Point) : DPoint { + override val size: Int get() = point.size + + override fun get(index: Int): T = point[index] + + override fun iterator(): Iterator = point.iterator() +} + + +/** + * Basic operations on matrices. Operates on [Matrix] + */ +inline class DMatrixContext>(val context: GenericMatrixContext) { + + inline fun Matrix.coerce(): DMatrix = + DMatrix.coerceUnsafe(this) + + /** + * Produce a matrix with this context and given dimensions + */ + inline fun produce(noinline initializer: (i: Int, j: Int) -> T): DMatrix { + val rows = Dimension.dim() + val cols = Dimension.dim() + return context.produce(rows.toInt(), cols.toInt(), initializer).coerce() + } + + inline fun point(noinline initializer: (Int) -> T): DPoint { + val size = Dimension.dim() + return DPoint.coerceUnsafe(context.point(size.toInt(), initializer)) + } + + inline infix fun DMatrix.dot( + other: DMatrix + ): DMatrix { + return context.run { this@dot dot other }.coerce() + } + + inline infix fun DMatrix.dot(vector: DPoint): DPoint { + return DPoint.coerceUnsafe(context.run { this@dot dot vector }) + } + + inline operator fun DMatrix.times(value: T): DMatrix { + return context.run { this@times.times(value) }.coerce() + } + + inline operator fun T.times(m: DMatrix): DMatrix = + m * this + + + inline operator fun DMatrix.plus(other: DMatrix): DMatrix { + return context.run { this@plus + other }.coerce() + } + + inline operator fun DMatrix.minus(other: DMatrix): DMatrix { + return context.run { this@minus + other }.coerce() + } + + inline operator fun DMatrix.unaryMinus(): DMatrix { + return context.run { this@unaryMinus.unaryMinus() }.coerce() + } + + inline fun DMatrix.transpose(): DMatrix { + return context.run { (this@transpose as Matrix).transpose() }.coerce() + } + + /** + * A square unit matrix + */ + inline fun one(): DMatrix = produce { i, j -> + if (i == j) context.elementContext.one else context.elementContext.zero + } + + inline fun zero(): DMatrix = produce { _, _ -> + context.elementContext.zero + } + + companion object { + val real = DMatrixContext(MatrixContext.real) + } +} \ No newline at end of file diff --git a/kmath-dimensions/src/commonTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt b/kmath-dimensions/src/commonTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt new file mode 100644 index 000000000..8d4bae810 --- /dev/null +++ b/kmath-dimensions/src/commonTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt @@ -0,0 +1,27 @@ +package scientifik.kmath.dimensions + +import kotlin.test.Test + + +class DMatrixContextTest { + @Test + fun testDimensionSafeMatrix() { + val res = DMatrixContext.real.run { + val m = produce { i, j -> (i + j).toDouble() } + + //The dimension of `one()` is inferred from type + (m + one()) + } + } + + @Test + fun testTypeCheck() { + val res = DMatrixContext.real.run { + val m1 = produce { i, j -> (i + j).toDouble() } + val m2 = produce { i, j -> (i + j).toDouble() } + + //Dimension-safe addition + m1.transpose() + m2 + } + } +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index 004b432fd..4afae465f 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -29,5 +29,6 @@ include( ":kmath-commons", ":kmath-koma", ":kmath-prob", + ":kmath-dimensions", ":examples" ) -- 2.34.1 From 22a3c6e4b9cc5ccb56972af66122ff0d86e77995 Mon Sep 17 00:00:00 2001 From: Peter Klimai Date: Sun, 14 Jul 2019 18:22:22 +0300 Subject: [PATCH 02/22] Initial implementation of kmath-for-real module based on https://github.com/thomasnield/numky --- kmath-for-real/build.gradle.kts | 9 ++ .../kmath/real/DoubleMatrixOperations.kt | 117 ++++++++++++++++++ settings.gradle.kts | 1 + 3 files changed, 127 insertions(+) create mode 100644 kmath-for-real/build.gradle.kts create mode 100644 kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt diff --git a/kmath-for-real/build.gradle.kts b/kmath-for-real/build.gradle.kts new file mode 100644 index 000000000..7eaa5e174 --- /dev/null +++ b/kmath-for-real/build.gradle.kts @@ -0,0 +1,9 @@ +plugins { + `npm-multiplatform` +} + +kotlin.sourceSets.commonMain { + dependencies { + api(project(":kmath-core")) + } +} \ No newline at end of file diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt new file mode 100644 index 000000000..54e711a7b --- /dev/null +++ b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt @@ -0,0 +1,117 @@ +package scientifik.kmath.real + +import scientifik.kmath.linear.MatrixContext +import scientifik.kmath.linear.RealMatrixContext.elementContext +import scientifik.kmath.linear.VirtualMatrix +import scientifik.kmath.operations.sum +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.asSequence +import kotlin.math.pow + +// Initial implementation of these functions is taken from: +// https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt + +fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) = MatrixContext.real.produce(rowNum, colNum, initializer) + +fun Sequence.toMatrix() = toList().let { + MatrixContext.real.produce(it.size,it[0].size) { row, col -> it[row][col] } +} + +operator fun Matrix.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@times[row, col] * double +} + +fun Matrix.square() = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@square[row,col].pow(2) +} + +operator fun Matrix.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@plus[row,col] + double +} + +operator fun Matrix.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@minus[row,col] - double +} + +operator fun Matrix.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@div[row,col] / double +} + +operator fun Double.times(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + matrix[row,col] * this +} + +operator fun Double.plus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + matrix[row,col] + this +} + +operator fun Double.minus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + matrix[row,col] - this +} + +operator fun Double.div(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> + matrix[row,col] / this +} + +operator fun Matrix.times(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@times[row,col] * other[row,col] +} + +operator fun Matrix.minus(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { row, col -> + this@minus[row,col] - other[row,col] +} + +operator fun Matrix.plus(other: Matrix) = MatrixContext.real.add(this,other) + +fun Matrix.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) { row, col -> + get(if (row == 0) 0 else row % rowNum, col) +} + +inline fun Matrix.appendColumn(crossinline mapper: (Buffer) -> Double) = + MatrixContext.real.produce(rowNum,colNum+1) { row,col -> + if (col < colNum) + this[row,col] + else + mapper(rows[row]) + } + +fun Matrix.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex) + +fun Matrix.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) { row, col -> + this@extractColumns[row, columnRange.start + col] +} + +fun Matrix.sumByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> + val column = columns[j] + with(elementContext) { + sum(column.asSequence()) + } +} + +fun Matrix.minByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> + val column = columns[j] + column.asSequence().min()?:throw Exception("Cannot produce min on empty column") +} + +fun Matrix.maxByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> + val column = columns[j] + column.asSequence().max()?:throw Exception("Cannot produce min on empty column") +} + +fun Matrix.averageByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> + val column = columns[j] + column.asSequence().average() +} + +fun Matrix.sum() = this.elements().map { (_,value) -> value }.sum() + +fun Matrix.min() = this.elements().map { (_,value) -> value }.min() + +fun Matrix.max() = this.elements().map { (_,value) -> value }.max() + +fun Matrix.average() = this.elements().map { (_,value) -> value }.average() + +fun Matrix.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) { i, j -> + this@pow[i,j].pow(n) +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index 004b432fd..e7f27714c 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -29,5 +29,6 @@ include( ":kmath-commons", ":kmath-koma", ":kmath-prob", + ":kmath-for-real", ":examples" ) -- 2.34.1 From c31c9a9532d86aa30768763c286b030230f8d61c Mon Sep 17 00:00:00 2001 From: Peter Klimai Date: Wed, 7 Aug 2019 21:47:50 +0300 Subject: [PATCH 03/22] Move RealVector to kmath-for-real and use scientifik plugin --- .../scientifik/kmath/linear/MatrixTest.kt | 30 ---------------- kmath-for-real/build.gradle.kts | 10 +++--- .../scientifik/kmath/linear/RealVector.kt | 0 .../scientifik/kmath/linear/VectorTest.kt | 36 +++++++++++++++++++ kmath-histograms/build.gradle.kts | 1 + 5 files changed, 43 insertions(+), 34 deletions(-) rename {kmath-core => kmath-for-real}/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt (100%) create mode 100644 kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt index de13c9888..d27aa665c 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt @@ -6,21 +6,6 @@ import kotlin.test.assertEquals class MatrixTest { - @Test - fun testSum() { - val vector1 = RealVector(5) { it.toDouble() } - val vector2 = RealVector(5) { 5 - it.toDouble() } - val sum = vector1 + vector2 - assertEquals(5.0, sum[2]) - } - - @Test - fun testVectorToMatrix() { - val vector = RealVector(5) { it.toDouble() } - val matrix = vector.asMatrix() - assertEquals(4.0, matrix[4, 0]) - } - @Test fun testTranspose() { val matrix = MatrixContext.real.one(3, 3) @@ -28,21 +13,6 @@ class MatrixTest { assertEquals(matrix, transposed) } - - @Test - fun testDot() { - val vector1 = RealVector(5) { it.toDouble() } - val vector2 = RealVector(5) { 5 - it.toDouble() } - - val matrix1 = vector1.asMatrix() - val matrix2 = vector2.asMatrix().transpose() - val product = MatrixContext.real.run { matrix1 dot matrix2 } - - - assertEquals(5.0, product[1, 0]) - assertEquals(6.0, product[2, 2]) - } - @Test fun testBuilder() { val matrix = Matrix.build(2, 3)( diff --git a/kmath-for-real/build.gradle.kts b/kmath-for-real/build.gradle.kts index 7eaa5e174..a8a8975bc 100644 --- a/kmath-for-real/build.gradle.kts +++ b/kmath-for-real/build.gradle.kts @@ -1,9 +1,11 @@ plugins { - `npm-multiplatform` + id("scientifik.mpp") } -kotlin.sourceSets.commonMain { - dependencies { - api(project(":kmath-core")) +kotlin.sourceSets { + commonMain { + dependencies { + api(project(":kmath-core")) + } } } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt similarity index 100% rename from kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt rename to kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt diff --git a/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt b/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt new file mode 100644 index 000000000..0e73ee4a5 --- /dev/null +++ b/kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt @@ -0,0 +1,36 @@ +package scientifik.kmath.linear + +import kotlin.test.Test +import kotlin.test.assertEquals + +class VectorTest { + @Test + fun testSum() { + val vector1 = RealVector(5) { it.toDouble() } + val vector2 = RealVector(5) { 5 - it.toDouble() } + val sum = vector1 + vector2 + assertEquals(5.0, sum[2]) + } + + @Test + fun testVectorToMatrix() { + val vector = RealVector(5) { it.toDouble() } + val matrix = vector.asMatrix() + assertEquals(4.0, matrix[4, 0]) + } + + @Test + fun testDot() { + val vector1 = RealVector(5) { it.toDouble() } + val vector2 = RealVector(5) { 5 - it.toDouble() } + + val matrix1 = vector1.asMatrix() + val matrix2 = vector2.asMatrix().transpose() + val product = MatrixContext.real.run { matrix1 dot matrix2 } + + + assertEquals(5.0, product[1, 0]) + assertEquals(6.0, product[2, 2]) + } + +} \ No newline at end of file diff --git a/kmath-histograms/build.gradle.kts b/kmath-histograms/build.gradle.kts index 39aa833ad..993bfed8e 100644 --- a/kmath-histograms/build.gradle.kts +++ b/kmath-histograms/build.gradle.kts @@ -5,5 +5,6 @@ plugins { kotlin.sourceSets.commonMain { dependencies { api(project(":kmath-core")) + api(project(":kmath-for-real")) } } \ No newline at end of file -- 2.34.1 From fbe7988bd03156f5dd59b023deb0772e041cacda Mon Sep 17 00:00:00 2001 From: Peter Klimai Date: Tue, 3 Sep 2019 20:35:50 +0300 Subject: [PATCH 04/22] Some refactoring, add test file --- .../kmath/real/DoubleMatrixOperations.kt | 137 +++++++++++------- .../scientific.kmath.real/RealMatrixTest.kt | 16 ++ kmath-for-real/src/jvmMain/kotlin/.gitkeep | 0 3 files changed, 99 insertions(+), 54 deletions(-) create mode 100644 kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt create mode 100644 kmath-for-real/src/jvmMain/kotlin/.gitkeep diff --git a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt index 54e711a7b..d9c6b5eab 100644 --- a/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt +++ b/kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt @@ -9,109 +9,138 @@ import scientifik.kmath.structures.Matrix import scientifik.kmath.structures.asSequence import kotlin.math.pow -// Initial implementation of these functions is taken from: -// https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt +/* + * Functions for convenient "numpy-like" operations with Double matrices. + * + * Initial implementation of these functions is taken from: + * https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt + * + */ -fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) = MatrixContext.real.produce(rowNum, colNum, initializer) +/* + * Functions that help create a real (Double) matrix + */ + +fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) = + MatrixContext.real.produce(rowNum, colNum, initializer) fun Sequence.toMatrix() = toList().let { - MatrixContext.real.produce(it.size,it[0].size) { row, col -> it[row][col] } + MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] } } -operator fun Matrix.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@times[row, col] * double +fun Matrix.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) { + row, col -> get(if (row == 0) 0 else row % rowNum, col) } -fun Matrix.square() = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@square[row,col].pow(2) +/* + * Operations for matrix and real number + */ + +operator fun Matrix.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] * double } -operator fun Matrix.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@plus[row,col] + double +operator fun Matrix.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] + double } -operator fun Matrix.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@minus[row,col] - double +operator fun Matrix.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] - double } -operator fun Matrix.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@div[row,col] / double +operator fun Matrix.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] / double } -operator fun Double.times(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> - matrix[row,col] * this +operator fun Double.times(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { + row, col -> this * matrix[row, col] } -operator fun Double.plus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> - matrix[row,col] + this +operator fun Double.plus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { + row, col -> this * matrix[row, col] } -operator fun Double.minus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> - matrix[row,col] - this +operator fun Double.minus(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { + row, col -> this - matrix[row, col] } -operator fun Double.div(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> - matrix[row,col] / this +// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'? +//operator fun Double.div(matrix: Matrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { +// row, col -> matrix[row, col] / this +//} + +/* + * Per-element (!) square and power operations + */ + +fun Matrix.square() = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col].pow(2) } -operator fun Matrix.times(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@times[row,col] * other[row,col] +fun Matrix.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) { + i, j -> this[i, j].pow(n) } -operator fun Matrix.minus(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this@minus[row,col] - other[row,col] +/* + * Operations on two matrices (per-element!) + */ + +operator fun Matrix.times(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row, col] * other[row, col] } -operator fun Matrix.plus(other: Matrix) = MatrixContext.real.add(this,other) +operator fun Matrix.plus(other: Matrix) = MatrixContext.real.add(this, other) -fun Matrix.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) { row, col -> - get(if (row == 0) 0 else row % rowNum, col) +operator fun Matrix.minus(other: Matrix) = MatrixContext.real.produce(rowNum, colNum) { + row, col -> this[row,col] - other[row,col] } -inline fun Matrix.appendColumn(crossinline mapper: (Buffer) -> Double) = - MatrixContext.real.produce(rowNum,colNum+1) { row,col -> +/* + * Operations on columns + */ + +inline fun Matrix.appendColumn(crossinline mapper: (Buffer) -> Double) = + MatrixContext.real.produce(rowNum,colNum+1) { + row, col -> if (col < colNum) - this[row,col] + this[row, col] else mapper(rows[row]) - } + } + +fun Matrix.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) { + row, col -> this[row, columnRange.first + col] +} fun Matrix.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex) -fun Matrix.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) { row, col -> - this@extractColumns[row, columnRange.start + col] -} - -fun Matrix.sumByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> +fun Matrix.sumByColumn() = MatrixContext.real.produce(1, colNum) { _, j -> val column = columns[j] with(elementContext) { sum(column.asSequence()) } } -fun Matrix.minByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> - val column = columns[j] - column.asSequence().min()?:throw Exception("Cannot produce min on empty column") +fun Matrix.minByColumn() = MatrixContext.real.produce(1, colNum) { + _, j -> columns[j].asSequence().min() ?: throw Exception("Cannot produce min on empty column") } -fun Matrix.maxByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> - val column = columns[j] - column.asSequence().max()?:throw Exception("Cannot produce min on empty column") +fun Matrix.maxByColumn() = MatrixContext.real.produce(1, colNum) { + _, j -> columns[j].asSequence().max() ?: throw Exception("Cannot produce min on empty column") } -fun Matrix.averageByColumn() = MatrixContext.real.produce(1, colNum) { i, j -> - val column = columns[j] - column.asSequence().average() +fun Matrix.averageByColumn() = MatrixContext.real.produce(1, colNum) { + _, j -> columns[j].asSequence().average() } -fun Matrix.sum() = this.elements().map { (_,value) -> value }.sum() +/* + * Operations processing all elements + */ -fun Matrix.min() = this.elements().map { (_,value) -> value }.min() +fun Matrix.sum() = elements().map { (_, value) -> value }.sum() -fun Matrix.max() = this.elements().map { (_,value) -> value }.max() +fun Matrix.min() = elements().map { (_, value) -> value }.min() -fun Matrix.average() = this.elements().map { (_,value) -> value }.average() +fun Matrix.max() = elements().map { (_, value) -> value }.max() -fun Matrix.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) { i, j -> - this@pow[i,j].pow(n) -} \ No newline at end of file +fun Matrix.average() = elements().map { (_, value) -> value }.average() diff --git a/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt b/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt new file mode 100644 index 000000000..808794442 --- /dev/null +++ b/kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt @@ -0,0 +1,16 @@ +package scientific.kmath.real + +import scientifik.kmath.real.average +import scientifik.kmath.real.realMatrix +import scientifik.kmath.real.sum +import kotlin.test.Test +import kotlin.test.assertEquals + +class RealMatrixTest { + @Test + fun testSum() { + val m = realMatrix(10, 10) { i, j -> (i + j).toDouble() } + assertEquals(m.sum(), 900.0) + assertEquals(m.average(), 9.0) + } +} diff --git a/kmath-for-real/src/jvmMain/kotlin/.gitkeep b/kmath-for-real/src/jvmMain/kotlin/.gitkeep new file mode 100644 index 000000000..e69de29bb -- 2.34.1 From f9836a64779c61a510716641b9f45c5ec658d80f Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 8 Dec 2019 15:48:25 +0300 Subject: [PATCH 05/22] Viktor prototype --- examples/build.gradle.kts | 3 +- .../scientifik/kmath/structures/NDField.kt | 30 +++----- .../scientifik/kmath/structures/VictorTest.kt | 42 +++++++++++ kmath-commons/build.gradle.kts | 2 - .../kmath/structures/ExtendedNDField.kt | 1 - kmath-viktor/build.gradle.kts | 10 +++ .../scientifik/kmath/viktor/ViktorBuffer.kt | 20 +++++ .../kmath/viktor/ViktorNDStructure.kt | 73 +++++++++++++++++++ settings.gradle.kts | 1 + 9 files changed, 160 insertions(+), 22 deletions(-) create mode 100644 examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt create mode 100644 kmath-viktor/build.gradle.kts create mode 100644 kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt create mode 100644 kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index ad59afc5c..dfd530618 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -4,7 +4,7 @@ import org.jetbrains.kotlin.gradle.tasks.KotlinCompile plugins { java kotlin("jvm") - kotlin("plugin.allopen") version "1.3.60" + kotlin("plugin.allopen") version "1.3.61" id("kotlinx.benchmark") version "0.2.0-dev-5" } @@ -29,6 +29,7 @@ dependencies { implementation(project(":kmath-coroutines")) implementation(project(":kmath-commons")) implementation(project(":kmath-koma")) + implementation(project(":kmath-viktor")) implementation("com.kyonifer:koma-core-ejml:0.12") implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:${Scientifik.ioVersion}") diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt b/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt index 7e8c433c0..cfd1206ff 100644 --- a/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt +++ b/examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt @@ -4,7 +4,13 @@ import kotlinx.coroutines.GlobalScope import scientifik.kmath.operations.RealField import kotlin.system.measureTimeMillis -fun main(args: Array) { +internal inline fun measureAndPrint(title: String, block: () -> Unit) { + val time = measureTimeMillis(block) + println("$title completed in $time millis") +} + + +fun main() { val dim = 1000 val n = 1000 @@ -15,8 +21,7 @@ fun main(args: Array) { //A generic boxing field. It should be used for objects, not primitives. val genericField = NDField.boxing(RealField, dim, dim) - - val autoTime = measureTimeMillis { + measureAndPrint("Automatic field addition") { autoField.run { var res = one repeat(n) { @@ -25,18 +30,14 @@ fun main(args: Array) { } } - println("Automatic field addition completed in $autoTime millis") - - val elementTime = measureTimeMillis { + measureAndPrint("Element addition"){ var res = genericField.one repeat(n) { res += 1.0 } } - println("Element addition completed in $elementTime millis") - - val specializedTime = measureTimeMillis { + measureAndPrint("Specialized addition") { specializedField.run { var res: NDBuffer = one repeat(n) { @@ -45,10 +46,7 @@ fun main(args: Array) { } } - println("Specialized addition completed in $specializedTime millis") - - - val lazyTime = measureTimeMillis { + measureAndPrint("Lazy addition") { val res = specializedField.one.mapAsync(GlobalScope) { var c = 0.0 repeat(n) { @@ -60,9 +58,7 @@ fun main(args: Array) { res.elements().forEach { it.second } } - println("Lazy addition completed in $lazyTime millis") - - val genericTime = measureTimeMillis { + measureAndPrint("Generic addition") { //genericField.run(action) genericField.run { var res: NDBuffer = one @@ -72,6 +68,4 @@ fun main(args: Array) { } } - println("Generic addition completed in $genericTime millis") - } \ No newline at end of file diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt b/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt new file mode 100644 index 000000000..da0c193ad --- /dev/null +++ b/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt @@ -0,0 +1,42 @@ +package scientifik.kmath.structures + +import org.jetbrains.bio.viktor.F64Array +import scientifik.kmath.operations.RealField +import scientifik.kmath.viktor.ViktorNDField + +fun main() { + val dim = 1000 + val n = 400 + + // automatically build context most suited for given type. + val autoField = NDField.auto(RealField, dim, dim) + + val viktorField = ViktorNDField(intArrayOf(dim, dim)) + + + measureAndPrint("Automatic field addition") { + autoField.run { + var res = one + repeat(n) { + res += 1.0 + } + } + } + + measureAndPrint("Raw Viktor") { + val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim)) + var res = one + repeat(n) { + res = res + one + } + } + + measureAndPrint("Viktor field addition") { + viktorField.run { + var res = one + repeat(n) { + res += 1.0 + } + } + } +} \ No newline at end of file diff --git a/kmath-commons/build.gradle.kts b/kmath-commons/build.gradle.kts index 9b446f872..9ac0f6e95 100644 --- a/kmath-commons/build.gradle.kts +++ b/kmath-commons/build.gradle.kts @@ -9,6 +9,4 @@ dependencies { api(project(":kmath-coroutines")) api(project(":kmath-prob")) api("org.apache.commons:commons-math3:3.6.1") - testImplementation("org.jetbrains.kotlin:kotlin-test") - testImplementation("org.jetbrains.kotlin:kotlin-test-junit") } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt index 370dc646e..3437644ff 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt @@ -2,7 +2,6 @@ package scientifik.kmath.structures import scientifik.kmath.operations.* - interface ExtendedNDField> : NDField, TrigonometricOperations, diff --git a/kmath-viktor/build.gradle.kts b/kmath-viktor/build.gradle.kts new file mode 100644 index 000000000..52ee7c497 --- /dev/null +++ b/kmath-viktor/build.gradle.kts @@ -0,0 +1,10 @@ +plugins { + id("scientifik.jvm") +} + +description = "Binding for https://github.com/JetBrains-Research/viktor" + +dependencies { + api(project(":kmath-core")) + api("org.jetbrains.bio:viktor:1.0.1") +} \ No newline at end of file diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt new file mode 100644 index 000000000..040eee951 --- /dev/null +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt @@ -0,0 +1,20 @@ +package scientifik.kmath.viktor + +import org.jetbrains.bio.viktor.F64FlatArray +import scientifik.kmath.structures.MutableBuffer + +@Suppress("NOTHING_TO_INLINE", "OVERRIDE_BY_INLINE") +inline class ViktorBuffer(val flatArray: F64FlatArray) : MutableBuffer { + override val size: Int get() = flatArray.size + + override inline fun get(index: Int): Double = flatArray[index] + override inline fun set(index: Int, value: Double) { + flatArray[index] = value + } + + override fun copy(): MutableBuffer { + return ViktorBuffer(flatArray.copy().flatten()) + } + + override fun iterator(): Iterator = flatArray.data.iterator() +} \ No newline at end of file diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt new file mode 100644 index 000000000..b572485ff --- /dev/null +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt @@ -0,0 +1,73 @@ +package scientifik.kmath.viktor + +import org.jetbrains.bio.viktor.F64Array +import scientifik.kmath.operations.RealField +import scientifik.kmath.structures.DefaultStrides +import scientifik.kmath.structures.MutableNDStructure +import scientifik.kmath.structures.NDField + +inline class ViktorNDStructure(val f64Buffer: F64Array) : MutableNDStructure { + + override val shape: IntArray get() = f64Buffer.shape + + override fun get(index: IntArray): Double = f64Buffer.get(*index) + + override fun set(index: IntArray, value: Double) { + f64Buffer.set(value = value, indices = *index) + } + + override fun elements(): Sequence> { + return DefaultStrides(shape).indices().map { it to get(it) } + } +} + +fun F64Array.asStructure(): ViktorNDStructure = ViktorNDStructure(this) + +class ViktorNDField(override val shape: IntArray) : NDField { + override val zero: ViktorNDStructure + get() = F64Array.full(init = 0.0, shape = *shape).asStructure() + override val one: ViktorNDStructure + get() = F64Array.full(init = 1.0, shape = *shape).asStructure() + + val strides = DefaultStrides(shape) + + override val elementContext: RealField get() = RealField + + override fun produce(initializer: RealField.(IntArray) -> Double): ViktorNDStructure = F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.initializer(index), indices = *index) + } + }.asStructure() + + override fun map(arg: ViktorNDStructure, transform: RealField.(Double) -> Double): ViktorNDStructure = + F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.transform(arg[index]), indices = *index) + } + }.asStructure() + + override fun mapIndexed( + arg: ViktorNDStructure, + transform: RealField.(index: IntArray, Double) -> Double + ): ViktorNDStructure = F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.transform(index, arg[index]), indices = *index) + } + }.asStructure() + + override fun combine( + a: ViktorNDStructure, + b: ViktorNDStructure, + transform: RealField.(Double, Double) -> Double + ): ViktorNDStructure = F64Array(*shape).apply { + this@ViktorNDField.strides.indices().forEach { index -> + set(value = RealField.transform(a[index], b[index]), indices = *index) + } + }.asStructure() + + override fun ViktorNDStructure.plus(b: ViktorNDStructure): ViktorNDStructure = + (f64Buffer + b.f64Buffer).asStructure() + + override fun ViktorNDStructure.minus(b: ViktorNDStructure): ViktorNDStructure = + (f64Buffer - b.f64Buffer).asStructure() +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index f52edfac5..8c0c4420b 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -33,6 +33,7 @@ include( ":kmath-coroutines", ":kmath-histograms", ":kmath-commons", + ":kmath-viktor", ":kmath-koma", ":kmath-prob", ":kmath-io", -- 2.34.1 From cbac4fca334bb81a434798d325747be3b03bb00f Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 8 Dec 2019 17:25:36 +0300 Subject: [PATCH 06/22] Viktor prototype with compiler bug --- .../scientifik/kmath/structures/VictorTest.kt | 31 +++++++++++++------ .../kmath/viktor/ViktorNDStructure.kt | 14 ++++++--- 2 files changed, 32 insertions(+), 13 deletions(-) diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt b/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt index da0c193ad..8d710c584 100644 --- a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt +++ b/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt @@ -13,6 +13,12 @@ fun main() { val viktorField = ViktorNDField(intArrayOf(dim, dim)) + autoField.run { + var res = one + repeat(n/2) { + res += 1.0 + } + } measureAndPrint("Automatic field addition") { autoField.run { @@ -23,6 +29,22 @@ fun main() { } } + viktorField.run { + var res = one + repeat(n/2) { + res += one + } + } + + measureAndPrint("Viktor field addition") { + viktorField.run { + var res = one + repeat(n) { + res += one + } + } + } + measureAndPrint("Raw Viktor") { val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim)) var res = one @@ -30,13 +52,4 @@ fun main() { res = res + one } } - - measureAndPrint("Viktor field addition") { - viktorField.run { - var res = one - repeat(n) { - res += 1.0 - } - } - } } \ No newline at end of file diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt index b572485ff..b8cf688f4 100644 --- a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt @@ -65,9 +65,15 @@ class ViktorNDField(override val shape: IntArray) : NDField Date: Sun, 8 Dec 2019 19:52:47 +0300 Subject: [PATCH 07/22] Viktor prototype with inline --- .../scientifik/kmath/viktor/ViktorNDStructure.kt | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt index b8cf688f4..9654ba283 100644 --- a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt @@ -69,9 +69,14 @@ class ViktorNDField(override val shape: IntArray) : NDField Date: Mon, 9 Dec 2019 17:37:17 +0300 Subject: [PATCH 08/22] Viktor prototype with inline --- .../scientifik/kmath/structures/VictorTest.kt | 21 +++++++++++++++ .../kmath/viktor/ViktorNDStructure.kt | 26 ++++++++++--------- 2 files changed, 35 insertions(+), 12 deletions(-) diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt b/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt index 8d710c584..19469ca72 100644 --- a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt +++ b/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt @@ -4,12 +4,14 @@ import org.jetbrains.bio.viktor.F64Array import scientifik.kmath.operations.RealField import scientifik.kmath.viktor.ViktorNDField + fun main() { val dim = 1000 val n = 400 // automatically build context most suited for given type. val autoField = NDField.auto(RealField, dim, dim) + val realField = NDField.real(dim,dim) val viktorField = ViktorNDField(intArrayOf(dim, dim)) @@ -52,4 +54,23 @@ fun main() { res = res + one } } + + measureAndPrint("Automatic field log") { + realField.run { + val fortyTwo = produce { 42.0 } + var res = one + + repeat(n) { + res = ln(fortyTwo) + } + } + } + + measureAndPrint("Raw Viktor log") { + val fortyTwo = F64Array.full(dim, dim, init = 42.0) + var res: F64Array + repeat(n) { + res = fortyTwo.log() + } + } } \ No newline at end of file diff --git a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt index 9654ba283..84e927721 100644 --- a/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt +++ b/kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt @@ -6,14 +6,15 @@ import scientifik.kmath.structures.DefaultStrides import scientifik.kmath.structures.MutableNDStructure import scientifik.kmath.structures.NDField +@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") inline class ViktorNDStructure(val f64Buffer: F64Array) : MutableNDStructure { override val shape: IntArray get() = f64Buffer.shape - override fun get(index: IntArray): Double = f64Buffer.get(*index) + override inline fun get(index: IntArray): Double = f64Buffer.get(*index) - override fun set(index: IntArray, value: Double) { - f64Buffer.set(value = value, indices = *index) + override inline fun set(index: IntArray, value: Double) { + f64Buffer.set(*index, value = value) } override fun elements(): Sequence> { @@ -23,6 +24,7 @@ inline class ViktorNDStructure(val f64Buffer: F64Array) : MutableNDStructure { override val zero: ViktorNDStructure get() = F64Array.full(init = 0.0, shape = *shape).asStructure() @@ -69,16 +71,16 @@ class ViktorNDField(override val shape: IntArray) : NDField Date: Mon, 9 Dec 2019 19:52:00 +0300 Subject: [PATCH 09/22] Examples for type-safe dimensions --- examples/build.gradle.kts | 11 ++--- .../kmath/structures/ViktorBenchmark.kt} | 47 +++++++++---------- .../kotlin/scientifik/kmath/utils/utils.kt | 8 ++++ .../kmath/structures/typeSafeDimensions.kt | 34 ++++++++++++++ kmath-dimensions/build.gradle.kts | 2 +- .../scientifik/kmath/dimensions/Dimensions.kt | 9 +++- .../scientifik/kmath/dimensions/Wrappers.kt | 26 ++++++++-- 7 files changed, 96 insertions(+), 41 deletions(-) rename examples/src/{main/kotlin/scientifik/kmath/structures/VictorTest.kt => benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt} (58%) create mode 100644 examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt create mode 100644 examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index dfd530618..47acaa5ba 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -5,7 +5,7 @@ plugins { java kotlin("jvm") kotlin("plugin.allopen") version "1.3.61" - id("kotlinx.benchmark") version "0.2.0-dev-5" + id("kotlinx.benchmark") version "0.2.0-dev-7" } configure { @@ -13,10 +13,7 @@ configure { } repositories { - maven("https://dl.bintray.com/kotlin/kotlin-eap") maven("http://dl.bintray.com/kyonifer/maven") - maven ("https://dl.bintray.com/orangy/maven") - mavenCentral() } @@ -30,12 +27,10 @@ dependencies { implementation(project(":kmath-commons")) implementation(project(":kmath-koma")) implementation(project(":kmath-viktor")) + implementation(project(":kmath-dimensions")) implementation("com.kyonifer:koma-core-ejml:0.12") implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:${Scientifik.ioVersion}") - - implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-2") - - + implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-7") "benchmarksCompile"(sourceSets.main.get().compileClasspath) } diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt b/examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt similarity index 58% rename from examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt rename to examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt index 19469ca72..be4115d81 100644 --- a/examples/src/main/kotlin/scientifik/kmath/structures/VictorTest.kt +++ b/examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt @@ -1,28 +1,26 @@ package scientifik.kmath.structures import org.jetbrains.bio.viktor.F64Array +import org.openjdk.jmh.annotations.Benchmark +import org.openjdk.jmh.annotations.Scope +import org.openjdk.jmh.annotations.State import scientifik.kmath.operations.RealField import scientifik.kmath.viktor.ViktorNDField -fun main() { - val dim = 1000 - val n = 400 +@State(Scope.Benchmark) +class ViktorBenchmark { + final val dim = 1000 + final val n = 100 // automatically build context most suited for given type. - val autoField = NDField.auto(RealField, dim, dim) - val realField = NDField.real(dim,dim) + final val autoField = NDField.auto(RealField, dim, dim) + final val realField = NDField.real(dim, dim) - val viktorField = ViktorNDField(intArrayOf(dim, dim)) + final val viktorField = ViktorNDField(intArrayOf(dim, dim)) - autoField.run { - var res = one - repeat(n/2) { - res += 1.0 - } - } - - measureAndPrint("Automatic field addition") { + @Benchmark + fun `Automatic field addition`() { autoField.run { var res = one repeat(n) { @@ -31,14 +29,8 @@ fun main() { } } - viktorField.run { - var res = one - repeat(n/2) { - res += one - } - } - - measureAndPrint("Viktor field addition") { + @Benchmark + fun `Viktor field addition`() { viktorField.run { var res = one repeat(n) { @@ -47,7 +39,8 @@ fun main() { } } - measureAndPrint("Raw Viktor") { + @Benchmark + fun `Raw Viktor`() { val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim)) var res = one repeat(n) { @@ -55,9 +48,10 @@ fun main() { } } - measureAndPrint("Automatic field log") { + @Benchmark + fun `Real field log`() { realField.run { - val fortyTwo = produce { 42.0 } + val fortyTwo = produce { 42.0 } var res = one repeat(n) { @@ -66,7 +60,8 @@ fun main() { } } - measureAndPrint("Raw Viktor log") { + @Benchmark + fun `Raw Viktor log`() { val fortyTwo = F64Array.full(dim, dim, init = 42.0) var res: F64Array repeat(n) { diff --git a/examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt b/examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt new file mode 100644 index 000000000..6ec9e9c17 --- /dev/null +++ b/examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt @@ -0,0 +1,8 @@ +package scientifik.kmath.utils + +import kotlin.system.measureTimeMillis + +internal inline fun measureAndPrint(title: String, block: () -> Unit) { + val time = measureTimeMillis(block) + println("$title completed in $time millis") +} \ No newline at end of file diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt b/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt new file mode 100644 index 000000000..8693308ca --- /dev/null +++ b/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt @@ -0,0 +1,34 @@ +package scientifik.kmath.structures + +import scientifik.kmath.dimensions.D2 +import scientifik.kmath.dimensions.D3 +import scientifik.kmath.dimensions.DMatrixContext +import scientifik.kmath.dimensions.Dimension +import scientifik.kmath.operations.RealField + +fun DMatrixContext.simple() { + val m1 = produce { i, j -> (i + j).toDouble() } + val m2 = produce { i, j -> (i + j).toDouble() } + + //Dimension-safe addition + m1.transpose() + m2 +} + +object D5: Dimension{ + override val dim: UInt = 5u +} + +fun DMatrixContext.custom() { + val m1 = produce { i, j -> (i+j).toDouble() } + val m2 = produce { i, j -> (i-j).toDouble() } + val m3 = produce { i, j -> (i-j).toDouble() } + + (m1 dot m2) + m3 +} + +fun main() { + DMatrixContext.real.run { + simple() + custom() + } +} \ No newline at end of file diff --git a/kmath-dimensions/build.gradle.kts b/kmath-dimensions/build.gradle.kts index 59a4c0fc3..50fb41391 100644 --- a/kmath-dimensions/build.gradle.kts +++ b/kmath-dimensions/build.gradle.kts @@ -1,5 +1,5 @@ plugins { - `npm-multiplatform` + id("scientifik.mpp") } description = "A proof of concept module for adding typ-safe dimensions to structures" diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt index 421ccb853..87511b885 100644 --- a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt @@ -1,9 +1,14 @@ package scientifik.kmath.dimensions +/** + * An interface which is not used in runtime. Designates a size of some structure. + * All descendants should be singleton objects. + */ interface Dimension { val dim: UInt companion object { + @Suppress("NOTHING_TO_INLINE") inline fun of(dim: UInt): Dimension { return when (dim) { 1u -> D1 @@ -15,8 +20,8 @@ interface Dimension { } } - inline fun dim(): UInt{ - return D::class.objectInstance!!.dim + inline fun dim(): UInt { + return D::class.objectInstance?.dim ?: error("Dimension object must be a singleton") } } } diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt index 3138558b8..9985888b9 100644 --- a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt @@ -8,14 +8,21 @@ import scientifik.kmath.operations.Ring import scientifik.kmath.structures.Matrix import scientifik.kmath.structures.Structure2D +/** + * A matrix with compile-time controlled dimension + */ interface DMatrix : Structure2D { companion object { /** * Coerces a regular matrix to a matrix with type-safe dimensions and throws a error if coercion failed */ inline fun coerce(structure: Structure2D): DMatrix { - if (structure.rowNum != Dimension.dim().toInt()) error("Row number mismatch: expected ${Dimension.dim()} but found ${structure.rowNum}") - if (structure.colNum != Dimension.dim().toInt()) error("Column number mismatch: expected ${Dimension.dim()} but found ${structure.colNum}") + if (structure.rowNum != Dimension.dim().toInt()) { + error("Row number mismatch: expected ${Dimension.dim()} but found ${structure.rowNum}") + } + if (structure.colNum != Dimension.dim().toInt()) { + error("Column number mismatch: expected ${Dimension.dim()} but found ${structure.colNum}") + } return DMatrixWrapper(structure) } @@ -28,6 +35,9 @@ interface DMatrix : Structure2D { } } +/** + * An inline wrapper for a Matrix + */ inline class DMatrixWrapper( val structure: Structure2D ) : DMatrix { @@ -35,10 +45,15 @@ inline class DMatrixWrapper( override fun get(i: Int, j: Int): T = structure[i, j] } +/** + * Dimension-safe point + */ interface DPoint : Point { companion object { inline fun coerce(point: Point): DPoint { - if (point.size != Dimension.dim().toInt()) error("Vector dimension mismatch: expected ${Dimension.dim()}, but found ${point.size}") + if (point.size != Dimension.dim().toInt()) { + error("Vector dimension mismatch: expected ${Dimension.dim()}, but found ${point.size}") + } return DPointWrapper(point) } @@ -48,6 +63,9 @@ interface DPoint : Point { } } +/** + * Dimension-safe point wrapper + */ inline class DPointWrapper(val point: Point) : DPoint { override val size: Int get() = point.size @@ -58,7 +76,7 @@ inline class DPointWrapper(val point: Point) : DPoint /** - * Basic operations on matrices. Operates on [Matrix] + * Basic operations on dimension-safe matrices. Operates on [Matrix] */ inline class DMatrixContext>(val context: GenericMatrixContext) { -- 2.34.1 From 09f2f377802f24cb522c741cf7eddf0be8ba1d0b Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Mon, 9 Dec 2019 19:52:48 +0300 Subject: [PATCH 10/22] Examples for type-safe dimensions --- buildSrc/build.gradle.kts | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 buildSrc/build.gradle.kts diff --git a/buildSrc/build.gradle.kts b/buildSrc/build.gradle.kts deleted file mode 100644 index e69de29bb..000000000 -- 2.34.1 From 169801060b6c1adc8fee15820b7350a693f2e36e Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Mon, 9 Dec 2019 19:57:26 +0300 Subject: [PATCH 11/22] Examples for type-safe dimensions --- README.md | 2 ++ .../scientifik/kmath/structures/typeSafeDimensions.kt | 8 ++++---- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 5678fc530..0df279889 100644 --- a/README.md +++ b/README.md @@ -40,6 +40,8 @@ can be used for a wide variety of purposes from high performance calculations to * **Streaming** Streaming operations on mathematical objects and objects buffers. +* **Type-safe dimensions** Type-safe dimensions for matrix operations. + * **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free to submit a feature request if you want something to be done first. diff --git a/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt b/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt index 8693308ca..9169c7d7b 100644 --- a/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt +++ b/examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt @@ -14,14 +14,14 @@ fun DMatrixContext.simple() { m1.transpose() + m2 } -object D5: Dimension{ +object D5 : Dimension { override val dim: UInt = 5u } fun DMatrixContext.custom() { - val m1 = produce { i, j -> (i+j).toDouble() } - val m2 = produce { i, j -> (i-j).toDouble() } - val m3 = produce { i, j -> (i-j).toDouble() } + val m1 = produce { i, j -> (i + j).toDouble() } + val m2 = produce { i, j -> (i - j).toDouble() } + val m3 = produce { i, j -> (i - j).toDouble() } (m1 dot m2) + m3 } -- 2.34.1 From 2d714523412b80e9047d11194aa6b27ab73117b1 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Mon, 9 Dec 2019 20:56:09 +0300 Subject: [PATCH 12/22] Workaround for JS reflection for dimenstions --- .../kotlin/scientifik/kmath/dimensions/Dimensions.kt | 6 ++---- .../kotlin/scientifik/kmath/dimensions/Wrappers.kt | 10 ++++++++-- .../jsMain/kotlin/scientifik/kmath/dimensions/dim.kt | 5 +++++ .../jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt | 4 ++++ .../scientifik/kmath/dimensions/DMatrixContextTest.kt | 0 5 files changed, 19 insertions(+), 6 deletions(-) create mode 100644 kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt create mode 100644 kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt rename kmath-dimensions/src/{commonTest => jvmTest}/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt (100%) diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt index 87511b885..2eef69a0c 100644 --- a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt @@ -19,13 +19,11 @@ interface Dimension { } } } - - inline fun dim(): UInt { - return D::class.objectInstance?.dim ?: error("Dimension object must be a singleton") - } } } +expect inline fun Dimension.Companion.dim(): UInt + object D1 : Dimension { override val dim: UInt = 1u } diff --git a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt index 9985888b9..c78018753 100644 --- a/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt +++ b/kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt @@ -66,7 +66,8 @@ interface DPoint : Point { /** * Dimension-safe point wrapper */ -inline class DPointWrapper(val point: Point) : DPoint { +inline class DPointWrapper(val point: Point) : + DPoint { override val size: Int get() = point.size override fun get(index: Int): T = point[index] @@ -94,7 +95,12 @@ inline class DMatrixContext>(val context: GenericMatrixCon inline fun point(noinline initializer: (Int) -> T): DPoint { val size = Dimension.dim() - return DPoint.coerceUnsafe(context.point(size.toInt(), initializer)) + return DPoint.coerceUnsafe( + context.point( + size.toInt(), + initializer + ) + ) } inline infix fun DMatrix.dot( diff --git a/kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt b/kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt new file mode 100644 index 000000000..557234d84 --- /dev/null +++ b/kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt @@ -0,0 +1,5 @@ +package scientifik.kmath.dimensions + +actual inline fun Dimension.Companion.dim(): UInt { + TODO("KClass::objectInstance does not work") +} \ No newline at end of file diff --git a/kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt b/kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt new file mode 100644 index 000000000..a16dd2266 --- /dev/null +++ b/kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt @@ -0,0 +1,4 @@ +package scientifik.kmath.dimensions + +actual inline fun Dimension.Companion.dim(): UInt = + D::class.objectInstance?.dim ?: error("Dimension object must be a singleton") \ No newline at end of file diff --git a/kmath-dimensions/src/commonTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt b/kmath-dimensions/src/jvmTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt similarity index 100% rename from kmath-dimensions/src/commonTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt rename to kmath-dimensions/src/jvmTest/kotlin/scientifik/kmath/dimensions/DMatrixContextTest.kt -- 2.34.1 From 74653e74c6a6974db3cf51b658e289173a98fbb4 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 11 Dec 2019 15:06:18 +0300 Subject: [PATCH 13/22] Direct memory for file reading. --- .../kotlin/scientifik/memory/ByteBufferMemory.kt | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt index 30a87228a..df2e9847a 100644 --- a/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt +++ b/kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt @@ -1,6 +1,10 @@ package scientifik.memory import java.nio.ByteBuffer +import java.nio.channels.FileChannel +import java.nio.file.Files +import java.nio.file.Path +import java.nio.file.StandardOpenOption /** @@ -11,7 +15,7 @@ actual fun Memory.Companion.allocate(length: Int): Memory { return ByteBufferMemory(buffer) } -class ByteBufferMemory( +private class ByteBufferMemory( val buffer: ByteBuffer, val startOffset: Int = 0, override val size: Int = buffer.limit() @@ -90,4 +94,13 @@ class ByteBufferMemory( } override fun writer(): MemoryWriter = writer +} + +/** + * Use direct memory-mapped buffer from file to read something and close it afterwards. + */ +fun Path.readAsMemory(position: Long = 0, size: Long = Files.size(this), block: Memory.() -> R): R { + return FileChannel.open(this, StandardOpenOption.READ).use { + ByteBufferMemory(it.map(FileChannel.MapMode.READ_ONLY, position, size)).block() + } } \ No newline at end of file -- 2.34.1 From 396b31d106e9fd7eea52c45c5c278d13bab066e7 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 12 Jan 2020 10:49:42 +0300 Subject: [PATCH 14/22] Linear interpolation --- README.md | 11 ++- build.gradle.kts | 2 +- .../scientifik/kmath/operations/Algebra.kt | 6 +- .../kmath/operations/AlgebraExtensions.kt | 13 ++- .../kmath/operations/OptionalOperations.kt | 4 +- kmath-functions/build.gradle.kts | 12 --- .../scientifik.kmath.functions/Piecewise.kt | 2 - .../scientifik.kmath.functions/Polynomial.kt | 94 ++++++++++++++++++- .../scientifik.kmath.functions/functions.kt | 53 ++++------- .../kmath/interpolation/Interpolator.kt | 18 +++- .../kmath/interpolation/LinearInterpolator.kt | 22 ++++- .../interpolation/LinearInterpolatorTest.kt | 25 ++++- settings.gradle.kts | 11 ++- 13 files changed, 198 insertions(+), 75 deletions(-) delete mode 100644 kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt diff --git a/README.md b/README.md index 0df279889..34761e838 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,12 @@ -Bintray: [ ![Download](https://api.bintray.com/packages/mipt-npm/scientifik/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) - -Bintray-dev: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) - +[![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub) [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382) +![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg) + +Bintray: [ ![Download](https://api.bintray.com/packages/mipt-npm/scientifik/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) + +Bintray-dev: [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion) + # KMath Could be pronounced as `key-math`. The Kotlin MATHematics library is intended as a Kotlin-based analog to Python's `numpy` library. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. diff --git a/build.gradle.kts b/build.gradle.kts index a63966dfe..f3a3c13d4 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,5 +1,5 @@ plugins { - id("scientifik.publish") version "0.2.6" apply false + id("scientifik.publish") version "0.3.1" apply false } val kmathVersion by extra("0.1.4-dev-1") diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt index 0ed769db8..485185526 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt @@ -6,14 +6,14 @@ annotation class KMathContext /** * Marker interface for any algebra */ -interface Algebra +interface Algebra -inline operator fun T.invoke(block: T.() -> R): R = run(block) +inline operator fun , R> T.invoke(block: T.() -> R): R = run(block) /** * Space-like operations without neutral element */ -interface SpaceOperations : Algebra { +interface SpaceOperations : Algebra { /** * Addition operation for two context elements */ diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt index 4e8dbf36d..1e0453f08 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt @@ -1,4 +1,13 @@ package scientifik.kmath.operations -fun Space.sum(data : Iterable): T = data.fold(zero) { left, right -> add(left,right) } -fun Space.sum(data : Sequence): T = data.fold(zero) { left, right -> add(left, right) } \ No newline at end of file +fun Space.sum(data: Iterable): T = data.fold(zero) { left, right -> add(left, right) } +fun Space.sum(data: Sequence): T = data.fold(zero) { left, right -> add(left, right) } + +//TODO optimized power operation +fun RingOperations.power(arg: T, power: Int): T { + var res = arg + repeat(power - 1) { + res *= arg + } + return res +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt index 2a77e36ef..bd83932e7 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt @@ -29,7 +29,7 @@ fun >> ctg(arg: T): T = arg.conte /** * A context extension to include power operations like square roots, etc */ -interface PowerOperations { +interface PowerOperations : Algebra { fun power(arg: T, pow: Number): T fun sqrt(arg: T) = power(arg, 0.5) @@ -42,7 +42,7 @@ fun >> sqr(arg: T): T = arg pow 2.0 /* Exponential */ -interface ExponentialOperations { +interface ExponentialOperations: Algebra { fun exp(arg: T): T fun ln(arg: T): T } diff --git a/kmath-functions/build.gradle.kts b/kmath-functions/build.gradle.kts index 373d9b8ac..4c158a32e 100644 --- a/kmath-functions/build.gradle.kts +++ b/kmath-functions/build.gradle.kts @@ -1,23 +1,11 @@ plugins { id("scientifik.mpp") - //id("scientifik.atomic") } kotlin.sourceSets { commonMain { dependencies { api(project(":kmath-core")) - api("org.jetbrains.kotlinx:kotlinx-coroutines-core-common:${Scientifik.coroutinesVersion}") - } - } - jvmMain { - dependencies { - api("org.jetbrains.kotlinx:kotlinx-coroutines-core:${Scientifik.coroutinesVersion}") - } - } - jsMain { - dependencies { - api("org.jetbrains.kotlinx:kotlinx-coroutines-core-js:${Scientifik.coroutinesVersion}") } } } diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt deleted file mode 100644 index d1263e58b..000000000 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt +++ /dev/null @@ -1,2 +0,0 @@ -package scientifik.kmath.functions - diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt index 855ed26cd..a05462863 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt @@ -1,4 +1,94 @@ package scientifik.kmath.functions -interface Polynomial { -} \ No newline at end of file +import scientifik.kmath.operations.RealField +import scientifik.kmath.operations.Ring +import scientifik.kmath.operations.Space +import kotlin.jvm.JvmName +import kotlin.math.max +import kotlin.math.pow + +/** + * Polynomial coefficients without fixation on specific context they are applied to + * @param coefficients constant is the leftmost coefficient + */ +inline class Polynomial(val coefficients: List) { + constructor(vararg coefficients: T) : this(coefficients.toList()) +} + +fun Polynomial.value() = + coefficients.reduceIndexed { index: Int, acc: Double, d: Double -> acc + d.pow(index) } + + +fun > Polynomial.value(ring: C, arg: T): T = ring.run { + if( coefficients.isEmpty()) return@run zero + var res = coefficients.first() + var powerArg = arg + for( index in 1 until coefficients.size){ + res += coefficients[index]*powerArg + //recalculating power on each step to avoid power costs on long polynomials + powerArg *= arg + } + return@run res +} + +/** + * Represent a polynomial as a context-dependent function + */ +fun > Polynomial.asMathFunction(): MathFunction = object : MathFunction { + override fun C.invoke(arg: T): T = value(this, arg) +} + +/** + * Represent the polynomial as a regular context-less function + */ +fun > Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } + +@JvmName("asRealUFunction") +fun Polynomial.asFunction(): (Double) -> Double = asFunction(RealField) + +/** + * An algebra for polynomials + */ +class PolynomialSpace>(val ring: C) : Space> { + + override fun add(a: Polynomial, b: Polynomial): Polynomial { + val dim = max(a.coefficients.size, b.coefficients.size) + ring.run { + return Polynomial(List(dim) { index -> + a.coefficients.getOrElse(index) { zero } + b.coefficients.getOrElse(index) { zero } + }) + } + } + + override fun multiply(a: Polynomial, k: Number): Polynomial { + ring.run { + return Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * k }) + } + } + + override val zero: Polynomial = Polynomial(emptyList()) + + operator fun Polynomial.invoke(arg: T): T = value(ring, arg) +} + +fun , R> C.polynomial(block: PolynomialSpace.() -> R): R { + return PolynomialSpace(this).run(block) +} + +class PiecewisePolynomial> internal constructor( + val lowerBoundary: T, + val pieces: List>> +) + +private fun > PiecewisePolynomial.findPiece(arg: T): Polynomial? { + if (arg < lowerBoundary || arg > pieces.last().first) return null + return pieces.first { arg < it.first }.second +} + +/** + * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. + */ +fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = + findPiece(arg)?.value(ring, arg) + +fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/functions.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/functions.kt index c26443926..2b822b3ba 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/functions.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/functions.kt @@ -1,54 +1,33 @@ -package scientifik.kmath.misc +package scientifik.kmath.functions +import scientifik.kmath.operations.Algebra import scientifik.kmath.operations.RealField -import scientifik.kmath.operations.SpaceOperations -import kotlin.jvm.JvmName /** * A regular function that could be called only inside specific algebra context + * @param T source type + * @param C source algebra constraint + * @param R result type */ -interface UFunction> { - operator fun C.invoke(arg: T): T +interface MathFunction, R> { + operator fun C.invoke(arg: T): R } +fun MathFunction.invoke(arg: Double): R = RealField.invoke(arg) + /** - * A suspendable univariate function defined in algebraic context + * A suspendable function defined in algebraic context */ -interface USFunction> { - suspend operator fun C.invoke(arg: T): T +interface SuspendableMathFunction, R> { + suspend operator fun C.invoke(arg: T): R } -suspend fun USFunction.invoke(arg: Double) = RealField.invoke(arg) +suspend fun SuspendableMathFunction.invoke(arg: Double) = RealField.invoke(arg) -interface MFunction> { - /** - * The input dimension of the function - */ - val dimension: UInt - - operator fun C.invoke(vararg args: T): T -} - /** - * A suspendable multivariate (N->1) function defined on algebraic context + * A parametric function with parameter */ -interface MSFunction> { - /** - * The input dimension of the function - */ - val dimension: UInt - - suspend operator fun C.invoke(vararg args: T): T -} - -suspend fun MSFunction.invoke(args: DoubleArray) = RealField.invoke(*args.toTypedArray()) -@JvmName("varargInvoke") -suspend fun MSFunction.invoke(vararg args: Double) = RealField.invoke(*args.toTypedArray()) - -/** - * A suspendable parametric function with parameter - */ -interface PSFunction> { - suspend operator fun C.invoke(arg: T, parameter: P): T +interface ParametricFunction> { + operator fun C.invoke(arg: T, parameter: P): T } diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt index 0b1ca7795..1df9df5d7 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt @@ -1,7 +1,21 @@ package scientifik.kmath.interpolation -import scientifik.kmath.functions.MathFunction +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.value +import scientifik.kmath.operations.Ring interface Interpolator { - fun interpolate(points: Collection>): MathFunction + fun interpolate(points: Collection>): (X) -> Y +} + +interface PolynomialInterpolator> : Interpolator { + val algebra: Ring + + fun getDefaultValue(): T = error("Out of bounds") + + fun interpolatePolynomials(points: Collection>): PiecewisePolynomial + + override fun interpolate(points: Collection>): (T) -> T = { x -> + interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue() + } } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt index 4db9da7b6..2fc8aaccc 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt @@ -1,4 +1,24 @@ package scientifik.kmath.interpolation -class LinearInterpolator { +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.Polynomial +import scientifik.kmath.operations.Field + +/** + * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java + */ +class LinearInterpolator>(override val algebra: Field) : PolynomialInterpolator { + + override fun interpolatePolynomials(points: Collection>): PiecewisePolynomial = algebra.run { + //sorting points + val sorted = points.sortedBy { it.first } + + val pairs: List>> = (0 until points.size - 1).map { i -> + val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first) + val const = sorted[i].second - slope * sorted[i].first + sorted[i + 1].first to Polynomial(const, slope) + } + + return PiecewisePolynomial(sorted.first().first, pairs) + } } \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt index 34202b388..71303107e 100644 --- a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt +++ b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt @@ -1,5 +1,26 @@ package scientifik.kmath.interpolation -import org.junit.Assert.* +import scientifik.kmath.functions.asFunction +import scientifik.kmath.operations.RealField +import kotlin.test.Test +import kotlin.test.assertEquals -class LinearInterpolatorTest \ No newline at end of file + +class LinearInterpolatorTest { + @Test + fun testInterpolation() { + val data = listOf( + 0.0 to 0.0, + 1.0 to 1.0, + 2.0 to 3.0, + 3.0 to 4.0 + ) + val polynomial = LinearInterpolator(RealField).interpolatePolynomials(data) + val function = polynomial.asFunction(RealField) + +// assertEquals(null, function(-1.0)) +// assertEquals(0.5, function(0.5)) + assertEquals(2.0, function(1.5)) + assertEquals(3.0, function(2.0)) + } +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index e85b32fd2..fdb140541 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -1,10 +1,10 @@ pluginManagement { plugins { - id("scientifik.mpp") version "0.2.5" - id("scientifik.jvm") version "0.2.5" - id("scientifik.atomic") version "0.2.5" - id("scientifik.publish") version "0.2.5" + id("scientifik.mpp") version "0.3.1" + id("scientifik.jvm") version "0.3.1" + id("scientifik.atomic") version "0.3.1" + id("scientifik.publish") version "0.3.1" } repositories { @@ -19,7 +19,7 @@ pluginManagement { resolutionStrategy { eachPlugin { when (requested.id.id) { - "scientifik.mpp", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}") + "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}") } } } @@ -29,6 +29,7 @@ rootProject.name = "kmath" include( ":kmath-memory", ":kmath-core", + ":kmath-functions", // ":kmath-io", ":kmath-coroutines", ":kmath-histograms", -- 2.34.1 From d56b4148bef72f6cffe39e7426e47c5bdbed6240 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 12 Jan 2020 20:51:16 +0300 Subject: [PATCH 15/22] Generalized linear interpolation --- .../scientifik.kmath.functions/Piecewise.kt | 55 +++++++++++++++++++ .../scientifik.kmath.functions/Polynomial.kt | 31 ++--------- .../kmath/interpolation/LinearInterpolator.kt | 16 ++++-- .../interpolation/LinearInterpolatorTest.kt | 4 +- 4 files changed, 71 insertions(+), 35 deletions(-) create mode 100644 kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt new file mode 100644 index 000000000..c05660d52 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt @@ -0,0 +1,55 @@ +package scientifik.kmath.functions + +import scientifik.kmath.operations.Ring + +interface Piecewise { + fun findPiece(arg: T): R? +} + +interface PiecewisePolynomial : Piecewise> + +/** + * Ordered list of pieces in piecewise function + */ +class OrderedPiecewisePolynomial>(left: T) : PiecewisePolynomial { + + private val delimiters: ArrayList = arrayListOf(left) + private val pieces: ArrayList> = ArrayList() + + /** + * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece) + * @param right new rightmost position. If is less then current rightmost position, a error is thrown. + */ + fun putRight(right: T, piece: Polynomial) { + require(right > delimiters.last()) { "New delimiter should be to the right of old one" } + delimiters.add(right) + pieces.add(piece) + } + + fun putLeft(left: T, piece: Polynomial) { + require(left < delimiters.first()) { "New delimiter should be to the left of old one" } + delimiters.add(0, left) + pieces.add(0, piece) + } + + override fun findPiece(arg: T): Polynomial? { + if (arg < delimiters.first() || arg >= delimiters.last()) { + return null + } else { + for (index in 1 until delimiters.size) { + if (arg < delimiters[index]) { + return pieces[index - 1] + } + } + error("Piece not found") + } + } +} + +/** + * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. + */ +fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = + findPiece(arg)?.value(ring, arg) + +fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt index a05462863..bd5062661 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt @@ -1,9 +1,7 @@ package scientifik.kmath.functions -import scientifik.kmath.operations.RealField import scientifik.kmath.operations.Ring import scientifik.kmath.operations.Space -import kotlin.jvm.JvmName import kotlin.math.max import kotlin.math.pow @@ -20,11 +18,11 @@ fun Polynomial.value() = fun > Polynomial.value(ring: C, arg: T): T = ring.run { - if( coefficients.isEmpty()) return@run zero + if (coefficients.isEmpty()) return@run zero var res = coefficients.first() var powerArg = arg - for( index in 1 until coefficients.size){ - res += coefficients[index]*powerArg + for (index in 1 until coefficients.size) { + res += coefficients[index] * powerArg //recalculating power on each step to avoid power costs on long polynomials powerArg *= arg } @@ -43,9 +41,6 @@ fun > Polynomial.asMathFunction(): MathFunction> Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } -@JvmName("asRealUFunction") -fun Polynomial.asFunction(): (Double) -> Double = asFunction(RealField) - /** * An algebra for polynomials */ @@ -73,22 +68,4 @@ class PolynomialSpace>(val ring: C) : Space> fun , R> C.polynomial(block: PolynomialSpace.() -> R): R { return PolynomialSpace(this).run(block) -} - -class PiecewisePolynomial> internal constructor( - val lowerBoundary: T, - val pieces: List>> -) - -private fun > PiecewisePolynomial.findPiece(arg: T): Polynomial? { - if (arg < lowerBoundary || arg > pieces.last().first) return null - return pieces.first { arg < it.first }.second -} - -/** - * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. - */ -fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = - findPiece(arg)?.value(ring, arg) - -fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt index 2fc8aaccc..720f1080b 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt @@ -1,5 +1,6 @@ package scientifik.kmath.interpolation +import scientifik.kmath.functions.OrderedPiecewisePolynomial import scientifik.kmath.functions.PiecewisePolynomial import scientifik.kmath.functions.Polynomial import scientifik.kmath.operations.Field @@ -10,15 +11,18 @@ import scientifik.kmath.operations.Field class LinearInterpolator>(override val algebra: Field) : PolynomialInterpolator { override fun interpolatePolynomials(points: Collection>): PiecewisePolynomial = algebra.run { + require(points.isNotEmpty()) { "Point array should not be empty" } + //sorting points val sorted = points.sortedBy { it.first } - val pairs: List>> = (0 until points.size - 1).map { i -> - val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first) - val const = sorted[i].second - slope * sorted[i].first - sorted[i + 1].first to Polynomial(const, slope) + return@run OrderedPiecewisePolynomial(points.first().first).apply { + for (i in 0 until points.size - 1) { + val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first) + val const = sorted[i].second - slope * sorted[i].first + val polynomial = Polynomial(const, slope) + putRight(sorted[i + 1].first, polynomial) + } } - - return PiecewisePolynomial(sorted.first().first, pairs) } } \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt index 71303107e..6d35d19f1 100644 --- a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt +++ b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt @@ -18,8 +18,8 @@ class LinearInterpolatorTest { val polynomial = LinearInterpolator(RealField).interpolatePolynomials(data) val function = polynomial.asFunction(RealField) -// assertEquals(null, function(-1.0)) -// assertEquals(0.5, function(0.5)) + assertEquals(null, function(-1.0)) + assertEquals(0.5, function(0.5)) assertEquals(2.0, function(1.5)) assertEquals(3.0, function(2.0)) } -- 2.34.1 From 73f40105c47630b544dd44157233b3f54e00d7d6 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 12 Feb 2020 21:57:21 +0300 Subject: [PATCH 16/22] Interpolation API --- kmath-commons/build.gradle.kts | 1 + .../scientifik/kmath/structures/Buffers.kt | 2 + .../kmath/functions}/Piecewise.kt | 8 +- .../kmath/functions}/Polynomial.kt | 6 +- .../kmath/functions}/functions.kt | 0 .../kmath/interpolation/Interpolator.kt | 6 +- .../kmath/interpolation/LinearInterpolator.kt | 16 +- .../kmath/interpolation/LoessInterpolator.kt | 296 ++++++++++++++++++ .../kmath/interpolation/SplineInterpolator.kt | 58 ++++ .../kmath/interpolation/XYPointSet.kt | 54 ++++ 10 files changed, 430 insertions(+), 17 deletions(-) rename kmath-functions/src/commonMain/kotlin/{scientifik.kmath.functions => scientifik/kmath/functions}/Piecewise.kt (85%) rename kmath-functions/src/commonMain/kotlin/{scientifik.kmath.functions => scientifik/kmath/functions}/Polynomial.kt (94%) rename kmath-functions/src/commonMain/kotlin/{scientifik.kmath.functions => scientifik/kmath/functions}/functions.kt (100%) create mode 100644 kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt create mode 100644 kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt create mode 100644 kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt diff --git a/kmath-commons/build.gradle.kts b/kmath-commons/build.gradle.kts index 9ac0f6e95..5ce1b935a 100644 --- a/kmath-commons/build.gradle.kts +++ b/kmath-commons/build.gradle.kts @@ -8,5 +8,6 @@ dependencies { api(project(":kmath-core")) api(project(":kmath-coroutines")) api(project(":kmath-prob")) + api(project(":kmath-functions")) api("org.apache.commons:commons-math3:3.6.1") } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt index c63384fb9..f02fd8dd0 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt @@ -73,6 +73,8 @@ fun Buffer.asSequence(): Sequence = Sequence(::iterator) fun Buffer.asIterable(): Iterable = asSequence().asIterable() +val Buffer<*>.indices: IntRange get() = IntRange(0, size - 1) + interface MutableBuffer : Buffer { operator fun set(index: Int, value: T) diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt similarity index 85% rename from kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt rename to kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt index c05660d52..16f8aa12b 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt @@ -6,14 +6,16 @@ interface Piecewise { fun findPiece(arg: T): R? } -interface PiecewisePolynomial : Piecewise> +interface PiecewisePolynomial : + Piecewise> /** * Ordered list of pieces in piecewise function */ -class OrderedPiecewisePolynomial>(left: T) : PiecewisePolynomial { +class OrderedPiecewisePolynomial>(delimeter: T) : + PiecewisePolynomial { - private val delimiters: ArrayList = arrayListOf(left) + private val delimiters: ArrayList = arrayListOf(delimeter) private val pieces: ArrayList> = ArrayList() /** diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt similarity index 94% rename from kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt rename to kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt index bd5062661..b747b521d 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt @@ -32,7 +32,8 @@ fun > Polynomial.value(ring: C, arg: T): T = ring.run { /** * Represent a polynomial as a context-dependent function */ -fun > Polynomial.asMathFunction(): MathFunction = object : MathFunction { +fun > Polynomial.asMathFunction(): MathFunction = object : + MathFunction { override fun C.invoke(arg: T): T = value(this, arg) } @@ -61,7 +62,8 @@ class PolynomialSpace>(val ring: C) : Space> } } - override val zero: Polynomial = Polynomial(emptyList()) + override val zero: Polynomial = + Polynomial(emptyList()) operator fun Polynomial.invoke(arg: T): T = value(ring, arg) } diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/functions.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt similarity index 100% rename from kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/functions.kt rename to kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt index 1df9df5d7..9b13c92e3 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt @@ -5,7 +5,7 @@ import scientifik.kmath.functions.value import scientifik.kmath.operations.Ring interface Interpolator { - fun interpolate(points: Collection>): (X) -> Y + fun interpolate(points: XYPointSet): (X) -> Y } interface PolynomialInterpolator> : Interpolator { @@ -13,9 +13,9 @@ interface PolynomialInterpolator> : Interpolator { fun getDefaultValue(): T = error("Out of bounds") - fun interpolatePolynomials(points: Collection>): PiecewisePolynomial + fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial - override fun interpolate(points: Collection>): (T) -> T = { x -> + override fun interpolate(points: XYPointSet): (T) -> T = { x -> interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue() } } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt index 720f1080b..9467da421 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt @@ -10,18 +10,16 @@ import scientifik.kmath.operations.Field */ class LinearInterpolator>(override val algebra: Field) : PolynomialInterpolator { - override fun interpolatePolynomials(points: Collection>): PiecewisePolynomial = algebra.run { - require(points.isNotEmpty()) { "Point array should not be empty" } + override fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial = algebra.run { + require(points.size > 0) { "Point array should not be empty" } + insureSorted(points) - //sorting points - val sorted = points.sortedBy { it.first } - - return@run OrderedPiecewisePolynomial(points.first().first).apply { + OrderedPiecewisePolynomial(points.x[0]).apply { for (i in 0 until points.size - 1) { - val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first) - val const = sorted[i].second - slope * sorted[i].first + val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i]) + val const = points.x[i] - slope * points.x[i] val polynomial = Polynomial(const, slope) - putRight(sorted[i + 1].first, polynomial) + putRight(points.x[i + 1], polynomial) } } } diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt new file mode 100644 index 000000000..c29549ed0 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt @@ -0,0 +1,296 @@ +//package scientifik.kmath.interpolation +// +//import scientifik.kmath.functions.PiecewisePolynomial +//import scientifik.kmath.operations.Ring +//import scientifik.kmath.structures.Buffer +//import kotlin.math.abs +//import kotlin.math.sqrt +// +// +///** +// * Original code: https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/LoessInterpolator.java +// */ +//class LoessInterpolator>(override val algebra: Ring) : PolynomialInterpolator { +// /** +// * The bandwidth parameter: when computing the loess fit at +// * a particular point, this fraction of source points closest +// * to the current point is taken into account for computing +// * a least-squares regression. +// * +// * +// * A sensible value is usually 0.25 to 0.5. +// */ +// private var bandwidth = 0.0 +// +// /** +// * The number of robustness iterations parameter: this many +// * robustness iterations are done. +// * +// * +// * A sensible value is usually 0 (just the initial fit without any +// * robustness iterations) to 4. +// */ +// private var robustnessIters = 0 +// +// /** +// * If the median residual at a certain robustness iteration +// * is less than this amount, no more iterations are done. +// */ +// private var accuracy = 0.0 +// +// /** +// * Constructs a new [LoessInterpolator] +// * with a bandwidth of [.DEFAULT_BANDWIDTH], +// * [.DEFAULT_ROBUSTNESS_ITERS] robustness iterations +// * and an accuracy of {#link #DEFAULT_ACCURACY}. +// * See [.LoessInterpolator] for an explanation of +// * the parameters. +// */ +// fun LoessInterpolator() { +// bandwidth = DEFAULT_BANDWIDTH +// robustnessIters = DEFAULT_ROBUSTNESS_ITERS +// accuracy = DEFAULT_ACCURACY +// } +// +// fun LoessInterpolator(bandwidth: Double, robustnessIters: Int) { +// this(bandwidth, robustnessIters, DEFAULT_ACCURACY) +// } +// +// fun LoessInterpolator(bandwidth: Double, robustnessIters: Int, accuracy: Double) { +// if (bandwidth < 0 || +// bandwidth > 1 +// ) { +// throw OutOfRangeException(LocalizedFormats.BANDWIDTH, bandwidth, 0, 1) +// } +// this.bandwidth = bandwidth +// if (robustnessIters < 0) { +// throw NotPositiveException(LocalizedFormats.ROBUSTNESS_ITERATIONS, robustnessIters) +// } +// this.robustnessIters = robustnessIters +// this.accuracy = accuracy +// } +// +// fun interpolate( +// xval: DoubleArray, +// yval: DoubleArray +// ): PolynomialSplineFunction { +// return SplineInterpolator().interpolate(xval, smooth(xval, yval)) +// } +// +// fun XYZPointSet.smooth(): XYPointSet { +// checkAllFiniteReal(x) +// checkAllFiniteReal(y) +// checkAllFiniteReal(z) +// MathArrays.checkOrder(xval) +// if (size == 1) { +// return doubleArrayOf(y[0]) +// } +// if (size == 2) { +// return doubleArrayOf(y[0], y[1]) +// } +// val bandwidthInPoints = (bandwidth * size).toInt() +// if (bandwidthInPoints < 2) { +// throw NumberIsTooSmallException( +// LocalizedFormats.BANDWIDTH, +// bandwidthInPoints, 2, true +// ) +// } +// val res = DoubleArray(size) +// val residuals = DoubleArray(size) +// val sortedResiduals = DoubleArray(size) +// val robustnessWeights = DoubleArray(size) +// // Do an initial fit and 'robustnessIters' robustness iterations. +// // This is equivalent to doing 'robustnessIters+1' robustness iterations +// // starting with all robustness weights set to 1. +// Arrays.fill(robustnessWeights, 1.0) +// for (iter in 0..robustnessIters) { +// val bandwidthInterval = intArrayOf(0, bandwidthInPoints - 1) +// // At each x, compute a local weighted linear regression +// for (i in 0 until size) { +//// val x = x[i] +// // Find out the interval of source points on which +// // a regression is to be made. +// if (i > 0) { +// updateBandwidthInterval(x, z, i, bandwidthInterval) +// } +// val ileft = bandwidthInterval[0] +// val iright = bandwidthInterval[1] +// // Compute the point of the bandwidth interval that is +// // farthest from x +// val edge: Int +// edge = if (x[i] - x[ileft] > x[iright] - x[i]) { +// ileft +// } else { +// iright +// } +// // Compute a least-squares linear fit weighted by +// // the product of robustness weights and the tricube +// // weight function. +// // See http://en.wikipedia.org/wiki/Linear_regression +// // (section "Univariate linear case") +// // and http://en.wikipedia.org/wiki/Weighted_least_squares +// // (section "Weighted least squares") +// var sumWeights = 0.0 +// var sumX = 0.0 +// var sumXSquared = 0.0 +// var sumY = 0.0 +// var sumXY = 0.0 +// val denom: Double = abs(1.0 / (x[edge] - x[i])) +// for (k in ileft..iright) { +// val xk = x[k] +// val yk = y[k] +// val dist = if (k < i) x - xk else xk - x[i] +// val w = tricube(dist * denom) * robustnessWeights[k] * z[k] +// val xkw = xk * w +// sumWeights += w +// sumX += xkw +// sumXSquared += xk * xkw +// sumY += yk * w +// sumXY += yk * xkw +// } +// val meanX = sumX / sumWeights +// val meanY = sumY / sumWeights +// val meanXY = sumXY / sumWeights +// val meanXSquared = sumXSquared / sumWeights +// val beta: Double +// beta = if (sqrt(abs(meanXSquared - meanX * meanX)) < accuracy) { +// 0.0 +// } else { +// (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX) +// } +// val alpha = meanY - beta * meanX +// res[i] = beta * x[i] + alpha +// residuals[i] = abs(y[i] - res[i]) +// } +// // No need to recompute the robustness weights at the last +// // iteration, they won't be needed anymore +// if (iter == robustnessIters) { +// break +// } +// // Recompute the robustness weights. +// // Find the median residual. +// // An arraycopy and a sort are completely tractable here, +// // because the preceding loop is a lot more expensive +// java.lang.System.arraycopy(residuals, 0, sortedResiduals, 0, size) +// Arrays.sort(sortedResiduals) +// val medianResidual = sortedResiduals[size / 2] +// if (abs(medianResidual) < accuracy) { +// break +// } +// for (i in 0 until size) { +// val arg = residuals[i] / (6 * medianResidual) +// if (arg >= 1) { +// robustnessWeights[i] = 0.0 +// } else { +// val w = 1 - arg * arg +// robustnessWeights[i] = w * w +// } +// } +// } +// return res +// } +// +// fun smooth(xval: DoubleArray, yval: DoubleArray): DoubleArray { +// if (xval.size != yval.size) { +// throw DimensionMismatchException(xval.size, yval.size) +// } +// val unitWeights = DoubleArray(xval.size) +// Arrays.fill(unitWeights, 1.0) +// return smooth(xval, yval, unitWeights) +// } +// +// /** +// * Given an index interval into xval that embraces a certain number of +// * points closest to `xval[i-1]`, update the interval so that it +// * embraces the same number of points closest to `xval[i]`, +// * ignoring zero weights. +// * +// * @param xval Arguments array. +// * @param weights Weights array. +// * @param i Index around which the new interval should be computed. +// * @param bandwidthInterval a two-element array {left, right} such that: +// * `(left==0 or xval[i] - xval[left-1] > xval[right] - xval[i])` +// * and +// * `(right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])`. +// * The array will be updated. +// */ +// private fun updateBandwidthInterval( +// xval: Buffer, weights: Buffer, +// i: Int, +// bandwidthInterval: IntArray +// ) { +// val left = bandwidthInterval[0] +// val right = bandwidthInterval[1] +// // The right edge should be adjusted if the next point to the right +// // is closer to xval[i] than the leftmost point of the current interval +// val nextRight = nextNonzero(weights, right) +// if (nextRight < xval.size && xval[nextRight] - xval[i] < xval[i] - xval[left]) { +// val nextLeft = nextNonzero(weights, bandwidthInterval[0]) +// bandwidthInterval[0] = nextLeft +// bandwidthInterval[1] = nextRight +// } +// } +// +// /** +// * Return the smallest index `j` such that +// * `j > i && (j == weights.length || weights[j] != 0)`. +// * +// * @param weights Weights array. +// * @param i Index from which to start search. +// * @return the smallest compliant index. +// */ +// private fun nextNonzero(weights: Buffer, i: Int): Int { +// var j = i + 1 +// while (j < weights.size && weights[j] == 0.0) { +// ++j +// } +// return j +// } +// +// /** +// * Compute the +// * [tricube](http://en.wikipedia.org/wiki/Local_regression#Weight_function) +// * weight function +// * +// * @param x Argument. +// * @return `(1 - |x|3)3` for |x| < 1, 0 otherwise. +// */ +// private fun tricube(x: Double): Double { +// val absX: Double = FastMath.abs(x) +// if (absX >= 1.0) { +// return 0.0 +// } +// val tmp = 1 - absX * absX * absX +// return tmp * tmp * tmp +// } +// +// /** +// * Check that all elements of an array are finite real numbers. +// * +// * @param values Values array. +// * @throws org.apache.commons.math4.exception.NotFiniteNumberException +// * if one of the values is not a finite real number. +// */ +// private fun checkAllFiniteReal(values: DoubleArray) { +// for (i in values.indices) { +// MathUtils.checkFinite(values[i]) +// } +// } +// +// override fun interpolatePolynomials(points: Collection>): PiecewisePolynomial { +// TODO("not implemented") //To change body of created functions use File | Settings | File Templates. +// } +// +// companion object { +// /** Default value of the bandwidth parameter. */ +// const val DEFAULT_BANDWIDTH = 0.3 +// +// /** Default value of the number of robustness iterations. */ +// const val DEFAULT_ROBUSTNESS_ITERS = 2 +// +// /** +// * Default value for accuracy. +// */ +// const val DEFAULT_ACCURACY = 1e-12 +// } +//} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt new file mode 100644 index 000000000..e1af0c1a2 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt @@ -0,0 +1,58 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.functions.OrderedPiecewisePolynomial +import scientifik.kmath.functions.PiecewisePolynomial +import scientifik.kmath.functions.Polynomial +import scientifik.kmath.operations.Field +import scientifik.kmath.structures.MutableBufferFactory + +/** + * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type specific ones. + * Based on https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java + */ +class SplineInterpolator>( + override val algebra: Field, + val bufferFactory: MutableBufferFactory +) : PolynomialInterpolator { + + //TODO possibly optimize zeroed buffers + + override fun interpolatePolynomials(points: XYPointSet): PiecewisePolynomial = algebra.run { + if (points.size < 3) { + error("Can't use spline interpolator with less than 3 points") + } + insureSorted(points) + + // Number of intervals. The number of data points is n + 1. + val n = points.size - 1 + // Differences between knot points + val h = bufferFactory(points.size) { i -> points.x[i + 1] - points.x[i] } + val mu = bufferFactory(points.size - 1) { zero } + val z = bufferFactory(points.size) { zero } + + for (i in 1 until n) { + val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1] + mu[i] = h[i] / g + z[i] = + (3.0 * (points.y[i + 1] * h[i - 1] - points.x[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) / (h[i - 1] * h[i]) + - h[i - 1] * z[i - 1]) / g + } + + // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) + + OrderedPiecewisePolynomial(points.x[points.size - 1]).apply { + var cOld = zero + for (j in n - 1 downTo 0) { + val c = z[j] - mu[j] * cOld + val a = points.y[j] + val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0 + val d = (cOld - c) / (3.0 * h[j]) + val polynomial = Polynomial(a, b, c, d) + cOld = c + putLeft(points.x[j], polynomial) + } + } + + } + +} diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt new file mode 100644 index 000000000..d8e10b880 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt @@ -0,0 +1,54 @@ +package scientifik.kmath.interpolation + +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.Structure2D + +interface XYPointSet { + val size: Int + val x: Buffer + val y: Buffer +} + +interface XYZPointSet : XYPointSet { + val z: Buffer +} + +internal fun > insureSorted(points: XYPointSet) { + for (i in 0 until points.size - 1) { + if (points.x[i + 1] <= points.x[i]) error("Input data is not sorted at index $i") + } +} + +class NDStructureColumn(val structure: Structure2D, val column: Int) : Buffer { + init { + require(column < structure.colNum) { "Column index is outside of structure column range" } + } + + override val size: Int get() = structure.rowNum + + override fun get(index: Int): T = structure[index, column] + + override fun iterator(): Iterator = sequence { + repeat(size) { + yield(get(it)) + } + }.iterator() +} + +class BufferXYPointSet(override val x: Buffer, override val y: Buffer) : XYPointSet { + init { + require(x.size == y.size) { "Sizes of x and y buffers should be the same" } + } + + override val size: Int + get() = x.size +} + +fun Structure2D.asXYPointSet(): XYPointSet { + require(shape[1] == 2) { "Structure second dimension should be of size 2" } + return object : XYPointSet { + override val size: Int get() = this@asXYPointSet.shape[0] + override val x: Buffer get() = NDStructureColumn(this@asXYPointSet, 0) + override val y: Buffer get() = NDStructureColumn(this@asXYPointSet, 1) + } +} \ No newline at end of file -- 2.34.1 From 068b90e7ab9e116bf6668dd960111f71da00b35e Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 12 Feb 2020 22:26:25 +0300 Subject: [PATCH 17/22] Add different ways to provide data for XY interpolator --- .../kmath/interpolation/Interpolator.kt | 24 +++++++++++++++++++ .../kmath/interpolation/LoessInterpolator.kt | 2 +- 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt index 9b13c92e3..8d83e4198 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt @@ -3,6 +3,8 @@ package scientifik.kmath.interpolation import scientifik.kmath.functions.PiecewisePolynomial import scientifik.kmath.functions.value import scientifik.kmath.operations.Ring +import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.asBuffer interface Interpolator { fun interpolate(points: XYPointSet): (X) -> Y @@ -18,4 +20,26 @@ interface PolynomialInterpolator> : Interpolator { override fun interpolate(points: XYPointSet): (T) -> T = { x -> interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue() } +} + +fun > PolynomialInterpolator.interpolatePolynomials( + x: Buffer, + y: Buffer +): PiecewisePolynomial { + val pointSet = BufferXYPointSet(x, y) + return interpolatePolynomials(pointSet) +} + +fun > PolynomialInterpolator.interpolatePolynomials( + data: Map +): PiecewisePolynomial { + val pointSet = BufferXYPointSet(data.keys.toList().asBuffer(), data.values.toList().asBuffer()) + return interpolatePolynomials(pointSet) +} + +fun > PolynomialInterpolator.interpolatePolynomials( + data: List> +): PiecewisePolynomial { + val pointSet = BufferXYPointSet(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer()) + return interpolatePolynomials(pointSet) } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt index c29549ed0..6707bd8bc 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt @@ -1,4 +1,4 @@ -//package scientifik.kmath.interpolation +package scientifik.kmath.interpolation // //import scientifik.kmath.functions.PiecewisePolynomial //import scientifik.kmath.operations.Ring -- 2.34.1 From e00a66ca428ccdd18b8836f6dc9c9fb200a494c4 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Thu, 13 Feb 2020 16:21:41 +0300 Subject: [PATCH 18/22] Fix linear interpolator const --- .../scientifik/kmath/interpolation/LinearInterpolator.kt | 2 +- .../scientifik/kmath/interpolation/LinearInterpolatorTest.kt | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt index 9467da421..98beb4391 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt @@ -17,7 +17,7 @@ class LinearInterpolator>(override val algebra: Field) : Po OrderedPiecewisePolynomial(points.x[0]).apply { for (i in 0 until points.size - 1) { val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i]) - val const = points.x[i] - slope * points.x[i] + val const = points.y[i] - slope * points.x[i] val polynomial = Polynomial(const, slope) putRight(points.x[i + 1], polynomial) } diff --git a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt index 6d35d19f1..23acd835c 100644 --- a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt +++ b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt @@ -1,5 +1,6 @@ package scientifik.kmath.interpolation +import scientifik.kmath.functions.PiecewisePolynomial import scientifik.kmath.functions.asFunction import scientifik.kmath.operations.RealField import kotlin.test.Test @@ -15,7 +16,7 @@ class LinearInterpolatorTest { 2.0 to 3.0, 3.0 to 4.0 ) - val polynomial = LinearInterpolator(RealField).interpolatePolynomials(data) + val polynomial: PiecewisePolynomial = LinearInterpolator(RealField).interpolatePolynomials(data) val function = polynomial.asFunction(RealField) assertEquals(null, function(-1.0)) -- 2.34.1 From 0e898ee1ea6d708c13a17ef7631144ca79195904 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 3 Apr 2020 19:01:58 +0300 Subject: [PATCH 19/22] Add strided matrix dot test to check #82 --- .../scientifik/kmath/linear/MatrixTest.kt | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt index d27aa665c..7d1209963 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt @@ -1,6 +1,8 @@ package scientifik.kmath.linear import scientifik.kmath.structures.Matrix +import scientifik.kmath.structures.NDStructure +import scientifik.kmath.structures.as2D import kotlin.test.Test import kotlin.test.assertEquals @@ -44,4 +46,20 @@ class MatrixTest { val toTenthPower = transitionMatrix pow 10 } + + @Test + fun test2DDot() { + val firstMatrix = NDStructure.auto(2,3){ (i, j) -> (i + j).toDouble() }.as2D() + val secondMatrix = NDStructure.auto(3,2){ (i, j) -> (i + j).toDouble() }.as2D() + MatrixContext.real.run { +// val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() } +// val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() } + val result = firstMatrix dot secondMatrix + assertEquals(2, result.rowNum) + assertEquals(2, result.colNum) + assertEquals(8.0, result[0,1]) + assertEquals(8.0, result[1,0]) + assertEquals(14.0, result[1,1]) + } + } } \ No newline at end of file -- 2.34.1 From 0cb53792b1de5ba5589f28d03e12e554f6ac311c Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 3 Apr 2020 21:44:07 +0300 Subject: [PATCH 20/22] Update plugin and fix median Statistic --- build.gradle.kts | 2 +- .../src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/build.gradle.kts b/build.gradle.kts index f3a3c13d4..1034f3f44 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -1,5 +1,5 @@ plugins { - id("scientifik.publish") version "0.3.1" apply false + id("scientifik.publish") version "0.4.1" apply false } val kmathVersion by extra("0.1.4-dev-1") diff --git a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt index 1af2570b0..804aed089 100644 --- a/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt +++ b/kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt @@ -11,6 +11,7 @@ import scientifik.kmath.coroutines.mapParallel import scientifik.kmath.operations.* import scientifik.kmath.structures.Buffer import scientifik.kmath.structures.asIterable +import scientifik.kmath.structures.asSequence /** * A function, that transforms a buffer of random quantities to some resulting value @@ -83,9 +84,9 @@ class Mean(val space: Space) : ComposableStatistic, T> { /** * Non-composable median */ -class Median(comparator: Comparator) : Statistic { +class Median(private val comparator: Comparator) : Statistic { override suspend fun invoke(data: Buffer): T { - return data.asIterable().toList()[data.size / 2] //TODO check if this is correct + return data.asSequence().sortedWith(comparator).toList()[data.size / 2] //TODO check if this is correct } companion object { -- 2.34.1 From a1dd71a74b2e232ca1b55ca79f7c5ac04085081c Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 3 Apr 2020 22:29:23 +0300 Subject: [PATCH 21/22] update build tools and dependencies --- build.gradle.kts | 2 +- examples/build.gradle.kts | 6 ++++-- settings.gradle.kts | 9 +++++---- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/build.gradle.kts b/build.gradle.kts index 1034f3f44..6c4536b4c 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -2,7 +2,7 @@ plugins { id("scientifik.publish") version "0.4.1" apply false } -val kmathVersion by extra("0.1.4-dev-1") +val kmathVersion by extra("0.1.4-dev-2") val bintrayRepo by extra("scientifik") val githubProject by extra("kmath") diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index 47acaa5ba..8853b78a5 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -4,7 +4,7 @@ import org.jetbrains.kotlin.gradle.tasks.KotlinCompile plugins { java kotlin("jvm") - kotlin("plugin.allopen") version "1.3.61" + kotlin("plugin.allopen") version "1.3.71" id("kotlinx.benchmark") version "0.2.0-dev-7" } @@ -14,6 +14,8 @@ configure { repositories { maven("http://dl.bintray.com/kyonifer/maven") + maven("https://dl.bintray.com/mipt-npm/scientifik") + maven("https://dl.bintray.com/mipt-npm/dev") mavenCentral() } @@ -29,7 +31,7 @@ dependencies { implementation(project(":kmath-viktor")) implementation(project(":kmath-dimensions")) implementation("com.kyonifer:koma-core-ejml:0.12") - implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:${Scientifik.ioVersion}") + implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:0.2.0-npm-dev-6") implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-7") "benchmarksCompile"(sourceSets.main.get().compileClasspath) } diff --git a/settings.gradle.kts b/settings.gradle.kts index d7bcb249f..671db444a 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -1,10 +1,10 @@ pluginManagement { plugins { - id("scientifik.mpp") version "0.3.1" - id("scientifik.jvm") version "0.3.1" - id("scientifik.atomic") version "0.3.1" - id("scientifik.publish") version "0.3.1" + id("scientifik.mpp") version "0.4.1" + id("scientifik.jvm") version "0.4.1" + id("scientifik.atomic") version "0.4.1" + id("scientifik.publish") version "0.4.1" } repositories { @@ -13,6 +13,7 @@ pluginManagement { gradlePluginPortal() maven("https://dl.bintray.com/kotlin/kotlin-eap") maven("https://dl.bintray.com/mipt-npm/scientifik") + maven("https://dl.bintray.com/mipt-npm/dev") maven("https://dl.bintray.com/kotlin/kotlinx") } -- 2.34.1 From cb1156fd9051df20250856cc22ae76cb6cbd01cc Mon Sep 17 00:00:00 2001 From: Andreas Radke Date: Fri, 10 Apr 2020 20:54:10 +0200 Subject: [PATCH 22/22] Fix two simple typos in file names: Rename LinearAlgrebra.kt to LinearAlgebra.kt and BigNumers.kt to BigNumbers.kt --- .../kmath/linear/{LinearAlgrebra.kt => LinearAlgebra.kt} | 0 .../scientifik/kmath/operations/{BigNumers.kt => BigNumbers.kt} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/{LinearAlgrebra.kt => LinearAlgebra.kt} (100%) rename kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/{BigNumers.kt => BigNumbers.kt} (100%) diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt similarity index 100% rename from kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgebra.kt diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumers.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt similarity index 100% rename from kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumers.kt rename to kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/BigNumbers.kt -- 2.34.1