Improved trigonometry operations #99

Closed
CommanderTvis wants to merge 22 commits from adv-expr-improved-trigonometry into adv-expr
19 changed files with 333 additions and 99 deletions
Showing only changes of commit 76ad027460 - Show all commits

View File

@ -0,0 +1,27 @@
plugins {
id("scientifik.mpp")
}
repositories{
maven("https://dl.bintray.com/hotkeytlt/maven")
}
kotlin.sourceSets {
commonMain {
dependencies {
api(project(":kmath-core"))
implementation("com.github.h0tk3y.betterParse:better-parse-multiplatform:0.4.0-alpha-3")
implementation("com.github.h0tk3y.betterParse:better-parse-multiplatform-metadata:0.4.0-alpha-3")
}
}
jvmMain{
dependencies{
implementation("com.github.h0tk3y.betterParse:better-parse-jvm:0.4.0-alpha-3")
}
}
jsMain{
dependencies{
implementation("com.github.h0tk3y.betterParse:better-parse-js:0.4.0-alpha-3")
}
}
}

View File

@ -0,0 +1,62 @@
package scientifik.kmath.ast
import scientifik.kmath.operations.NumericAlgebra
import scientifik.kmath.operations.RealField
/**
* A Mathematical Syntax Tree node for mathematical expressions
*/
sealed class MST {
/**
* A node containing unparsed string
*/
data class Singular(val value: String) : MST()
/**
* A node containing a number
*/
data class Numeric(val value: Number) : MST()
/**
* A node containing an unary operation
*/
data class Unary(val operation: String, val value: MST) : MST() {
companion object {
const val ABS_OPERATION = "abs"
//TODO add operations
}
}
/**
* A node containing binary operation
*/
data class Binary(val operation: String, val left: MST, val right: MST) : MST() {
companion object
}
}
//TODO add a function with positional arguments
//TODO add a function with named arguments
fun <T> NumericAlgebra<T>.evaluate(node: MST): T {
return when (node) {
is MST.Numeric -> number(node.value)
is MST.Singular -> symbol(node.value)
is MST.Unary -> unaryOperation(node.operation, evaluate(node.value))
is MST.Binary -> when {
node.left is MST.Numeric && node.right is MST.Numeric -> {
val number = RealField.binaryOperation(
node.operation,
node.left.value.toDouble(),
node.right.value.toDouble()
)
number(number)
}
node.left is MST.Numeric -> leftSideNumberOperation(node.operation, node.left.value, evaluate(node.right))
node.right is MST.Numeric -> rightSideNumberOperation(node.operation, evaluate(node.left), node.right.value)
else -> binaryOperation(node.operation, evaluate(node.left), evaluate(node.right))
}
}
}

View File

@ -0,0 +1,19 @@
package scientifik.kmath.ast
import scientifik.kmath.expressions.Expression
import scientifik.kmath.operations.NumericAlgebra
/**
* The expression evaluates MST on-flight. Should be much faster than functional expression, but slower than ASM-generated expressions.
*/
class MSTExpression<T>(val algebra: NumericAlgebra<T>, val mst: MST) : Expression<T> {
/**
* Substitute algebra raw value
*/
private inner class InnerAlgebra(val arguments: Map<String, T>) : NumericAlgebra<T> by algebra {
override fun symbol(value: String): T = arguments[value] ?: super.symbol(value)
}
override fun invoke(arguments: Map<String, T>): T = InnerAlgebra(arguments).evaluate(mst)
}

View File

@ -0,0 +1,59 @@
package scientifik.kmath.ast
import com.github.h0tk3y.betterParse.combinators.*
import com.github.h0tk3y.betterParse.grammar.Grammar
import com.github.h0tk3y.betterParse.grammar.parseToEnd
import com.github.h0tk3y.betterParse.grammar.parser
import com.github.h0tk3y.betterParse.grammar.tryParseToEnd
import com.github.h0tk3y.betterParse.parser.ParseResult
import com.github.h0tk3y.betterParse.parser.Parser
import scientifik.kmath.operations.FieldOperations
import scientifik.kmath.operations.PowerOperations
import scientifik.kmath.operations.RingOperations
import scientifik.kmath.operations.SpaceOperations
/**
* TODO move to common
*/
private object ArithmeticsEvaluator : Grammar<MST>() {
val num by token("-?[\\d.]+(?:[eE]-?\\d+)?")
val lpar by token("\\(")
val rpar by token("\\)")
val mul by token("\\*")
val pow by token("\\^")
val div by token("/")
val minus by token("-")
val plus by token("\\+")
val ws by token("\\s+", ignore = true)
val number: Parser<MST> by num use { MST.Numeric(text.toDouble()) }
val term: Parser<MST> by number or
(skip(minus) and parser(this::term) map { MST.Unary(SpaceOperations.MINUS_OPERATION, it) }) or
(skip(lpar) and parser(this::rootParser) and skip(rpar))
val powChain by leftAssociative(term, pow) { a, _, b ->
MST.Binary(PowerOperations.POW_OPERATION, a, b)
}
val divMulChain: Parser<MST> by leftAssociative(powChain, div or mul use { type }) { a, op, b ->
if (op == div) {
MST.Binary(FieldOperations.DIV_OPERATION, a, b)
} else {
MST.Binary(RingOperations.TIMES_OPERATION, a, b)
}
}
val subSumChain: Parser<MST> by leftAssociative(divMulChain, plus or minus use { type }) { a, op, b ->
if (op == plus) {
MST.Binary(SpaceOperations.PLUS_OPERATION, a, b)
} else {
MST.Binary(SpaceOperations.MINUS_OPERATION, a, b)
}
}
override val rootParser: Parser<MST> by subSumChain
}
fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryParseToEnd(this)
fun String.parseMath(): MST = ArithmeticsEvaluator.parseToEnd(this)

View File

@ -0,0 +1,23 @@
package scientifik.kmath.ast
import scientifik.kmath.expressions.Expression
import scientifik.kmath.operations.Algebra
import scientifik.kmath.operations.NumericAlgebra
//TODO stubs for asm generation
interface AsmExpression<T>
interface AsmExpressionAlgebra<T, A : Algebra<T>> : NumericAlgebra<AsmExpression<T>> {
val algebra: A
}
fun <T> AsmExpression<T>.compile(): Expression<T> = TODO()
//TODO add converter for functional expressions
inline fun <reified T : Any, A : Algebra<T>> A.asm(
block: AsmExpressionAlgebra<T, A>.() -> AsmExpression<T>
): Expression<T> = TODO()
inline fun <reified T : Any, A : Algebra<T>> A.asm(ast: MST): Expression<T> = asm { evaluate(ast) }

View File

@ -0,0 +1,17 @@
package scietifik.kmath.ast
import org.junit.jupiter.api.Assertions.assertEquals
import org.junit.jupiter.api.Test
import scientifik.kmath.ast.evaluate
import scientifik.kmath.ast.parseMath
import scientifik.kmath.operations.Complex
import scientifik.kmath.operations.ComplexField
internal class ParserTest{
@Test
fun parsedExpression(){
val mst = "2+2*(2+2)".parseMath()
val res = ComplexField.evaluate(mst)
assertEquals(Complex(10.0,0.0), res)
}
}

View File

@ -2,7 +2,7 @@ package scientifik.kmath.commons.expressions
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
import scientifik.kmath.expressions.Expression import scientifik.kmath.expressions.Expression
import scientifik.kmath.expressions.ExpressionContext import scientifik.kmath.expressions.ExpressionAlgebra
import scientifik.kmath.operations.ExtendedField import scientifik.kmath.operations.ExtendedField
import scientifik.kmath.operations.Field import scientifik.kmath.operations.Field
import kotlin.properties.ReadOnlyProperty import kotlin.properties.ReadOnlyProperty
@ -115,7 +115,7 @@ fun DiffExpression.derivative(name: String) = derivative(name to 1)
/** /**
* A context for [DiffExpression] (not to be confused with [DerivativeStructure]) * A context for [DiffExpression] (not to be confused with [DerivativeStructure])
*/ */
object DiffExpressionContext : ExpressionContext<Double, DiffExpression>, Field<DiffExpression> { object DiffExpressionAlgebra : ExpressionAlgebra<Double, DiffExpression>, Field<DiffExpression> {
override fun variable(name: String, default: Double?) = override fun variable(name: String, default: Double?) =
DiffExpression { variable(name, default?.const()) } DiffExpression { variable(name, default?.const()) }

View File

@ -7,6 +7,12 @@ import scientifik.kmath.operations.Algebra
*/ */
interface Expression<T> { interface Expression<T> {
operator fun invoke(arguments: Map<String, T>): T operator fun invoke(arguments: Map<String, T>): T
companion object {
operator fun <T> invoke(block: (Map<String, T>) -> T): Expression<T> = object : Expression<T> {
override fun invoke(arguments: Map<String, T>): T = block(arguments)
}
}
} }
operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs)) operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs))
@ -14,7 +20,7 @@ operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke
/** /**
* A context for expression construction * A context for expression construction
*/ */
interface ExpressionContext<T, E> : Algebra<E> { interface ExpressionAlgebra<T, E> : Algebra<E> {
/** /**
* Introduce a variable into expression context * Introduce a variable into expression context
*/ */

View File

@ -1,52 +0,0 @@
package scientifik.kmath.expressions
import scientifik.kmath.operations.NumericAlgebra
/**
* A syntax tree node for mathematical expressions
*/
sealed class SyntaxTreeNode
/**
* A node containing unparsed string
*/
data class SingularNode(val value: String) : SyntaxTreeNode()
/**
* A node containing a number
*/
data class NumberNode(val value: Number) : SyntaxTreeNode()
/**
* A node containing an unary operation
*/
data class UnaryNode(val operation: String, val value: SyntaxTreeNode) : SyntaxTreeNode() {
companion object {
const val ABS_OPERATION = "abs"
const val SIN_OPERATION = "sin"
const val COS_OPERATION = "cos"
const val EXP_OPERATION = "exp"
const val LN_OPERATION = "ln"
//TODO add operations
}
}
/**
* A node containing binary operation
*/
data class BinaryNode(val operation: String, val left: SyntaxTreeNode, val right: SyntaxTreeNode) : SyntaxTreeNode() {
companion object
}
//TODO add a function with positional arguments
//TODO add a function with named arguments
fun <T> NumericAlgebra<T>.compile(node: SyntaxTreeNode): T{
return when (node) {
is NumberNode -> number(node.value)
is SingularNode -> raw(node.value)
is UnaryNode -> unaryOperation(node.operation, compile(node.value))
is BinaryNode -> binaryOperation(node.operation, compile(node.left), compile(node.right))
}
}

View File

@ -21,8 +21,11 @@ internal class SumExpression<T>(
override fun invoke(arguments: Map<String, T>): T = context.add(first.invoke(arguments), second.invoke(arguments)) override fun invoke(arguments: Map<String, T>): T = context.add(first.invoke(arguments), second.invoke(arguments))
} }
internal class ProductExpression<T>(val context: Ring<T>, val first: Expression<T>, val second: Expression<T>) : internal class ProductExpression<T>(
Expression<T> { val context: Ring<T>,
val first: Expression<T>,
val second: Expression<T>
) : Expression<T> {
override fun invoke(arguments: Map<String, T>): T = override fun invoke(arguments: Map<String, T>): T =
context.multiply(first.invoke(arguments), second.invoke(arguments)) context.multiply(first.invoke(arguments), second.invoke(arguments))
} }
@ -39,7 +42,7 @@ internal class DivExpession<T>(val context: Field<T>, val expr: Expression<T>, v
open class FunctionalExpressionSpace<T>( open class FunctionalExpressionSpace<T>(
val space: Space<T> val space: Space<T>
) : Space<Expression<T>>, ExpressionContext<T,Expression<T>> { ) : Space<Expression<T>>, ExpressionAlgebra<T, Expression<T>> {
override val zero: Expression<T> = ConstantExpression(space.zero) override val zero: Expression<T> = ConstantExpression(space.zero)
@ -61,10 +64,9 @@ open class FunctionalExpressionSpace<T>(
open class FunctionalExpressionField<T>( open class FunctionalExpressionField<T>(
val field: Field<T> val field: Field<T>
) : Field<Expression<T>>, ExpressionContext<T,Expression<T>>, FunctionalExpressionSpace<T>(field) { ) : Field<Expression<T>>, ExpressionAlgebra<T, Expression<T>>, FunctionalExpressionSpace<T>(field) {
override val one: Expression<T> override val one: Expression<T> = ConstantExpression(this.field.one)
get() = const(this.field.one)
fun const(value: Double): Expression<T> = const(field.run { one * value }) fun const(value: Double): Expression<T> = const(field.run { one * value })
@ -78,3 +80,15 @@ open class FunctionalExpressionField<T>(
operator fun T.times(arg: Expression<T>) = arg * this operator fun T.times(arg: Expression<T>) = arg * this
operator fun T.div(arg: Expression<T>) = arg / this operator fun T.div(arg: Expression<T>) = arg / this
} }
/**
* Create a functional expression on this [Space]
*/
fun <T> Space<T>.buildExpression(block: FunctionalExpressionSpace<T>.() -> Expression<T>): Expression<T> =
FunctionalExpressionSpace(this).run(block)
/**
* Create a functional expression on this [Field]
*/
fun <T> Field<T>.buildExpression(block: FunctionalExpressionField<T>.() -> Expression<T>): Expression<T> =
FunctionalExpressionField(this).run(block)

View File

@ -128,14 +128,14 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
luRow[col] = sum luRow[col] = sum
// maintain best permutation choice // maintain best permutation choice
if (abs(sum) > largest) { if (this@lup.abs(sum) > largest) {
largest = abs(sum) largest = this@lup.abs(sum)
max = row max = row
} }
} }
// Singularity check // Singularity check
if (checkSingular(abs(lu[max, col]))) { if (checkSingular(this@lup.abs(lu[max, col]))) {
error("The matrix is singular") error("The matrix is singular")
} }

View File

@ -1,5 +1,7 @@
package scientifik.kmath.misc package scientifik.kmath.misc
import kotlin.math.abs
/** /**
* Convert double range to sequence. * Convert double range to sequence.
* *
@ -8,8 +10,7 @@ package scientifik.kmath.misc
* *
* If step is negative, the same goes from upper boundary downwards * If step is negative, the same goes from upper boundary downwards
*/ */
fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double> = fun ClosedFloatingPointRange<Double>.toSequenceWithStep(step: Double): Sequence<Double> = when {
when {
step == 0.0 -> error("Zero step in double progression") step == 0.0 -> error("Zero step in double progression")
step > 0 -> sequence { step > 0 -> sequence {
var current = start var current = start
@ -27,9 +28,18 @@ fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double>
} }
} }
/**
* Convert double range to sequence with the fixed number of points
*/
fun ClosedFloatingPointRange<Double>.toSequenceWithPoints(numPoints: Int): Sequence<Double> {
require(numPoints > 1) { "The number of points should be more than 2" }
return toSequenceWithStep(abs(endInclusive - start) / (numPoints - 1))
}
/** /**
* Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints] * Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints]
*/ */
@Deprecated("Replace by 'toSequenceWithPoints'")
fun ClosedFloatingPointRange<Double>.toGrid(numPoints: Int): DoubleArray { fun ClosedFloatingPointRange<Double>.toGrid(numPoints: Int): DoubleArray {
if (numPoints < 2) error("Can't create generic grid with less than two points") if (numPoints < 2) error("Can't create generic grid with less than two points")
return DoubleArray(numPoints) { i -> start + (endInclusive - start) / (numPoints - 1) * i } return DoubleArray(numPoints) { i -> start + (endInclusive - start) / (numPoints - 1) * i }

View File

@ -10,7 +10,10 @@ interface Algebra<T> {
/** /**
* Wrap raw string or variable * Wrap raw string or variable
*/ */
fun raw(value: String): T = error("Wrapping of '$value' is not supported in $this") fun symbol(value: String): T = error("Wrapping of '$value' is not supported in $this")
@Deprecated("Symbol is more concise",replaceWith = ReplaceWith("symbol"))
fun raw(value: String): T = symbol(value)
/** /**
* Dynamic call of unary operation with name [operation] on [arg] * Dynamic call of unary operation with name [operation] on [arg]
@ -31,6 +34,12 @@ interface NumericAlgebra<T> : Algebra<T> {
* Wrap a number * Wrap a number
*/ */
fun number(value: Number): T fun number(value: Number): T
fun leftSideNumberOperation(operation: String, left: Number, right: T): T =
binaryOperation(operation, number(left), right)
fun rightSideNumberOperation(operation: String, left: T, right: Number): T =
leftSideNumberOperation(operation, right, left)
} }
/** /**
@ -128,7 +137,13 @@ interface Ring<T> : Space<T>, RingOperations<T>, NumericAlgebra<T> {
override fun number(value: Number): T = one * value.toDouble() override fun number(value: Number): T = one * value.toDouble()
// those operators are blocked by type conflict in RealField override fun leftSideNumberOperation(operation: String, left: Number, right: T): T = when (operation) {
RingOperations.TIMES_OPERATION -> left * right
else -> super.leftSideNumberOperation(operation, left, right)
}
//TODO those operators are blocked by type conflict in RealField
// operator fun T.plus(b: Number) = this.plus(b * one) // operator fun T.plus(b: Number) = this.plus(b * one)
// operator fun Number.plus(b: T) = b + this // operator fun Number.plus(b: T) = b + this
// //

View File

@ -13,7 +13,7 @@ private val PI_DIV_2 = Complex(PI / 2, 0)
/** /**
* A field for complex numbers * A field for complex numbers
*/ */
object ComplexField : ExtendedFieldOperations<Complex>, Field<Complex> { object ComplexField : ExtendedField<Complex> {
override val zero: Complex = Complex(0.0, 0.0) override val zero: Complex = Complex(0.0, 0.0)
override val one: Complex = Complex(1.0, 0.0) override val one: Complex = Complex(1.0, 0.0)
@ -54,6 +54,12 @@ object ComplexField : ExtendedFieldOperations<Complex>, Field<Complex> {
operator fun Complex.minus(d: Double) = add(this, -d.toComplex()) operator fun Complex.minus(d: Double) = add(this, -d.toComplex())
operator fun Double.times(c: Complex) = Complex(c.re * this, c.im * this) operator fun Double.times(c: Complex) = Complex(c.re * this, c.im * this)
override fun symbol(value: String): Complex = if (value == "i") {
i
} else {
super.symbol(value)
}
} }
/** /**

View File

@ -10,10 +10,30 @@ interface ExtendedFieldOperations<T> :
FieldOperations<T>, FieldOperations<T>,
InverseTrigonometricOperations<T>, InverseTrigonometricOperations<T>,
PowerOperations<T>, PowerOperations<T>,
ExponentialOperations<T> ExponentialOperations<T> {
interface ExtendedField<T> : ExtendedFieldOperations<T>, Field<T> override fun tan(arg: T): T = sin(arg) / cos(arg)
override fun unaryOperation(operation: String, arg: T): T = when (operation) {
TrigonometricOperations.COS_OPERATION -> cos(arg)
TrigonometricOperations.SIN_OPERATION -> sin(arg)
PowerOperations.SQRT_OPERATION -> sqrt(arg)
ExponentialOperations.EXP_OPERATION -> exp(arg)
ExponentialOperations.LN_OPERATION -> ln(arg)
else -> super.unaryOperation(operation, arg)
}
}
interface ExtendedField<T> : ExtendedFieldOperations<T>, Field<T> {
override fun rightSideNumberOperation(operation: String, left: T, right: Number): T {
return when (operation) {
PowerOperations.POW_OPERATION -> power(left, right)
else -> super.rightSideNumberOperation(operation, left, right)
}
}
}
/** /**
* Real field element wrapping double. * Real field element wrapping double.
* *

View File

@ -49,9 +49,11 @@ fun <T : MathElement<out InverseTrigonometricOperations<T>>> atan(arg: T): T = a
interface PowerOperations<T> : Algebra<T> { interface PowerOperations<T> : Algebra<T> {
fun power(arg: T, pow: Number): T fun power(arg: T, pow: Number): T
fun sqrt(arg: T) = power(arg, 0.5) fun sqrt(arg: T) = power(arg, 0.5)
infix fun T.pow(pow: Number) = power(this, pow) infix fun T.pow(pow: Number) = power(this, pow)
companion object { companion object {
const val POW_OPERATION = "pow"
const val SQRT_OPERATION = "sqrt" const val SQRT_OPERATION = "sqrt"
} }
} }
@ -65,6 +67,11 @@ fun <T : MathElement<out PowerOperations<T>>> sqr(arg: T): T = arg pow 2.0
interface ExponentialOperations<T>: Algebra<T> { interface ExponentialOperations<T>: Algebra<T> {
fun exp(arg: T): T fun exp(arg: T): T
fun ln(arg: T): T fun ln(arg: T): T
companion object {
const val EXP_OPERATION = "exp"
const val LN_OPERATION = "ln"
}
} }
fun <T : MathElement<out ExponentialOperations<T>>> exp(arg: T): T = arg.context.exp(arg) fun <T : MathElement<out ExponentialOperations<T>>> exp(arg: T): T = arg.context.exp(arg)

View File

@ -1,13 +1,8 @@
package scientifik.kmath.structures package scientifik.kmath.structures
import scientifik.kmath.operations.* import scientifik.kmath.operations.ExtendedField
interface ExtendedNDField<T : Any, F, N : NDStructure<T>> : interface ExtendedNDField<T : Any, F : ExtendedField<T>, N : NDStructure<T>> : NDField<T, F, N>, ExtendedField<N>
NDField<T, F, N>,
InverseTrigonometricOperations<N>,
PowerOperations<N>,
ExponentialOperations<N>
where F : ExtendedFieldOperations<T>, F : Field<T>
///** ///**

View File

@ -1,10 +1,12 @@
package scientifik.kmath.operations package scientifik.kmath.operations
import scientifik.kmath.structures.*
import java.math.BigDecimal import java.math.BigDecimal
import java.math.BigInteger import java.math.BigInteger
import java.math.MathContext import java.math.MathContext
/**
* A field wrapper for Java [BigInteger]
*/
object JBigIntegerField : Field<BigInteger> { object JBigIntegerField : Field<BigInteger> {
override val zero: BigInteger = BigInteger.ZERO override val zero: BigInteger = BigInteger.ZERO
override val one: BigInteger = BigInteger.ONE override val one: BigInteger = BigInteger.ONE
@ -18,6 +20,9 @@ object JBigIntegerField : Field<BigInteger> {
override fun divide(a: BigInteger, b: BigInteger): BigInteger = a.div(b) override fun divide(a: BigInteger, b: BigInteger): BigInteger = a.div(b)
} }
/**
* A Field wrapper for Java [BigDecimal]
*/
class JBigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Field<BigDecimal> { class JBigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Field<BigDecimal> {
override val zero: BigDecimal = BigDecimal.ZERO override val zero: BigDecimal = BigDecimal.ZERO
override val one: BigDecimal = BigDecimal.ONE override val one: BigDecimal = BigDecimal.ONE

View File

@ -44,5 +44,6 @@ include(
":kmath-dimensions", ":kmath-dimensions",
":kmath-for-real", ":kmath-for-real",
":kmath-geometry", ":kmath-geometry",
":kmath-ast",
":examples" ":examples"
) )