diff --git a/CHANGELOG.md b/CHANGELOG.md index 6733c1211..857ed060b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,7 @@ - `BigInt` operation performance improvement and fixes by @zhelenskiy (#328) - Integration between `MST` and Symja `IExpr` - Complex power +- Separate methods for UInt, Int and Number powers. NaN safety. ### Changed - Exponential operations merged with hyperbolic functions diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt index 64f9b5dff..9203c269e 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt @@ -13,6 +13,7 @@ import space.kscience.kmath.commons.linear.CMLinearSpace import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM import space.kscience.kmath.linear.invoke import space.kscience.kmath.linear.linearSpace +import space.kscience.kmath.multik.multikAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.Buffer import kotlin.random.Random @@ -58,6 +59,16 @@ internal class DotBenchmark { blackhole.consume(matrix1 dot matrix2) } +// @Benchmark +// fun tensorDot(blackhole: Blackhole) = with(Double.tensorAlgebra) { +// blackhole.consume(matrix1 dot matrix2) +// } + + @Benchmark + fun multikDot(blackhole: Blackhole) = with(Double.multikAlgebra) { + blackhole.consume(matrix1 dot matrix2) + } + @Benchmark fun bufferedDot(blackhole: Blackhole) = with(DoubleField.linearSpace(Buffer.Companion::auto)) { blackhole.consume(matrix1 dot matrix2) diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/neuralNetwork.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/neuralNetwork.kt index cdfc06922..1a2a94534 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/tensors/neuralNetwork.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/neuralNetwork.kt @@ -163,7 +163,7 @@ class NeuralNetwork(private val layers: List) { for ((xBatch, yBatch) in iterBatch(xTrain, yTrain)) { train(xBatch, yBatch) } - println("Accuracy:${accuracy(yTrain, predict(xTrain).argMax(1, true))}") + println("Accuracy:${accuracy(yTrain, predict(xTrain).argMax(1, true).asDouble())}") } } @@ -230,7 +230,7 @@ fun main() = BroadcastDoubleTensorAlgebra { val prediction = model.predict(xTest) // process raw prediction via argMax - val predictionLabels = prediction.argMax(1, true) + val predictionLabels = prediction.argMax(1, true).asDouble() // find out accuracy val acc = accuracy(yTest, predictionLabels) diff --git a/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt b/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt index 9d5b1cddd..2eaa17ded 100644 --- a/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt +++ b/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt @@ -18,7 +18,7 @@ import kotlin.contracts.contract */ @OptIn(UnstableKMathAPI::class) public sealed class ComplexFieldOpsND : BufferedFieldOpsND(ComplexField.bufferAlgebra), - ScaleOperations>, ExtendedFieldOps> { + ScaleOperations>, ExtendedFieldOps>, PowerOperations> { override fun StructureND.toBufferND(): BufferND = when (this) { is BufferND -> this @@ -33,9 +33,6 @@ public sealed class ComplexFieldOpsND : BufferedFieldOpsND, value: Double): BufferND = mapInline(a.toBufferND()) { it * value } - override fun power(arg: StructureND, pow: Number): BufferND = - mapInline(arg.toBufferND()) { power(it, pow) } - override fun exp(arg: StructureND): BufferND = mapInline(arg.toBufferND()) { exp(it) } override fun ln(arg: StructureND): BufferND = mapInline(arg.toBufferND()) { ln(it) } @@ -53,6 +50,9 @@ public sealed class ComplexFieldOpsND : BufferedFieldOpsND): BufferND = mapInline(arg.toBufferND()) { acosh(it) } override fun atanh(arg: StructureND): BufferND = mapInline(arg.toBufferND()) { atanh(it) } + override fun power(arg: StructureND, pow: Number): StructureND = + mapInline(arg.toBufferND()) { power(it,pow) } + public companion object : ComplexFieldOpsND() } @@ -63,7 +63,8 @@ public val ComplexField.bufferAlgebra: BufferFieldOps @OptIn(UnstableKMathAPI::class) public class ComplexFieldND(override val shape: Shape) : - ComplexFieldOpsND(), FieldND, NumbersAddOps> { + ComplexFieldOpsND(), FieldND, + NumbersAddOps> { override fun number(value: Number): BufferND { val d = value.toDouble() // minimize conversions diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt index 704c4edd8..aa9dd01ce 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt @@ -252,7 +252,7 @@ public class SimpleAutoDiffExpression>( * Generate [AutoDiffProcessor] for [SimpleAutoDiffExpression] */ public fun > simpleAutoDiff( - field: F + field: F, ): AutoDiffProcessor, SimpleAutoDiffField> = AutoDiffProcessor { function -> SimpleAutoDiffExpression(field, function) @@ -272,8 +272,8 @@ public fun > SimpleAutoDiffField.sqrt(x: Aut public fun > SimpleAutoDiffField.pow( x: AutoDiffValue, y: Double, -): AutoDiffValue = derive(const { power(x.value, y) }) { z -> - x.d += z.d * y * power(x.value, y - 1) +): AutoDiffValue = derive(const { x.value.pow(y)}) { z -> + x.d += z.d * y * x.value.pow(y - 1) } public fun > SimpleAutoDiffField.pow( diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt index 4e52c8ba9..113bd4c52 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt @@ -35,7 +35,7 @@ public interface WithShape { * @param T the type of ND-structure element. * @param C the type of the element context. */ -public interface AlgebraND> { +public interface AlgebraND>: Algebra> { /** * The algebra over elements of ND structure. */ diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt index 0e094a8c7..d25b455f4 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt @@ -53,21 +53,28 @@ public interface BufferAlgebraND> : AlgebraND { public inline fun > BufferAlgebraND.mapInline( arg: BufferND, - crossinline transform: A.(T) -> T + crossinline transform: A.(T) -> T, ): BufferND { val indexes = arg.indices - return BufferND(indexes, bufferAlgebra.mapInline(arg.buffer, transform)) + val buffer = arg.buffer + return BufferND( + indexes, + bufferAlgebra.run { + bufferFactory(buffer.size) { elementAlgebra.transform(buffer[it]) } + } + ) } internal inline fun > BufferAlgebraND.mapIndexedInline( arg: BufferND, - crossinline transform: A.(index: IntArray, arg: T) -> T + crossinline transform: A.(index: IntArray, arg: T) -> T, ): BufferND { val indexes = arg.indices + val buffer = arg.buffer return BufferND( indexes, - bufferAlgebra.mapIndexedInline(arg.buffer) { offset, value -> - transform(indexes.index(offset), value) + bufferAlgebra.run { + bufferFactory(buffer.size) { elementAlgebra.transform(indexes.index(it), buffer[it]) } } ) } @@ -75,35 +82,42 @@ internal inline fun > BufferAlgebraND.mapIndexedInline( internal inline fun > BufferAlgebraND.zipInline( l: BufferND, r: BufferND, - crossinline block: A.(l: T, r: T) -> T + crossinline block: A.(l: T, r: T) -> T, ): BufferND { require(l.indices == r.indices) { "Zip requires the same shapes, but found ${l.shape} on the left and ${r.shape} on the right" } val indexes = l.indices - return BufferND(indexes, bufferAlgebra.zipInline(l.buffer, r.buffer, block)) + val lbuffer = l.buffer + val rbuffer = r.buffer + return BufferND( + indexes, + bufferAlgebra.run { + bufferFactory(lbuffer.size) { elementAlgebra.block(lbuffer[it], rbuffer[it]) } + } + ) } @OptIn(PerformancePitfall::class) public open class BufferedGroupNDOps>( override val bufferAlgebra: BufferAlgebra, - override val indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder + override val indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder, ) : GroupOpsND, BufferAlgebraND { override fun StructureND.unaryMinus(): StructureND = map { -it } } public open class BufferedRingOpsND>( bufferAlgebra: BufferAlgebra, - indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder + indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder, ) : BufferedGroupNDOps(bufferAlgebra, indexerBuilder), RingOpsND public open class BufferedFieldOpsND>( bufferAlgebra: BufferAlgebra, - indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder + indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder, ) : BufferedRingOpsND(bufferAlgebra, indexerBuilder), FieldOpsND { public constructor( elementAlgebra: A, bufferFactory: BufferFactory, - indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder + indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder, ) : this(BufferFieldOps(elementAlgebra, bufferFactory), indexerBuilder) @OptIn(PerformancePitfall::class) @@ -117,11 +131,11 @@ public val > BufferAlgebra.nd: BufferedFieldOpsND ge public fun > BufferAlgebraND.structureND( vararg shape: Int, - initializer: A.(IntArray) -> T + initializer: A.(IntArray) -> T, ): BufferND = structureND(shape, initializer) public fun , A> A.structureND( - initializer: EA.(IntArray) -> T + initializer: EA.(IntArray) -> T, ): BufferND where A : BufferAlgebraND, A : WithShape = structureND(shape, initializer) //// group factories diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt index 19924616d..515988159 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt @@ -5,7 +5,6 @@ package space.kscience.kmath.nd -import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.BufferFactory import space.kscience.kmath.structures.MutableBuffer @@ -27,11 +26,6 @@ public open class BufferND( override val shape: IntArray get() = indices.shape - @PerformancePitfall - override fun elements(): Sequence> = indices.asSequence().map { - it to this[it] - } - override fun toString(): String = StructureND.toString(this) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt index 8cc5472a9..8baeac21f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt @@ -11,7 +11,7 @@ import space.kscience.kmath.operations.* import space.kscience.kmath.structures.DoubleBuffer import kotlin.contracts.InvocationKind import kotlin.contracts.contract -import kotlin.math.pow +import kotlin.math.pow as kpow public class DoubleBufferND( indexes: ShapeIndexer, @@ -30,9 +30,9 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D } } - private inline fun mapInline( + protected inline fun mapInline( arg: DoubleBufferND, - transform: (Double) -> Double + transform: (Double) -> Double, ): DoubleBufferND { val indexes = arg.indices val array = arg.buffer.array @@ -42,7 +42,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D private inline fun zipInline( l: DoubleBufferND, r: DoubleBufferND, - block: (l: Double, r: Double) -> Double + block: (l: Double, r: Double) -> Double, ): DoubleBufferND { require(l.indices == r.indices) { "Zip requires the same shapes, but found ${l.shape} on the left and ${r.shape} on the right" } val indexes = l.indices @@ -60,7 +60,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D override fun zip( left: StructureND, right: StructureND, - transform: DoubleField.(Double, Double) -> Double + transform: DoubleField.(Double, Double) -> Double, ): BufferND = zipInline(left.toBufferND(), right.toBufferND()) { l, r -> DoubleField.transform(l, r) } override fun structureND(shape: Shape, initializer: DoubleField.(IntArray) -> Double): DoubleBufferND { @@ -81,8 +81,8 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D override fun StructureND.unaryMinus(): DoubleBufferND = mapInline(toBufferND()) { -it } - override fun StructureND.div(other: StructureND): DoubleBufferND = - zipInline(toBufferND(), other.toBufferND()) { l, r -> l / r } + override fun StructureND.div(arg: StructureND): DoubleBufferND = + zipInline(toBufferND(), arg.toBufferND()) { l, r -> l / r } override fun divide(left: StructureND, right: StructureND): DoubleBufferND = zipInline(left.toBufferND(), right.toBufferND()) { l: Double, r: Double -> l / r } @@ -101,8 +101,8 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D override fun StructureND.minus(arg: StructureND): DoubleBufferND = zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l - r } - override fun StructureND.times(other: StructureND): DoubleBufferND = - zipInline(toBufferND(), other.toBufferND()) { l: Double, r: Double -> l * r } + override fun StructureND.times(arg: StructureND): DoubleBufferND = + zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l * r } override fun StructureND.times(k: Number): DoubleBufferND = mapInline(toBufferND()) { it * k.toDouble() } @@ -123,9 +123,6 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D override fun scale(a: StructureND, value: Double): DoubleBufferND = mapInline(a.toBufferND()) { it * value } - override fun power(arg: StructureND, pow: Number): DoubleBufferND = - mapInline(arg.toBufferND()) { it.pow(pow.toDouble()) } - override fun exp(arg: StructureND): DoubleBufferND = mapInline(arg.toBufferND()) { kotlin.math.exp(it) } @@ -173,7 +170,38 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND(D @OptIn(UnstableKMathAPI::class) public class DoubleFieldND(override val shape: Shape) : - DoubleFieldOpsND(), FieldND, NumbersAddOps> { + DoubleFieldOpsND(), FieldND, NumbersAddOps>, + ExtendedField> { + + override fun power(arg: StructureND, pow: UInt): DoubleBufferND = mapInline(arg.toBufferND()) { + it.kpow(pow.toInt()) + } + + override fun power(arg: StructureND, pow: Int): DoubleBufferND = mapInline(arg.toBufferND()) { + it.kpow(pow) + } + + override fun power(arg: StructureND, pow: Number): DoubleBufferND = if(pow.isInteger()){ + power(arg, pow.toInt()) + } else { + val dpow = pow.toDouble() + mapInline(arg.toBufferND()) { + if (it < 0) throw IllegalArgumentException("Can't raise negative $it to a fractional power") + else it.kpow(dpow) + } + } + + override fun sinh(arg: StructureND): DoubleBufferND = super.sinh(arg) + + override fun cosh(arg: StructureND): DoubleBufferND = super.cosh(arg) + + override fun tanh(arg: StructureND): DoubleBufferND = super.tan(arg) + + override fun asinh(arg: StructureND): DoubleBufferND = super.asinh(arg) + + override fun acosh(arg: StructureND): DoubleBufferND = super.acosh(arg) + + override fun atanh(arg: StructureND): DoubleBufferND = super.atanh(arg) override fun number(value: Number): DoubleBufferND { val d = value.toDouble() // minimize conversions diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt index c6af75237..244b9fea7 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt @@ -6,6 +6,8 @@ package space.kscience.kmath.operations import space.kscience.kmath.expressions.Symbol +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.Ring.Companion.optimizedPower /** * Stub for DSL the [Algebra] is. @@ -99,6 +101,14 @@ public interface Algebra { */ public fun binaryOperation(operation: String, left: T, right: T): T = binaryOperationFunction(operation)(left, right) + + /** + * Export an algebra element, so it could be accessed even after algebra scope is closed. + * This method must be used on algebras where data is stored externally or any local algebra state is used. + * By default (if not overridden), exports the object itself. + */ + @UnstableKMathAPI + public fun export(arg: T): T = arg } public fun Algebra.bindSymbolOrNull(symbol: Symbol): T? = bindSymbolOrNull(symbol.identity) @@ -162,6 +172,7 @@ public interface GroupOps : Algebra { * @return the difference. */ public operator fun T.minus(arg: T): T = add(this, -arg) + // Dynamic dispatch of operations override fun unaryOperationFunction(operation: String): (arg: T) -> T = when (operation) { PLUS_OPERATION -> { arg -> +arg } @@ -247,6 +258,40 @@ public interface Ring : Group, RingOps { * The neutral element of multiplication */ public val one: T + + /** + * Raises [arg] to the integer power [pow]. + */ + public fun power(arg: T, pow: UInt): T = optimizedPower(arg, pow) + + public companion object{ + /** + * Raises [arg] to the non-negative integer power [exponent]. + * + * Special case: 0 ^ 0 is 1. + * + * @receiver the algebra to provide multiplication. + * @param arg the base. + * @param exponent the exponent. + * @return the base raised to the power. + * @author Evgeniy Zhelenskiy + */ + internal fun Ring.optimizedPower(arg: T, exponent: UInt): T = when { + arg == zero && exponent > 0U -> zero + arg == one -> arg + arg == -one -> powWithoutOptimization(arg, exponent % 2U) + else -> powWithoutOptimization(arg, exponent) + } + + private fun Ring.powWithoutOptimization(base: T, exponent: UInt): T = when (exponent) { + 0U -> one + 1U -> base + else -> { + val pre = powWithoutOptimization(base, exponent shr 1).let { it * it } + if (exponent and 1U == 0U) pre else pre * base + } + } + } } /** @@ -270,10 +315,10 @@ public interface FieldOps : RingOps { * Division of two elements. * * @receiver the dividend. - * @param other the divisor. + * @param arg the divisor. * @return the quotient. */ - public operator fun T.div(other: T): T = divide(this, other) + public operator fun T.div(arg: T): T = divide(this, arg) override fun binaryOperationFunction(operation: String): (left: T, right: T) -> T = when (operation) { DIV_OPERATION -> ::divide @@ -297,4 +342,24 @@ public interface FieldOps : RingOps { */ public interface Field : Ring, FieldOps, ScaleOperations, NumericAlgebra { override fun number(value: Number): T = scale(one, value.toDouble()) + + public fun power(arg: T, pow: Int): T = optimizedPower(arg, pow) + + public companion object{ + /** + * Raises [arg] to the integer power [exponent]. + * + * Special case: 0 ^ 0 is 1. + * + * @receiver the algebra to provide multiplication and division. + * @param arg the base. + * @param exponent the exponent. + * @return the base raised to the power. + * @author Iaroslav Postovalov, Evgeniy Zhelenskiy + */ + private fun Field.optimizedPower(arg: T, exponent: Int): T = when { + exponent < 0 -> one / (this as Ring).optimizedPower(arg, if (exponent == Int.MIN_VALUE) Int.MAX_VALUE.toUInt().inc() else (-exponent).toUInt()) + else -> (this as Ring).optimizedPower(arg, exponent.toUInt()) + } + } } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt index bc05f3904..634a115c7 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.operations +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.BufferFactory import space.kscience.kmath.structures.DoubleBuffer @@ -34,11 +35,13 @@ public interface BufferAlgebra> : Algebra> { public fun Buffer.zip(other: Buffer, block: A.(left: T, right: T) -> T): Buffer = zipInline(this, other, block) + @UnstableKMathAPI override fun unaryOperationFunction(operation: String): (arg: Buffer) -> Buffer { val operationFunction = elementAlgebra.unaryOperationFunction(operation) return { arg -> bufferFactory(arg.size) { operationFunction(arg[it]) } } } + @UnstableKMathAPI override fun binaryOperationFunction(operation: String): (left: Buffer, right: Buffer) -> Buffer { val operationFunction = elementAlgebra.binaryOperationFunction(operation) return { left, right -> @@ -50,7 +53,7 @@ public interface BufferAlgebra> : Algebra> { /** * Inline map */ -public inline fun > BufferAlgebra.mapInline( +private inline fun > BufferAlgebra.mapInline( buffer: Buffer, crossinline block: A.(T) -> T ): Buffer = bufferFactory(buffer.size) { elementAlgebra.block(buffer[it]) } @@ -58,7 +61,7 @@ public inline fun > BufferAlgebra.mapInline( /** * Inline map */ -public inline fun > BufferAlgebra.mapIndexedInline( +private inline fun > BufferAlgebra.mapIndexedInline( buffer: Buffer, crossinline block: A.(index: Int, arg: T) -> T ): Buffer = bufferFactory(buffer.size) { elementAlgebra.block(it, buffer[it]) } @@ -66,7 +69,7 @@ public inline fun > BufferAlgebra.mapIndexedInline( /** * Inline zip */ -public inline fun > BufferAlgebra.zipInline( +private inline fun > BufferAlgebra.zipInline( l: Buffer, r: Buffer, crossinline block: A.(l: T, r: T) -> T @@ -126,7 +129,7 @@ public fun > BufferAlgebra.atanh(arg: Buff mapInline(arg) { atanh(it) } public fun > BufferAlgebra.pow(arg: Buffer, pow: Number): Buffer = - mapInline(arg) { power(it, pow) } + mapInline(arg) {it.pow(pow) } public open class BufferRingOps>( @@ -138,9 +141,11 @@ public open class BufferRingOps>( override fun multiply(left: Buffer, right: Buffer): Buffer = zipInline(left, right) { l, r -> l * r } override fun Buffer.unaryMinus(): Buffer = map { -it } + @UnstableKMathAPI override fun unaryOperationFunction(operation: String): (arg: Buffer) -> Buffer = super.unaryOperationFunction(operation) + @UnstableKMathAPI override fun binaryOperationFunction(operation: String): (left: Buffer, right: Buffer) -> Buffer = super.binaryOperationFunction(operation) } @@ -160,6 +165,7 @@ public open class BufferFieldOps>( override fun scale(a: Buffer, value: Double): Buffer = a.map { scale(it, value) } override fun Buffer.unaryMinus(): Buffer = map { -it } + @UnstableKMathAPI override fun binaryOperationFunction(operation: String): (left: Buffer, right: Buffer) -> Buffer = super.binaryOperationFunction(operation) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt index 060ea5a7e..0deb647a3 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt @@ -5,6 +5,8 @@ package space.kscience.kmath.operations +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.DoubleField.pow import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer @@ -27,7 +29,20 @@ public class DoubleBufferField(public val size: Int) : ExtendedField): DoubleBuffer = super.acosh(arg) - override fun atanh(arg: Buffer): DoubleBuffer= super.atanh(arg) + override fun atanh(arg: Buffer): DoubleBuffer = super.atanh(arg) + + override fun power(arg: Buffer, pow: Number): DoubleBuffer = if (pow.isInteger()) { + arg.mapInline { it.pow(pow.toInt()) } + } else { + arg.mapInline { + if(it<0) throw IllegalArgumentException("Negative argument $it could not be raised to the fractional power") + it.pow(pow.toDouble()) + } + } + + @UnstableKMathAPI + override fun unaryOperationFunction(operation: String): (arg: Buffer) -> Buffer = + super.unaryOperationFunction(operation) // override fun number(value: Number): Buffer = DoubleBuffer(size) { value.toDouble() } // diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOps.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOps.kt index 29b25aae8..b0cce91d3 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOps.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOps.kt @@ -6,21 +6,35 @@ package space.kscience.kmath.operations import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.BufferFactory import space.kscience.kmath.structures.DoubleBuffer +import space.kscience.kmath.structures.asBuffer import kotlin.math.* /** * [ExtendedFieldOps] over [DoubleBuffer]. */ -public abstract class DoubleBufferOps : ExtendedFieldOps>, Norm, Double> { +public abstract class DoubleBufferOps : BufferAlgebra, ExtendedFieldOps>, + Norm, Double> { - override fun Buffer.unaryMinus(): DoubleBuffer = if (this is DoubleBuffer) { - DoubleBuffer(size) { -array[it] } - } else { - DoubleBuffer(size) { -get(it) } - } + override val elementAlgebra: DoubleField get() = DoubleField + override val bufferFactory: BufferFactory get() = ::DoubleBuffer + + override fun Buffer.map(block: DoubleField.(Double) -> Double): DoubleBuffer = + mapInline { DoubleField.block(it) } + + @UnstableKMathAPI + override fun unaryOperationFunction(operation: String): (arg: Buffer) -> Buffer = + super.unaryOperationFunction(operation) + + @UnstableKMathAPI + override fun binaryOperationFunction(operation: String): (left: Buffer, right: Buffer) -> Buffer = + super.binaryOperationFunction(operation) + + override fun Buffer.unaryMinus(): DoubleBuffer = mapInline { -it } override fun add(left: Buffer, right: Buffer): DoubleBuffer { require(right.size == left.size) { @@ -34,18 +48,18 @@ public abstract class DoubleBufferOps : ExtendedFieldOps>, Norm.plus(other: Buffer): DoubleBuffer = add(this, other) + override fun Buffer.plus(arg: Buffer): DoubleBuffer = add(this, arg) - override fun Buffer.minus(other: Buffer): DoubleBuffer { - require(other.size == this.size) { - "The size of the first buffer ${this.size} should be the same as for second one: ${other.size} " + override fun Buffer.minus(arg: Buffer): DoubleBuffer { + require(arg.size == this.size) { + "The size of the first buffer ${this.size} should be the same as for second one: ${arg.size} " } - return if (this is DoubleBuffer && other is DoubleBuffer) { + return if (this is DoubleBuffer && arg is DoubleBuffer) { val aArray = this.array - val bArray = other.array + val bArray = arg.array DoubleBuffer(DoubleArray(this.size) { aArray[it] - bArray[it] }) - } else DoubleBuffer(DoubleArray(this.size) { this[it] - other[it] }) + } else DoubleBuffer(DoubleArray(this.size) { this[it] - arg[it] }) } // @@ -76,8 +90,7 @@ public abstract class DoubleBufferOps : ExtendedFieldOps>, Norm, right: Buffer): DoubleBuffer { @@ -92,101 +105,46 @@ public abstract class DoubleBufferOps : ExtendedFieldOps>, Norm): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { sin(array[it]) }) - } else DoubleBuffer(DoubleArray(arg.size) { sin(arg[it]) }) + override fun sin(arg: Buffer): DoubleBuffer = arg.mapInline(::sin) - override fun cos(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { cos(array[it]) }) - } else DoubleBuffer(DoubleArray(arg.size) { cos(arg[it]) }) + override fun cos(arg: Buffer): DoubleBuffer = arg.mapInline(::cos) - override fun tan(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { tan(array[it]) }) - } else DoubleBuffer(DoubleArray(arg.size) { tan(arg[it]) }) + override fun tan(arg: Buffer): DoubleBuffer = arg.mapInline(::tan) - override fun asin(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { asin(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { asin(arg[it]) }) + override fun asin(arg: Buffer): DoubleBuffer = arg.mapInline(::asin) - override fun acos(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { acos(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { acos(arg[it]) }) + override fun acos(arg: Buffer): DoubleBuffer = arg.mapInline(::acos) - override fun atan(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { atan(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { atan(arg[it]) }) + override fun atan(arg: Buffer): DoubleBuffer = arg.mapInline(::atan) - override fun sinh(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { sinh(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { sinh(arg[it]) }) + override fun sinh(arg: Buffer): DoubleBuffer = arg.mapInline(::sinh) - override fun cosh(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { cosh(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { cosh(arg[it]) }) + override fun cosh(arg: Buffer): DoubleBuffer = arg.mapInline(::cosh) - override fun tanh(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { tanh(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { tanh(arg[it]) }) + override fun tanh(arg: Buffer): DoubleBuffer = arg.mapInline(::tanh) - override fun asinh(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { asinh(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { asinh(arg[it]) }) + override fun asinh(arg: Buffer): DoubleBuffer = arg.mapInline(::asinh) - override fun acosh(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { acosh(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { acosh(arg[it]) }) + override fun acosh(arg: Buffer): DoubleBuffer = arg.mapInline(::acosh) - override fun atanh(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { atanh(array[it]) }) - } else - DoubleBuffer(DoubleArray(arg.size) { atanh(arg[it]) }) + override fun atanh(arg: Buffer): DoubleBuffer = arg.mapInline(::atanh) - override fun power(arg: Buffer, pow: Number): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { array[it].pow(pow.toDouble()) }) - } else - DoubleBuffer(DoubleArray(arg.size) { arg[it].pow(pow.toDouble()) }) + override fun exp(arg: Buffer): DoubleBuffer = arg.mapInline(::exp) - override fun exp(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { exp(array[it]) }) - } else DoubleBuffer(DoubleArray(arg.size) { exp(arg[it]) }) - - override fun ln(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { - val array = arg.array - DoubleBuffer(DoubleArray(arg.size) { ln(array[it]) }) - } else { - DoubleBuffer(DoubleArray(arg.size) { ln(arg[it]) }) - } + override fun ln(arg: Buffer): DoubleBuffer = arg.mapInline(::ln) override fun norm(arg: Buffer): Double = DoubleL2Norm.norm(arg) - override fun scale(a: Buffer, value: Double): DoubleBuffer = if (a is DoubleBuffer) { - val aArray = a.array - DoubleBuffer(DoubleArray(a.size) { aArray[it] * value }) - } else DoubleBuffer(DoubleArray(a.size) { a[it] * value }) + override fun scale(a: Buffer, value: Double): DoubleBuffer = a.mapInline { it * value } - public companion object : DoubleBufferOps() + public companion object : DoubleBufferOps() { + public inline fun Buffer.mapInline(block: (Double) -> Double): DoubleBuffer = + if (this is DoubleBuffer) { + DoubleArray(size) { block(array[it]) }.asBuffer() + } else { + DoubleArray(size) { block(get(it)) }.asBuffer() + } + } } public object DoubleL2Norm : Norm, Double> { diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/OptionalOperations.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/OptionalOperations.kt index d32e03533..332617158 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/OptionalOperations.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/OptionalOperations.kt @@ -74,14 +74,21 @@ public interface TrigonometricOperations : Algebra { } } +/** + * Check if number is an integer from platform point of view + */ +public expect fun Number.isInteger(): Boolean + /** * A context extension to include power operations based on exponentiation. * * @param T the type of element of this structure. */ -public interface PowerOperations : Algebra { +public interface PowerOperations : FieldOps { + /** - * Raises [arg] to the power [pow]. + * Raises [arg] to a power if possible (negative number could not be raised to a fractional power). + * Throws [IllegalArgumentException] if not possible. */ public fun power(arg: T, pow: Number): T @@ -108,6 +115,7 @@ public interface PowerOperations : Algebra { } } + /** * A container for operations related to `exp` and `ln` functions. * diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt index b26ebb2ea..493d90d2f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt @@ -96,46 +96,3 @@ public fun Iterable.averageWith(space: S): T where S : Ring, S : Sc public fun Sequence.averageWith(space: S): T where S : Ring, S : ScaleOperations = space.average(this) -/** - * Raises [arg] to the non-negative integer power [exponent]. - * - * Special case: 0 ^ 0 is 1. - * - * @receiver the algebra to provide multiplication. - * @param arg the base. - * @param exponent the exponent. - * @return the base raised to the power. - * @author Evgeniy Zhelenskiy - */ -public fun Ring.power(arg: T, exponent: UInt): T = when { - arg == zero && exponent > 0U -> zero - arg == one -> arg - arg == -one -> powWithoutOptimization(arg, exponent % 2U) - else -> powWithoutOptimization(arg, exponent) -} - -private fun Ring.powWithoutOptimization(base: T, exponent: UInt): T = when (exponent) { - 0U -> one - 1U -> base - else -> { - val pre = powWithoutOptimization(base, exponent shr 1).let { it * it } - if (exponent and 1U == 0U) pre else pre * base - } -} - - -/** - * Raises [arg] to the integer power [exponent]. - * - * Special case: 0 ^ 0 is 1. - * - * @receiver the algebra to provide multiplication and division. - * @param arg the base. - * @param exponent the exponent. - * @return the base raised to the power. - * @author Iaroslav Postovalov, Evgeniy Zhelenskiy - */ -public fun Field.power(arg: T, exponent: Int): T = when { - exponent < 0 -> one / (this as Ring).power(arg, if (exponent == Int.MIN_VALUE) Int.MAX_VALUE.toUInt().inc() else (-exponent).toUInt()) - else -> (this as Ring).power(arg, exponent.toUInt()) -} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt index 0b111349c..7c8030168 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt @@ -13,9 +13,8 @@ import kotlin.math.pow as kpow public interface ExtendedFieldOps : FieldOps, TrigonometricOperations, - PowerOperations, ExponentialOperations, - ScaleOperations { + ScaleOperations { override fun tan(arg: T): T = sin(arg) / cos(arg) override fun tanh(arg: T): T = sinh(arg) / cosh(arg) @@ -26,7 +25,6 @@ public interface ExtendedFieldOps : TrigonometricOperations.ACOS_OPERATION -> ::acos TrigonometricOperations.ASIN_OPERATION -> ::asin TrigonometricOperations.ATAN_OPERATION -> ::atan - PowerOperations.SQRT_OPERATION -> ::sqrt ExponentialOperations.EXP_OPERATION -> ::exp ExponentialOperations.LN_OPERATION -> ::ln ExponentialOperations.COSH_OPERATION -> ::cosh @@ -42,7 +40,7 @@ public interface ExtendedFieldOps : /** * Advanced Number-like field that implements basic operations. */ -public interface ExtendedField : ExtendedFieldOps, Field, NumericAlgebra{ +public interface ExtendedField : ExtendedFieldOps, Field, PowerOperations, NumericAlgebra { override fun sinh(arg: T): T = (exp(arg) - exp(-arg)) / 2.0 override fun cosh(arg: T): T = (exp(arg) + exp(-arg)) / 2.0 override fun tanh(arg: T): T = (exp(arg) - exp(-arg)) / (exp(-arg) + exp(arg)) @@ -50,6 +48,11 @@ public interface ExtendedField : ExtendedFieldOps, Field, NumericAlgebr override fun acosh(arg: T): T = ln(arg + sqrt((arg - one) * (arg + one))) override fun atanh(arg: T): T = (ln(arg + one) - ln(one - arg)) / 2.0 + override fun unaryOperationFunction(operation: String): (arg: T) -> T { + return if (operation == PowerOperations.SQRT_OPERATION) ::sqrt + else super.unaryOperationFunction(operation) + } + override fun rightSideNumberOperationFunction(operation: String): (left: T, right: Number) -> T = when (operation) { PowerOperations.POW_OPERATION -> ::power @@ -69,7 +72,7 @@ public object DoubleField : ExtendedField, Norm, ScaleOp override fun binaryOperationFunction(operation: String): (left: Double, right: Double) -> Double = when (operation) { - PowerOperations.POW_OPERATION -> ::power + PowerOperations.POW_OPERATION -> { l, r -> l.kpow(r) } else -> super.binaryOperationFunction(operation) } @@ -94,18 +97,23 @@ public object DoubleField : ExtendedField, Norm, ScaleOp override inline fun acosh(arg: Double): Double = kotlin.math.acosh(arg) override inline fun atanh(arg: Double): Double = kotlin.math.atanh(arg) - override inline fun sqrt(arg: Double): Double = kotlin.math.sqrt(arg) - override inline fun power(arg: Double, pow: Number): Double = arg.kpow(pow.toDouble()) + override fun sqrt(arg: Double): Double = kotlin.math.sqrt(arg) + override fun power(arg: Double, pow: Number): Double = when { + pow.isInteger() -> arg.kpow(pow.toInt()) + arg < 0 -> throw IllegalArgumentException("Can't raise negative $arg to a fractional power $pow") + else -> arg.kpow(pow.toDouble()) + } + override inline fun exp(arg: Double): Double = kotlin.math.exp(arg) override inline fun ln(arg: Double): Double = kotlin.math.ln(arg) override inline fun norm(arg: Double): Double = abs(arg) override inline fun Double.unaryMinus(): Double = -this - override inline fun Double.plus(other: Double): Double = this + other - override inline fun Double.minus(other: Double): Double = this - other - override inline fun Double.times(other: Double): Double = this * other - override inline fun Double.div(other: Double): Double = this / other + override inline fun Double.plus(arg: Double): Double = this + arg + override inline fun Double.minus(arg: Double): Double = this - arg + override inline fun Double.times(arg: Double): Double = this * arg + override inline fun Double.div(arg: Double): Double = this / arg } public val Double.Companion.algebra: DoubleField get() = DoubleField @@ -122,7 +130,7 @@ public object FloatField : ExtendedField, Norm { override fun binaryOperationFunction(operation: String): (left: Float, right: Float) -> Float = when (operation) { - PowerOperations.POW_OPERATION -> ::power + PowerOperations.POW_OPERATION -> { l, r -> l.kpow(r) } else -> super.binaryOperationFunction(operation) } @@ -149,16 +157,17 @@ public object FloatField : ExtendedField, Norm { override inline fun sqrt(arg: Float): Float = kotlin.math.sqrt(arg) override inline fun power(arg: Float, pow: Number): Float = arg.kpow(pow.toFloat()) + override inline fun exp(arg: Float): Float = kotlin.math.exp(arg) override inline fun ln(arg: Float): Float = kotlin.math.ln(arg) override inline fun norm(arg: Float): Float = abs(arg) override inline fun Float.unaryMinus(): Float = -this - override inline fun Float.plus(other: Float): Float = this + other - override inline fun Float.minus(other: Float): Float = this - other - override inline fun Float.times(other: Float): Float = this * other - override inline fun Float.div(other: Float): Float = this / other + override inline fun Float.plus(arg: Float): Float = this + arg + override inline fun Float.minus(arg: Float): Float = this - arg + override inline fun Float.times(arg: Float): Float = this * arg + override inline fun Float.div(arg: Float): Float = this / arg } public val Float.Companion.algebra: FloatField get() = FloatField @@ -180,9 +189,9 @@ public object IntRing : Ring, Norm, NumericAlgebra { override inline fun norm(arg: Int): Int = abs(arg) override inline fun Int.unaryMinus(): Int = -this - override inline fun Int.plus(other: Int): Int = this + other - override inline fun Int.minus(other: Int): Int = this - other - override inline fun Int.times(other: Int): Int = this * other + override inline fun Int.plus(arg: Int): Int = this + arg + override inline fun Int.minus(arg: Int): Int = this - arg + override inline fun Int.times(arg: Int): Int = this * arg } public val Int.Companion.algebra: IntRing get() = IntRing @@ -204,9 +213,9 @@ public object ShortRing : Ring, Norm, NumericAlgebra override fun norm(arg: Short): Short = if (arg > 0) arg else (-arg).toShort() override inline fun Short.unaryMinus(): Short = (-this).toShort() - override inline fun Short.plus(other: Short): Short = (this + other).toShort() - override inline fun Short.minus(other: Short): Short = (this - other).toShort() - override inline fun Short.times(other: Short): Short = (this * other).toShort() + override inline fun Short.plus(arg: Short): Short = (this + arg).toShort() + override inline fun Short.minus(arg: Short): Short = (this - arg).toShort() + override inline fun Short.times(arg: Short): Short = (this * arg).toShort() } public val Short.Companion.algebra: ShortRing get() = ShortRing @@ -230,7 +239,7 @@ public object ByteRing : Ring, Norm, NumericAlgebra { override inline fun Byte.unaryMinus(): Byte = (-this).toByte() override inline fun Byte.plus(arg: Byte): Byte = (this + arg).toByte() override inline fun Byte.minus(arg: Byte): Byte = (this - arg).toByte() - override inline fun Byte.times(other: Byte): Byte = (this * other).toByte() + override inline fun Byte.times(arg: Byte): Byte = (this * arg).toByte() } public val Byte.Companion.algebra: ByteRing get() = ByteRing @@ -252,9 +261,9 @@ public object LongRing : Ring, Norm, NumericAlgebra { override fun norm(arg: Long): Long = abs(arg) override inline fun Long.unaryMinus(): Long = (-this) - override inline fun Long.plus(other: Long): Long = (this + other) - override inline fun Long.minus(other: Long): Long = (this - other) - override inline fun Long.times(other: Long): Long = (this * other) + override inline fun Long.plus(arg: Long): Long = (this + arg) + override inline fun Long.minus(arg: Long): Long = (this - arg) + override inline fun Long.times(arg: Long): Long = (this * arg) } public val Long.Companion.algebra: LongRing get() = LongRing diff --git a/kmath-core/src/jsMain/kotlin/space/kscience/kmath/operations/isInteger.kt b/kmath-core/src/jsMain/kotlin/space/kscience/kmath/operations/isInteger.kt new file mode 100644 index 000000000..c15669145 --- /dev/null +++ b/kmath-core/src/jsMain/kotlin/space/kscience/kmath/operations/isInteger.kt @@ -0,0 +1,6 @@ +package space.kscience.kmath.operations + +/** + * Check if number is an integer + */ +public actual fun Number.isInteger(): Boolean = js("Number").isInteger(this) as Boolean \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/operations/isInteger.kt b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/operations/isInteger.kt new file mode 100644 index 000000000..746d1e530 --- /dev/null +++ b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/operations/isInteger.kt @@ -0,0 +1,6 @@ +package space.kscience.kmath.operations + +/** + * Check if number is an integer + */ +public actual fun Number.isInteger(): Boolean = (this is Int) || (this is Long) || (this is Short) || (this.toDouble() % 1 == 0.0) \ No newline at end of file diff --git a/kmath-core/src/nativeMain/kotlin/space/kscience/kmath/operations/isInteger.kt b/kmath-core/src/nativeMain/kotlin/space/kscience/kmath/operations/isInteger.kt new file mode 100644 index 000000000..746d1e530 --- /dev/null +++ b/kmath-core/src/nativeMain/kotlin/space/kscience/kmath/operations/isInteger.kt @@ -0,0 +1,6 @@ +package space.kscience.kmath.operations + +/** + * Check if number is an integer + */ +public actual fun Number.isInteger(): Boolean = (this is Int) || (this is Long) || (this is Short) || (this.toDouble() % 1 == 0.0) \ No newline at end of file diff --git a/kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/KMathJafama.kt b/kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/KMathJafama.kt index 9ff7ffc9c..64a935705 100644 --- a/kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/KMathJafama.kt +++ b/kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/KMathJafama.kt @@ -59,8 +59,8 @@ public object JafamaDoubleField : ExtendedField, Norm, S override inline fun Double.unaryMinus(): Double = -this override inline fun Double.plus(arg: Double): Double = this + arg override inline fun Double.minus(arg: Double): Double = this - arg - override inline fun Double.times(other: Double): Double = this * other - override inline fun Double.div(other: Double): Double = this / other + override inline fun Double.times(arg: Double): Double = this * arg + override inline fun Double.div(arg: Double): Double = this / arg } /** @@ -108,8 +108,8 @@ public object StrictJafamaDoubleField : ExtendedField, Norm(), + TrigonometricOperations>, ExponentialOperations> { + override val elementAlgebra: DoubleField get() = DoubleField + override val type: DataType get() = DataType.DoubleDataType + + override fun sin(arg: StructureND): MultikTensor = multikMath.mathEx.sin(arg.asMultik().array).wrap() + + override fun cos(arg: StructureND): MultikTensor = multikMath.mathEx.cos(arg.asMultik().array).wrap() + + override fun tan(arg: StructureND): MultikTensor = sin(arg) / cos(arg) + + override fun asin(arg: StructureND): MultikTensor = arg.map { asin(it) } + + override fun acos(arg: StructureND): MultikTensor = arg.map { acos(it) } + + override fun atan(arg: StructureND): MultikTensor = arg.map { atan(it) } + + override fun exp(arg: StructureND): MultikTensor = multikMath.mathEx.exp(arg.asMultik().array).wrap() + + override fun ln(arg: StructureND): MultikTensor = multikMath.mathEx.log(arg.asMultik().array).wrap() + + override fun sinh(arg: StructureND): MultikTensor = (exp(arg) - exp(-arg)) / 2.0 + + override fun cosh(arg: StructureND): MultikTensor = (exp(arg) + exp(-arg)) / 2.0 + + override fun tanh(arg: StructureND): MultikTensor { + val expPlus = exp(arg) + val expMinus = exp(-arg) + return (expPlus - expMinus) / (expPlus + expMinus) + } + + override fun asinh(arg: StructureND): MultikTensor = arg.map { asinh(it) } + + override fun acosh(arg: StructureND): MultikTensor = arg.map { acosh(it) } + + override fun atanh(arg: StructureND): MultikTensor = arg.map { atanh(it) } +} + +public val Double.Companion.multikAlgebra: MultikTensorAlgebra get() = MultikDoubleAlgebra +public val DoubleField.multikAlgebra: MultikTensorAlgebra get() = MultikDoubleAlgebra + diff --git a/kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikTensorAlgebra.kt b/kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikTensorAlgebra.kt index 02959d4fb..0c3cf6b14 100644 --- a/kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikTensorAlgebra.kt +++ b/kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikTensorAlgebra.kt @@ -7,11 +7,9 @@ package space.kscience.kmath.multik -import org.jetbrains.kotlinx.multik.api.Multik -import org.jetbrains.kotlinx.multik.api.linalg.dot -import org.jetbrains.kotlinx.multik.api.mk -import org.jetbrains.kotlinx.multik.api.ndarrayOf -import org.jetbrains.kotlinx.multik.api.zeros +import org.jetbrains.kotlinx.multik.api.* +import org.jetbrains.kotlinx.multik.api.linalg.LinAlg +import org.jetbrains.kotlinx.multik.api.math.Math import org.jetbrains.kotlinx.multik.ndarray.data.* import org.jetbrains.kotlinx.multik.ndarray.operations.* import space.kscience.kmath.misc.PerformancePitfall @@ -52,10 +50,16 @@ private fun MultiArray.asD2Array(): D2Array { else throw ClassCastException("Cannot cast MultiArray to NDArray.") } -public abstract class MultikTensorAlgebra> : TensorAlgebra where T : Number, T : Comparable { +public abstract class MultikTensorAlgebra> : TensorAlgebra + where T : Number, T : Comparable { public abstract val type: DataType + protected val multikMath: Math = mk.math + protected val multikLinAl: LinAlg = mk.linalg + protected val multikStat: Statistics = mk.stat + + override fun structureND(shape: Shape, initializer: A.(IntArray) -> T): MultikTensor { val strides = DefaultStrides(shape) val memoryView = initMemoryView(strides.linearSize, type) @@ -65,6 +69,7 @@ public abstract class MultikTensorAlgebra> : TensorAlgebra return MultikTensor(NDArray(memoryView, shape = shape, dim = DN(shape.size))) } + @OptIn(PerformancePitfall::class) override fun StructureND.map(transform: A.(T) -> T): MultikTensor = if (this is MultikTensor) { val data = initMemoryView(array.size, type) var count = 0 @@ -76,6 +81,7 @@ public abstract class MultikTensorAlgebra> : TensorAlgebra } } + @OptIn(PerformancePitfall::class) override fun StructureND.mapIndexed(transform: A.(index: IntArray, T) -> T): MultikTensor = if (this is MultikTensor) { val array = asMultik().array @@ -96,6 +102,7 @@ public abstract class MultikTensorAlgebra> : TensorAlgebra } } + @OptIn(PerformancePitfall::class) override fun zip(left: StructureND, right: StructureND, transform: A.(T, T) -> T): MultikTensor { require(left.shape.contentEquals(right.shape)) { "ND array shape mismatch" } //TODO replace by ShapeMismatchException val leftArray = left.asMultik().array @@ -208,9 +215,9 @@ public abstract class MultikTensorAlgebra> : TensorAlgebra override fun StructureND.unaryMinus(): MultikTensor = asMultik().array.unaryMinus().wrap() - override fun StructureND.get(i: Int): MultikTensor = asMultik().array.mutableView(i).wrap() + override fun Tensor.get(i: Int): MultikTensor = asMultik().array.mutableView(i).wrap() - override fun StructureND.transpose(i: Int, j: Int): MultikTensor = asMultik().array.transpose(i, j).wrap() + override fun Tensor.transpose(i: Int, j: Int): MultikTensor = asMultik().array.transpose(i, j).wrap() override fun Tensor.view(shape: IntArray): MultikTensor { require(shape.all { it > 0 }) @@ -236,12 +243,12 @@ public abstract class MultikTensorAlgebra> : TensorAlgebra override fun StructureND.dot(other: StructureND): MultikTensor = if (this.shape.size == 1 && other.shape.size == 1) { Multik.ndarrayOf( - asMultik().array.asD1Array() dot other.asMultik().array.asD1Array() - ).asDNArray().wrap() + multikLinAl.linAlgEx.dotVV(asMultik().array.asD1Array(), other.asMultik().array.asD1Array()) + ).wrap() } else if (this.shape.size == 2 && other.shape.size == 2) { - (asMultik().array.asD2Array() dot other.asMultik().array.asD2Array()).asDNArray().wrap() + multikLinAl.linAlgEx.dotMM(asMultik().array.asD2Array(), other.asMultik().array.asD2Array()).wrap() } else if (this.shape.size == 2 && other.shape.size == 1) { - (asMultik().array.asD2Array() dot other.asMultik().array.asD1Array()).asDNArray().wrap() + multikLinAl.linAlgEx.dotMV(asMultik().array.asD2Array(), other.asMultik().array.asD1Array()).wrap() } else { TODO("Not implemented for broadcasting") } @@ -270,7 +277,7 @@ public abstract class MultikTensorAlgebra> : TensorAlgebra TODO("Not yet implemented") } - override fun StructureND.argMax(dim: Int, keepDim: Boolean): Tensor { + override fun StructureND.argMax(dim: Int, keepDim: Boolean): Tensor { TODO("Not yet implemented") } } @@ -294,23 +301,15 @@ public abstract class MultikDivisionTensorAlgebra> } } - override fun Tensor.divAssign(other: StructureND) { + override fun Tensor.divAssign(arg: StructureND) { if (this is MultikTensor) { - array.divAssign(other.asMultik().array) + array.divAssign(arg.asMultik().array) } else { - mapInPlace { index, t -> elementAlgebra.divide(t, other[index]) } + mapInPlace { index, t -> elementAlgebra.divide(t, arg[index]) } } } } -public object MultikDoubleAlgebra : MultikDivisionTensorAlgebra() { - override val elementAlgebra: DoubleField get() = DoubleField - override val type: DataType get() = DataType.DoubleDataType -} - -public val Double.Companion.multikAlgebra: MultikTensorAlgebra get() = MultikDoubleAlgebra -public val DoubleField.multikAlgebra: MultikTensorAlgebra get() = MultikDoubleAlgebra - public object MultikFloatAlgebra : MultikDivisionTensorAlgebra() { override val elementAlgebra: FloatField get() = FloatField override val type: DataType get() = DataType.FloatDataType diff --git a/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayAlgebra.kt b/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayAlgebra.kt index b1cc1f834..dd27bc817 100644 --- a/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayAlgebra.kt +++ b/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayAlgebra.kt @@ -155,7 +155,7 @@ public sealed interface Nd4jArrayField> : FieldOpsND, * Represents intersection of [ExtendedField] and [Field] over [Nd4jArrayStructure]. */ public sealed interface Nd4jArrayExtendedFieldOps> : - ExtendedFieldOps>, Nd4jArrayField { + ExtendedFieldOps>, Nd4jArrayField, PowerOperations> { override fun sin(arg: StructureND): StructureND = Transforms.sin(arg.ndArray).wrap() override fun cos(arg: StructureND): StructureND = Transforms.cos(arg.ndArray).wrap() diff --git a/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jTensorAlgebra.kt b/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jTensorAlgebra.kt index 33765b40a..d7dd6e71b 100644 --- a/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jTensorAlgebra.kt +++ b/kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jTensorAlgebra.kt @@ -92,8 +92,8 @@ public sealed interface Nd4jTensorAlgebra> : AnalyticTe } override fun StructureND.unaryMinus(): Nd4jArrayStructure = ndArray.neg().wrap() - override fun StructureND.get(i: Int): Nd4jArrayStructure = ndArray.slice(i.toLong()).wrap() - override fun StructureND.transpose(i: Int, j: Int): Nd4jArrayStructure = ndArray.swapAxes(i, j).wrap() + override fun Tensor.get(i: Int): Nd4jArrayStructure = ndArray.slice(i.toLong()).wrap() + override fun Tensor.transpose(i: Int, j: Int): Nd4jArrayStructure = ndArray.swapAxes(i, j).wrap() override fun StructureND.dot(other: StructureND): Nd4jArrayStructure = ndArray.mmul(other.ndArray).wrap() override fun StructureND.min(dim: Int, keepDim: Boolean): Nd4jArrayStructure = @@ -108,8 +108,8 @@ public sealed interface Nd4jTensorAlgebra> : AnalyticTe override fun Tensor.view(shape: IntArray): Nd4jArrayStructure = ndArray.reshape(shape).wrap() override fun Tensor.viewAs(other: StructureND): Nd4jArrayStructure = view(other.shape) - override fun StructureND.argMax(dim: Int, keepDim: Boolean): Nd4jArrayStructure = - ndBase.get().argmax(ndArray, keepDim, dim).wrap() + override fun StructureND.argMax(dim: Int, keepDim: Boolean): Tensor = + ndBase.get().argmax(ndArray, keepDim, dim).asIntStructure() override fun StructureND.mean(dim: Int, keepDim: Boolean): Nd4jArrayStructure = ndArray.mean(keepDim, dim).wrap() @@ -148,8 +148,8 @@ public sealed interface Nd4jTensorAlgebra> : AnalyticTe ndArray.divi(value) } - override fun Tensor.divAssign(other: StructureND) { - ndArray.divi(other.ndArray) + override fun Tensor.divAssign(arg: StructureND) { + ndArray.divi(arg.ndArray) } override fun StructureND.variance(dim: Int, keepDim: Boolean): Nd4jArrayStructure = diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/AnalyticTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/AnalyticTensorAlgebra.kt index debfb3ef0..743105fdf 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/AnalyticTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/AnalyticTensorAlgebra.kt @@ -6,6 +6,7 @@ package space.kscience.kmath.tensors.api import space.kscience.kmath.nd.StructureND +import space.kscience.kmath.operations.ExtendedFieldOps import space.kscience.kmath.operations.Field @@ -14,7 +15,8 @@ import space.kscience.kmath.operations.Field * * @param T the type of items closed under analytic functions in the tensors. */ -public interface AnalyticTensorAlgebra> : TensorPartialDivisionAlgebra { +public interface AnalyticTensorAlgebra> : + TensorPartialDivisionAlgebra, ExtendedFieldOps> { /** * @return the mean of all elements in the input tensor. @@ -121,4 +123,27 @@ public interface AnalyticTensorAlgebra> : TensorPartialDivisionA //For information: https://pytorch.org/docs/stable/generated/torch.floor.html#torch.floor public fun StructureND.floor(): Tensor + override fun sin(arg: StructureND): StructureND = arg.sin() + + override fun cos(arg: StructureND): StructureND = arg.cos() + + override fun asin(arg: StructureND): StructureND = arg.asin() + + override fun acos(arg: StructureND): StructureND = arg.acos() + + override fun atan(arg: StructureND): StructureND = arg.atan() + + override fun exp(arg: StructureND): StructureND = arg.exp() + + override fun ln(arg: StructureND): StructureND = arg.ln() + + override fun sinh(arg: StructureND): StructureND = arg.sinh() + + override fun cosh(arg: StructureND): StructureND = arg.cosh() + + override fun asinh(arg: StructureND): StructureND = arg.asinh() + + override fun acosh(arg: StructureND): StructureND = arg.acosh() + + override fun atanh(arg: StructureND): StructureND = arg.atanh() } \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorAlgebra.kt index 8c445cf2d..60fc470fb 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorAlgebra.kt @@ -166,7 +166,7 @@ public interface TensorAlgebra> : RingOpsND { * @param i index of the extractable tensor * @return subtensor of the original tensor with index [i] */ - public operator fun StructureND.get(i: Int): Tensor + public operator fun Tensor.get(i: Int): Tensor /** * Returns a tensor that is a transposed version of this tensor. The given dimensions [i] and [j] are swapped. @@ -176,7 +176,7 @@ public interface TensorAlgebra> : RingOpsND { * @param j the second dimension to be transposed * @return transposed tensor */ - public fun StructureND.transpose(i: Int = -2, j: Int = -1): Tensor + public fun Tensor.transpose(i: Int = -2, j: Int = -1): Tensor /** * Returns a new tensor with the same data as the self tensor but of a different shape. @@ -324,7 +324,7 @@ public interface TensorAlgebra> : RingOpsND { * @param keepDim whether the output tensor has [dim] retained or not. * @return the index of maximum value of each row of the input tensor in the given dimension [dim]. */ - public fun StructureND.argMax(dim: Int, keepDim: Boolean): Tensor + public fun StructureND.argMax(dim: Int, keepDim: Boolean): Tensor override fun add(left: StructureND, right: StructureND): Tensor = left + right diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorPartialDivisionAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorPartialDivisionAlgebra.kt index 304cc66cc..0a1e09081 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorPartialDivisionAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorPartialDivisionAlgebra.kt @@ -53,9 +53,9 @@ public interface TensorPartialDivisionAlgebra> : TensorAlgebra.divAssign(value: T) /** - * Each element of this tensor is divided by each element of the [other] tensor. + * Each element of this tensor is divided by each element of the [arg] tensor. * - * @param other tensor to be divided by. + * @param arg tensor to be divided by. */ - public operator fun Tensor.divAssign(other: StructureND) + public operator fun Tensor.divAssign(arg: StructureND) } diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt index 83d1f7b35..10c747777 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.tensors.core +import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.StructureND import space.kscience.kmath.tensors.api.Tensor @@ -17,6 +18,8 @@ import space.kscience.kmath.tensors.core.internal.tensor * Basic linear algebra operations implemented with broadcasting. * For more information: https://pytorch.org/docs/stable/notes/broadcasting.html */ + +@PerformancePitfall public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { override fun StructureND.plus(arg: StructureND): DoubleTensor { @@ -74,8 +77,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { } } - override fun StructureND.div(other: StructureND): DoubleTensor { - val broadcast = broadcastTensors(tensor, other.tensor) + override fun StructureND.div(arg: StructureND): DoubleTensor { + val broadcast = broadcastTensors(tensor, arg.tensor) val newThis = broadcast[0] val newOther = broadcast[1] val resBuffer = DoubleArray(newThis.indices.linearSize) { i -> @@ -85,8 +88,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { return DoubleTensor(newThis.shape, resBuffer) } - override fun Tensor.divAssign(other: StructureND) { - val newOther = broadcastTo(other.tensor, tensor.shape) + override fun Tensor.divAssign(arg: StructureND) { + val newOther = broadcastTo(arg.tensor, tensor.shape) for (i in 0 until tensor.indices.linearSize) { tensor.mutableBuffer.array()[tensor.bufferStart + i] /= newOther.mutableBuffer.array()[tensor.bufferStart + i] @@ -99,5 +102,6 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() { * Compute a value using broadcast double tensor algebra */ @UnstableKMathAPI +@PerformancePitfall public fun DoubleTensorAlgebra.withBroadcast(block: BroadcastDoubleTensorAlgebra.() -> R): R = BroadcastDoubleTensorAlgebra.block() \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt index eab4d3d5b..54d8f54dc 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt @@ -9,7 +9,6 @@ import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.Strides import space.kscience.kmath.structures.MutableBuffer import space.kscience.kmath.tensors.api.Tensor -import space.kscience.kmath.tensors.core.internal.TensorLinearStructure /** * Represents [Tensor] over a [MutableBuffer] intended to be used through [DoubleTensor] and [IntTensor] diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt index 0856b5c54..def0b43ba 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt @@ -3,13 +3,18 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. */ + +@file:OptIn(PerformancePitfall::class) + package space.kscience.kmath.tensors.core +import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.nd.StructureND import space.kscience.kmath.nd.as1D import space.kscience.kmath.nd.as2D import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.structures.MutableBuffer import space.kscience.kmath.structures.indices import space.kscience.kmath.tensors.api.AnalyticTensorAlgebra import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra @@ -38,6 +43,7 @@ public open class DoubleTensorAlgebra : * @param transform the function to be applied to each element of the tensor. * @return the resulting tensor after applying the function. */ + @PerformancePitfall @Suppress("OVERRIDE_BY_INLINE") final override inline fun StructureND.map(transform: DoubleField.(Double) -> Double): DoubleTensor { val tensor = this.tensor @@ -51,6 +57,7 @@ public open class DoubleTensorAlgebra : ) } + @PerformancePitfall @Suppress("OVERRIDE_BY_INLINE") final override inline fun StructureND.mapIndexed(transform: DoubleField.(index: IntArray, Double) -> Double): DoubleTensor { val tensor = this.tensor @@ -64,12 +71,13 @@ public open class DoubleTensorAlgebra : ) } + @PerformancePitfall override fun zip( left: StructureND, right: StructureND, - transform: DoubleField.(Double, Double) -> Double + transform: DoubleField.(Double, Double) -> Double, ): DoubleTensor { - require(left.shape.contentEquals(right.shape)){ + require(left.shape.contentEquals(right.shape)) { "The shapes in zip are not equal: left - ${left.shape}, right - ${right.shape}" } val leftTensor = left.tensor @@ -115,7 +123,7 @@ public open class DoubleTensorAlgebra : TensorLinearStructure(shape).asSequence().map { DoubleField.initializer(it) }.toMutableList().toDoubleArray() ) - override operator fun StructureND.get(i: Int): DoubleTensor { + override operator fun Tensor.get(i: Int): DoubleTensor { val lastShape = tensor.shape.drop(1).toIntArray() val newShape = if (lastShape.isNotEmpty()) lastShape else intArrayOf(1) val newStart = newShape.reduce(Int::times) * i + tensor.bufferStart @@ -314,11 +322,11 @@ public open class DoubleTensorAlgebra : return DoubleTensor(shape, resBuffer) } - override fun StructureND.div(other: StructureND): DoubleTensor { - checkShapesCompatible(tensor, other) + override fun StructureND.div(arg: StructureND): DoubleTensor { + checkShapesCompatible(tensor, arg) val resBuffer = DoubleArray(tensor.numElements) { i -> - tensor.mutableBuffer.array()[other.tensor.bufferStart + i] / - other.tensor.mutableBuffer.array()[other.tensor.bufferStart + i] + tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] / + arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] } return DoubleTensor(tensor.shape, resBuffer) } @@ -329,11 +337,11 @@ public open class DoubleTensorAlgebra : } } - override fun Tensor.divAssign(other: StructureND) { - checkShapesCompatible(tensor, other) + override fun Tensor.divAssign(arg: StructureND) { + checkShapesCompatible(tensor, arg) for (i in 0 until tensor.numElements) { tensor.mutableBuffer.array()[tensor.bufferStart + i] /= - other.tensor.mutableBuffer.array()[tensor.bufferStart + i] + arg.tensor.mutableBuffer.array()[tensor.bufferStart + i] } } @@ -344,7 +352,7 @@ public open class DoubleTensorAlgebra : return DoubleTensor(tensor.shape, resBuffer) } - override fun StructureND.transpose(i: Int, j: Int): DoubleTensor { + override fun Tensor.transpose(i: Int, j: Int): DoubleTensor { val ii = tensor.minusIndex(i) val jj = tensor.minusIndex(j) checkTranspose(tensor.dimension, ii, jj) @@ -376,6 +384,7 @@ public open class DoubleTensorAlgebra : override fun Tensor.viewAs(other: StructureND): DoubleTensor = tensor.view(other.shape) + @PerformancePitfall override infix fun StructureND.dot(other: StructureND): DoubleTensor { if (tensor.shape.size == 1 && other.shape.size == 1) { return DoubleTensor(intArrayOf(1), doubleArrayOf(tensor.times(other).tensor.mutableBuffer.array().sum())) @@ -413,14 +422,11 @@ public open class DoubleTensorAlgebra : for ((res, ab) in resTensor.matrixSequence().zip(newThis.matrixSequence().zip(newOther.matrixSequence()))) { val (a, b) = ab - dotHelper(a.as2D(), b.as2D(), res.as2D(), l, m1, n) + dotTo(a.as2D(), b.as2D(), res.as2D(), l, m1, n) } return if (penultimateDim) { - resTensor.view( - resTensor.shape.dropLast(2).toIntArray() + - intArrayOf(resTensor.shape.last()) - ) + resTensor.view(resTensor.shape.dropLast(2).toIntArray() + intArrayOf(resTensor.shape.last())) } else if (lastDim) { resTensor.view(resTensor.shape.dropLast(1).toIntArray()) } else { @@ -432,7 +438,7 @@ public open class DoubleTensorAlgebra : diagonalEntries: Tensor, offset: Int, dim1: Int, - dim2: Int + dim2: Int, ): DoubleTensor { val n = diagonalEntries.shape.size val d1 = minusIndexFrom(n + 1, dim1) @@ -568,14 +574,14 @@ public open class DoubleTensorAlgebra : */ public fun Tensor.rowsByIndices(indices: IntArray): DoubleTensor = stack(indices.map { this[it] }) - internal inline fun StructureND.fold(foldFunction: (DoubleArray) -> Double): Double = + private inline fun StructureND.fold(foldFunction: (DoubleArray) -> Double): Double = foldFunction(tensor.copyArray()) - internal inline fun StructureND.foldDim( - foldFunction: (DoubleArray) -> Double, + private inline fun StructureND.foldDim( dim: Int, keepDim: Boolean, - ): DoubleTensor { + foldFunction: (DoubleArray) -> R, + ): BufferedTensor { check(dim < dimension) { "Dimension $dim out of range $dimension" } val resShape = if (keepDim) { shape.take(dim).toIntArray() + intArrayOf(1) + shape.takeLast(dimension - dim - 1).toIntArray() @@ -583,80 +589,75 @@ public open class DoubleTensorAlgebra : shape.take(dim).toIntArray() + shape.takeLast(dimension - dim - 1).toIntArray() } val resNumElements = resShape.reduce(Int::times) - val resTensor = DoubleTensor(resShape, DoubleArray(resNumElements) { 0.0 }, 0) - for (index in resTensor.indices.asSequence()) { + val init = foldFunction(DoubleArray(1) { 0.0 }) + val resTensor = BufferedTensor(resShape, + MutableBuffer.auto(resNumElements) { init }, 0) + for (index in resTensor.indices) { val prefix = index.take(dim).toIntArray() val suffix = index.takeLast(dimension - dim - 1).toIntArray() resTensor[index] = foldFunction(DoubleArray(shape[dim]) { i -> tensor[prefix + intArrayOf(i) + suffix] }) } - return resTensor } override fun StructureND.sum(): Double = tensor.fold { it.sum() } override fun StructureND.sum(dim: Int, keepDim: Boolean): DoubleTensor = - foldDim({ x -> x.sum() }, dim, keepDim) + foldDim(dim, keepDim) { x -> x.sum() }.toDoubleTensor() override fun StructureND.min(): Double = this.fold { it.minOrNull()!! } override fun StructureND.min(dim: Int, keepDim: Boolean): DoubleTensor = - foldDim({ x -> x.minOrNull()!! }, dim, keepDim) + foldDim(dim, keepDim) { x -> x.minOrNull()!! }.toDoubleTensor() override fun StructureND.max(): Double = this.fold { it.maxOrNull()!! } override fun StructureND.max(dim: Int, keepDim: Boolean): DoubleTensor = - foldDim({ x -> x.maxOrNull()!! }, dim, keepDim) + foldDim(dim, keepDim) { x -> x.maxOrNull()!! }.toDoubleTensor() - override fun StructureND.argMax(dim: Int, keepDim: Boolean): DoubleTensor = - foldDim({ x -> - x.withIndex().maxByOrNull { it.value }?.index!!.toDouble() - }, dim, keepDim) + + override fun StructureND.argMax(dim: Int, keepDim: Boolean): IntTensor = + foldDim(dim, keepDim) { x -> + x.withIndex().maxByOrNull { it.value }?.index!! + }.toIntTensor() override fun StructureND.mean(): Double = this.fold { it.sum() / tensor.numElements } - override fun StructureND.mean(dim: Int, keepDim: Boolean): DoubleTensor = - foldDim( - { arr -> - check(dim < dimension) { "Dimension $dim out of range $dimension" } - arr.sum() / shape[dim] - }, - dim, - keepDim - ) + override fun StructureND.mean(dim: Int, keepDim: Boolean): DoubleTensor = foldDim(dim, keepDim) { arr -> + check(dim < dimension) { "Dimension $dim out of range $dimension" } + arr.sum() / shape[dim] + }.toDoubleTensor() - override fun StructureND.std(): Double = this.fold { arr -> + override fun StructureND.std(): Double = fold { arr -> val mean = arr.sum() / tensor.numElements sqrt(arr.sumOf { (it - mean) * (it - mean) } / (tensor.numElements - 1)) } override fun StructureND.std(dim: Int, keepDim: Boolean): DoubleTensor = foldDim( - { arr -> - check(dim < dimension) { "Dimension $dim out of range $dimension" } - val mean = arr.sum() / shape[dim] - sqrt(arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1)) - }, dim, keepDim - ) + ) { arr -> + check(dim < dimension) { "Dimension $dim out of range $dimension" } + val mean = arr.sum() / shape[dim] + sqrt(arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1)) + }.toDoubleTensor() - override fun StructureND.variance(): Double = this.fold { arr -> + override fun StructureND.variance(): Double = fold { arr -> val mean = arr.sum() / tensor.numElements arr.sumOf { (it - mean) * (it - mean) } / (tensor.numElements - 1) } override fun StructureND.variance(dim: Int, keepDim: Boolean): DoubleTensor = foldDim( - { arr -> - check(dim < dimension) { "Dimension $dim out of range $dimension" } - val mean = arr.sum() / shape[dim] - arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1) - }, dim, keepDim - ) + ) { arr -> + check(dim < dimension) { "Dimension $dim out of range $dimension" } + val mean = arr.sum() / shape[dim] + arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1) + }.toDoubleTensor() private fun cov(x: DoubleTensor, y: DoubleTensor): Double { val n = x.shape[0] diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/IntTensor.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/IntTensor.kt index dd9f9c0c1..715d9035f 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/IntTensor.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/IntTensor.kt @@ -6,6 +6,7 @@ package space.kscience.kmath.tensors.core import space.kscience.kmath.structures.IntBuffer +import space.kscience.kmath.tensors.core.internal.array /** * Default [BufferedTensor] implementation for [Int] values @@ -14,4 +15,7 @@ public class IntTensor internal constructor( shape: IntArray, buffer: IntArray, offset: Int = 0 -) : BufferedTensor(shape, IntBuffer(buffer), offset) +) : BufferedTensor(shape, IntBuffer(buffer), offset){ + public fun asDouble() : DoubleTensor = + DoubleTensor(shape, mutableBuffer.array().map{ it.toDouble()}.toDoubleArray(), bufferStart) +} diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/TensorLinearStructure.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/TensorLinearStructure.kt new file mode 100644 index 000000000..19eefc2f8 --- /dev/null +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/TensorLinearStructure.kt @@ -0,0 +1,74 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.kmath.tensors.core + +import space.kscience.kmath.nd.Strides +import kotlin.math.max + +/** + * This [Strides] implementation follows the last dimension first convention + * For more information: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.strides.html + * + * @param shape the shape of the tensor. + */ +public class TensorLinearStructure(override val shape: IntArray) : Strides() { + override val strides: IntArray + get() = stridesFromShape(shape) + + override fun index(offset: Int): IntArray = + indexFromOffset(offset, strides, shape.size) + + override val linearSize: Int + get() = shape.reduce(Int::times) + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other == null || this::class != other::class) return false + + other as TensorLinearStructure + + if (!shape.contentEquals(other.shape)) return false + + return true + } + + override fun hashCode(): Int { + return shape.contentHashCode() + } + + public companion object { + + public fun stridesFromShape(shape: IntArray): IntArray { + val nDim = shape.size + val res = IntArray(nDim) + if (nDim == 0) + return res + + var current = nDim - 1 + res[current] = 1 + + while (current > 0) { + res[current - 1] = max(1, shape[current]) * res[current] + current-- + } + return res + } + + public fun indexFromOffset(offset: Int, strides: IntArray, nDim: Int): IntArray { + val res = IntArray(nDim) + var current = offset + var strideIndex = 0 + + while (strideIndex < nDim) { + res[strideIndex] = (current / strides[strideIndex]) + current %= strides[strideIndex] + strideIndex++ + } + return res + } + } + +} \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/TensorLinearStructure.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/TensorLinearStructure.kt deleted file mode 100644 index 57668722a..000000000 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/TensorLinearStructure.kt +++ /dev/null @@ -1,71 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. - */ - -package space.kscience.kmath.tensors.core.internal - -import space.kscience.kmath.nd.Strides -import kotlin.math.max - - -internal fun stridesFromShape(shape: IntArray): IntArray { - val nDim = shape.size - val res = IntArray(nDim) - if (nDim == 0) - return res - - var current = nDim - 1 - res[current] = 1 - - while (current > 0) { - res[current - 1] = max(1, shape[current]) * res[current] - current-- - } - return res -} - -internal fun indexFromOffset(offset: Int, strides: IntArray, nDim: Int): IntArray { - val res = IntArray(nDim) - var current = offset - var strideIndex = 0 - - while (strideIndex < nDim) { - res[strideIndex] = (current / strides[strideIndex]) - current %= strides[strideIndex] - strideIndex++ - } - return res -} - -/** - * This [Strides] implementation follows the last dimension first convention - * For more information: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.strides.html - * - * @param shape the shape of the tensor. - */ -internal class TensorLinearStructure(override val shape: IntArray) : Strides() { - override val strides: IntArray - get() = stridesFromShape(shape) - - override fun index(offset: Int): IntArray = - indexFromOffset(offset, strides, shape.size) - - override val linearSize: Int - get() = shape.reduce(Int::times) - - override fun equals(other: Any?): Boolean { - if (this === other) return true - if (other == null || this::class != other::class) return false - - other as TensorLinearStructure - - if (!shape.contentEquals(other.shape)) return false - - return true - } - - override fun hashCode(): Int { - return shape.contentHashCode() - } -} diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt index d31e02677..290809cfd 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt @@ -53,7 +53,7 @@ internal val BufferedTensor.matrices: VirtualBuffer> internal fun BufferedTensor.matrixSequence(): Sequence> = matrices.asSequence() -internal fun dotHelper( +internal fun dotTo( a: MutableStructure2D, b: MutableStructure2D, res: MutableStructure2D, diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/tensorCastsUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/tensorCastsUtils.kt index d3d327f66..602430b03 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/tensorCastsUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/tensorCastsUtils.kt @@ -12,6 +12,7 @@ import space.kscience.kmath.tensors.api.Tensor import space.kscience.kmath.tensors.core.BufferedTensor import space.kscience.kmath.tensors.core.DoubleTensor import space.kscience.kmath.tensors.core.IntTensor +import space.kscience.kmath.tensors.core.TensorLinearStructure internal fun BufferedTensor.asTensor(): IntTensor = IntTensor(this.shape, this.mutableBuffer.array(), this.bufferStart) diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorAlgebraExtensions.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorAlgebraExtensions.kt index 916388ba9..1e6dfd52e 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorAlgebraExtensions.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorAlgebraExtensions.kt @@ -3,8 +3,11 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. */ +@file:OptIn(PerformancePitfall::class) + package space.kscience.kmath.tensors.core +import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.Shape import kotlin.jvm.JvmName diff --git a/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorFieldOpsND.kt b/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorFieldOpsND.kt index c50404c9c..e43bbbc6f 100644 --- a/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorFieldOpsND.kt +++ b/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorFieldOpsND.kt @@ -6,17 +6,20 @@ package space.kscience.kmath.viktor import org.jetbrains.bio.viktor.F64Array +import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.* import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.ExtendedFieldOps import space.kscience.kmath.operations.NumbersAddOps +import space.kscience.kmath.operations.PowerOperations @OptIn(UnstableKMathAPI::class) @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") public open class ViktorFieldOpsND : FieldOpsND, - ExtendedFieldOps> { + ExtendedFieldOps>, + PowerOperations> { public val StructureND.f64Buffer: F64Array get() = when (this) { @@ -35,6 +38,7 @@ public open class ViktorFieldOpsND : override fun StructureND.unaryMinus(): StructureND = -1 * this + @PerformancePitfall override fun StructureND.map(transform: DoubleField.(Double) -> Double): ViktorStructureND = F64Array(*shape).apply { DefaultStrides(shape).asSequence().forEach { index -> @@ -42,6 +46,7 @@ public open class ViktorFieldOpsND : } }.asStructure() + @PerformancePitfall override fun StructureND.mapIndexed( transform: DoubleField.(index: IntArray, Double) -> Double, ): ViktorStructureND = F64Array(*shape).apply { @@ -50,6 +55,7 @@ public open class ViktorFieldOpsND : } }.asStructure() + @PerformancePitfall override fun zip( left: StructureND, right: StructureND, @@ -110,7 +116,7 @@ public open class ViktorFieldOpsND : public val DoubleField.viktorAlgebra: ViktorFieldOpsND get() = ViktorFieldOpsND public open class ViktorFieldND( - override val shape: Shape + override val shape: Shape, ) : ViktorFieldOpsND(), FieldND, NumbersAddOps> { override val zero: ViktorStructureND by lazy { F64Array.full(init = 0.0, shape = shape).asStructure() } override val one: ViktorStructureND by lazy { F64Array.full(init = 1.0, shape = shape).asStructure() }