Merge pull request #63 from shahriariravanian/main

upgrade to ver 1.2 with the addition of symbolic and definite integration.

Author: shahriariravanian

Project.toml

@@ -1,7 +1,7 @@
 name = "SymbolicNumericIntegration"
 uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 authors = ["Shahriar Iravanian <[email protected]>"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"

README.md

@@ -21,16 +21,19 @@ the documentation, which contains the unreleased features.
 ## Example
 
 ```julia
-using Symbolics
-using SymbolicNumericIntegration
+julia> using SymbolicNumericIntegration
+julia> using Symbolics
 
-@variables x
+julia> @variables x a b
 
-integrate(3x^3 + 2x - 5)
-```
+julia> integrate(3x^3 + 2x - 5)
+(x^2 + (3//4)*(x^4) - (5//1)*x, 0, 0)
 
-```
-(x^2 + (3//4)*(x^4) - (5x), 0, 0)
+julia> integrate(exp(a * x), x; symbolic = true)
+(exp(a*x) / a, 0, 0)
+
+julia> integrate(sin(a * x) * cos(b * x), x; symbolic = true, detailed = false)
+(-a*cos(a*x)*cos(b*x) - b*sin(a*x)*sin(b*x)) / (a^2 - (b^2))
 ```
 
 # Citation

docs/make.jl

@@ -6,21 +6,14 @@ using Documenter, SymbolicNumericIntegration
 include("pages.jl")
 
 makedocs(sitename = "SymbolicNumericIntegration.jl",
-         authors = "Shahriar Iravanian",
-         modules = [SymbolicNumericIntegration],
-         clean = true, doctest = false, linkcheck = true,
-         strict = [
-             :doctest,
-             :linkcheck,
-             :parse_error,
-             :example_block,
-             # Other available options are
-             # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
-         ],
-         format = Documenter.HTML(analytics = "UA-90474609-3",
-                                  assets = ["assets/favicon.ico"],
-                                  canonical = "https://docs.sciml.ai/SymbolicNumericIntegration/stable/"),
-         pages = pages)
+    authors = "Shahriar Iravanian",
+    modules = [SymbolicNumericIntegration],
+    clean = true, doctest = false, linkcheck = true,
+    warnonly = true,
+    format = Documenter.HTML(analytics = "UA-90474609-3",
+        assets = ["assets/favicon.ico"],
+        canonical = "https://docs.sciml.ai/SymbolicNumericIntegration/stable/"),
+    pages = pages)
 
 deploydocs(repo = "github.com/SciML/SymbolicNumericIntegration.jl.git";
-           push_preview = true)
+    push_preview = true)

docs/src/index.md

@@ -2,10 +2,13 @@
 
 **SymbolicNumericIntegration.jl** is a hybrid symbolic/numerical integration package that works on the [Julia Symbolics](https://docs.sciml.ai/Symbolics/stable/) expressions.
 
-**SymbolicNumericIntegration.jl** uses a randomized algorithm based on a hybrid of the *method of undetermined coefficients* and *sparse regression* and can solve a large subset of basic standard integrals (polynomials, exponential/logarithmic, trigonometric and hyperbolic, inverse trigonometric and hyperbolic, rational and square root).
-The basis of how it works and the theory of integration using the Symbolic-Numeric methods refer to [Basis of Symbolic-Numeric Integration](https://github.com/SciML/SymbolicNumericIntegration.jl/blob/main/docs/theory.ipynb).
+**SymbolicNumericIntegration.jl** uses a randomized algorithm based on a hybrid of the *method of undetermined coefficients* and *sparse regression* and can solve a large subset of basic standard integrals (polynomials, exponential/logarithmic, trigonometric and hyperbolic, inverse trigonometric and hyperbolic, rational and square root). The symbolic part of the algorithm is similar (but not identical) to Risch-Bronstein's poor man's integrator and generates a list of ansatzes (candidate terms). The numerical part uses sparse regression adopted from the Sparse identification of nonlinear dynamics (SINDy) algorithm to prune down the ansatzes and find the corresponding coefficients. The basis of how it works and the theory of integration using the Symbolic-Numeric methods refer to [Basis of Symbolic-Numeric Integration](https://github.com/SciML/SymbolicNumericIntegration.jl/blob/main/docs/theory.ipynb).
 
-Function `integrate` returns the integral of a univariate expression with *constant* real or complex coefficients. `integrate` returns a tuple with three values. The first one is the solved integral, the second one is the sum of the unsolved terms, and the third value is the residual error. If `integrate` is successful, the unsolved portion is reported as 0.
+[hyint](https://github.com/siravan/hyint) is the python counterpart of **SymbolicNumericIntegration.jl** that works with **sympy** expressions.
+
+Originally, **SymbolicNumericIntegration.jl** expected only univariate expression with *constant* real or complex coefficients as input. As of version 1.2, `integrate` function accepts symbolic constants for a subset of solvable integrals.
+
+`integrate` returns a tuple with three values. The first one is the solved integral, the second one is the sum of the unsolved terms, and the third value is the residual error. If `integrate` is successful, the unsolved portion is reported as 0. If we pass `detailed=false` to `integrate', the output is simplified to only the resulting integrals. In this case, `nothing` is returned if the integral cannot be solved. Note that the simplified output will become the default in a future version.
 
 ## Installation
 
@@ -19,57 +22,53 @@ Pkg.add("SymbolicNumericIntegration")
 Examples:
 
 ```julia
-using Symbolics
-using SymbolicNumericIntegration
+julia> using SymbolicNumericIntegration
+julia> using Symbolics
 
-@variables x
-```
+julia> @variables x a b
+
+# if `detailed = true` (default), the output is a tuple of (solution, unsolved portion, err)
 
-```julia
 julia> integrate(3x^3 + 2x - 5)
 (x^2 + (3//4)*(x^4) - (5x), 0, 0)
 
 julia> integrate((5 + 2x)^-1)
 ((1//2)*log((5//2) + x), 0, 0.0)
 
-julia> integrate(1 / (6 + x^2 - (5x)))
-(log(x - 3) - log(x - 2), 0, 3.339372764128952e-16)
+# `detailed = false` simplifies the output to just the resulting integral
 
-julia> integrate(1 / (x^2 - 16))
-((1//8)*log(x - 4) - ((1//8)*log(4 + x)), 0, 1.546926788028958e-16)
+julia> integrate(x^2 / (16 + x^2); detailed = false)
+x + 4atan((-1//4)*x)
 
-julia> integrate(x^2 / (16 + x^2))
-(x + 4atan((-1//4)*x), 0, 1.3318788420751984e-16)
+julia> integrate(x^2 * log(x); detailed = false)
+(1//3)*(x^3)*log(x) - (1//9)*(x^3)
 
-julia> integrate(x^2 / sqrt(4 + x^2))
-((1//2)*x*((4 + x^2)^0.5) - ((2//1)*log(x + sqrt(4 + x^2))), 0, 8.702422633074313e-17)
+julia> integrate(sec(x) * tan(x); detailed = false)
+sec(x)
 
-julia> integrate(x^2 * log(x))
-((1//3)*log(x)*(x^3) - ((1//9)*(x^3)), 0, 0)
+# Symbolic integration. Here, a is a symbolic constant; therefore, we need 
+# to explicitly define the independent variable (say, x). Also, we set
+# `symbolic = true` to force using the symbolic solver
 
-julia> integrate(x^2 * exp(x))
-(2exp(x) + exp(x)*(x^2) - (2x*exp(x)), 0, 0)
+julia> integrate(sin(a * x), x; detailed = false, symbolic = true)
+(-cos(a*x)) / a
 
-julia> integrate(tan(2x))
-((-1//2)*log(cos(2x)), 0, 0)
+julia> integrate(x^2 * cos(a * x), x; detailed = false, symbolic = true)
+((a^2)*(x^2)*sin(a*x) + 2.0a*x*cos(a*x) - 2.0sin(a*x)) / (a^3)
 
-julia> integrate(sec(x) * tan(x))
-(cos(x)^-1, 0, 0)
+julia> integrate(log(log(a * x)) / x, x; detailed = false, symbolic = true)
+log(a*x)*log(log(a*x)) - log(a*x)
 
-julia> integrate(cosh(2x) * exp(x))
-((2//3)*exp(x)*sinh(2x) - ((1//3)*exp(x)*cosh(2x)), 0, 7.073930088880992e-8)
+# multiple symbolic constants
 
-julia> integrate(cosh(x) * sin(x))
-((1//2)*sin(x)*sinh(x) - ((1//2)*cos(x)*cosh(x)), 0, 4.8956233716268386e-17)
+julia> integrate(cosh(a * x) * exp(b * x), x; detailed = false, symbolic = true)
+(a*sinh(a*x)*exp(b*x) - b*cosh(a*x)*exp(b*x)) / (a^2 - (b^2))
 
-julia> integrate(cosh(2x) * sin(3x))
-(0.153845sinh(2x)*sin(3x) - (0.23077cosh(2x)*cos(3x)), 0, 4.9807620877373405e-6)
+# definite integration, passing a tuple of (x, lower bound, higher bound) in the 
+# second argument
 
-julia> integrate(log(log(x)) * (x^-1))
-(log(x)*log(log(x)) - log(x), 0, 0)
-
-julia> integrate(exp(x^2))
-(0, exp(x^2), Inf)    # as expected!
+julia> integrate(x * sin(a * x), (x, 0, 1); symbolic = true, detailed = false)
+(sin(a) - a*cos(a)) / (a^2)
 ```
 
 SymbolicNumericIntegration.jl exports some special integral functions (defined over Complex numbers) and uses them in solving integrals:
@@ -81,15 +80,15 @@ SymbolicNumericIntegration.jl exports some special integral functions (defined o
 
 For examples:
 
-```
+```julia
 julia> integrate(exp(x + 1) / (x + 1))
-(SymbolicNumericIntegration.Ei(1 + x), 0, 1.1796119636642288e-16)
+(SymbolicNumericIntegration.Ei(1 + x), 0, 0.0)
 
-julia> integrate(x * cos(x^2 - 1) / (x^2 - 1))
-((1//2)*SymbolicNumericIntegration.Ci(x^2 - 1), 0, 2.7755575615628914e-17)
+julia> integrate(x * cos(a*x^2 - 1) / (a*x^2 - 1), x; detailed=false, symbolic=true)
+((1//2)*SymbolicNumericIntegration.Ci(a*(x^2) - 1)) / a
 
-julia> integrate(1 / (x*log(log(x))))
-(SymbolicNumericIntegration.Li(log(x)), 0, 1.1102230246251565e-16)
+julia> integrate(1 / (x*log(log(x))), x; detailed=false, symbolic=true)
+SymbolicNumericIntegration.Li(log(x))
 ```
 
 ```@docs
@@ -102,7 +101,7 @@ integrate(eq, x; kwargs...)
 
 Additionally, 12 test suites from the *Rule-based Integrator* ([Rubi](https://rulebasedintegration.org/)) are included in the `/test` directory. For example, we can test the first one as below. *Axiom* refers to the format of the test files)
 
-```julia
+```
 using SymbolicNumericIntegration
 include("test/axiom.jl")  # note, you may need to use the correct path
 
@@ -198,33 +197,3 @@ Pkg.status(; mode = PKGMODE_MANIFEST) # hide
 ```@raw html
 </details>
 ```
-
-```@raw html
-You can also download the 
-<a href="
-```
-
-```@eval
-using TOML
-version = TOML.parse(read("../../Project.toml", String))["version"]
-name = TOML.parse(read("../../Project.toml", String))["name"]
-link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
-       "/assets/Manifest.toml"
-```
-
-```@raw html
-">manifest</a> file and the
-<a href="
-```
-
-```@eval
-using TOML
-version = TOML.parse(read("../../Project.toml", String))["version"]
-name = TOML.parse(read("../../Project.toml", String))["name"]
-link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
-       "/assets/Project.toml"
-```
-
-```@raw html
-">project</a> file.
-```

src/SymbolicNumericIntegration.jl

@@ -10,6 +10,7 @@ using SymbolicUtils: exprtype, BasicSymbolic
 using DataDrivenDiffEq, DataDrivenSparse
 
 include("utils.jl")
+include("tree.jl")
 include("special.jl")
 
 include("cache.jl")
@@ -29,6 +30,5 @@ include("integral.jl")
 export integrate, generate_basis
 
 include("symbolic.jl")
-include("logger.jl")
 
 end # module

src/candidates.jl

@@ -64,72 +64,6 @@ function closure(eq, x; max_terms = 50)
     unique([[s for s in S]; [s * x for s in S]])
 end
 
-function candidate_pow_minus(p, k; abstol = 1e-8)
-    if isnan(poly_deg(p))
-        return [p^k, p^(k + 1), log(p)]
-    end
-
-    x = var(p)
-    d = poly_deg(p)
-
-    for j in 1:10  # will try 10 times to find the roots
-        r, s = find_roots(p, x)
-        if length(r) + length(s) >= d
-            break
-        end
-    end
-    s = s[1:2:end]
-    r = nice_parameter.(r)
-    s = nice_parameter.(s)
-
-    # ∫ 1 / ((x-z₁)(x-z₂)) dx = ... + c₁ * log(x-z₁) + c₂ * log(x-z₂)
-    q = Any[log(x - u) for u in r]
-    for i in eachindex(s)
-        β = s[i]
-        if abs(imag(β)) > abstol
-            push!(q, atan((x - real(β)) / imag(β)))
-            push!(q, (1 + x) * log(x^2 - 2 * real(β) * x + abs2(β)))
-        else
-            push!(q, log(x - real(β)))
-        end
-    end
-    q = unique(q)
-
-    # return [[p^k, p^(k+1)]; candidates(q₁, x)]
-    if k ≈ -1
-        return [[p^k]; q]
-    else
-        return [[p^k, p^(k + 1)]; q]
-    end
-end
-
-function candidate_sqrt(p, k)
-    x = var(p)
-
-    h = Any[p^k, p^(k + 1)]
-
-    if poly_deg(p) == 2
-        r, s = find_roots(p, x)
-        l = leading(p, x)
-        if length(r) == 2
-            if sum(r) ≈ 0
-                r₁ = abs(r[1])
-                if l > 0
-                    push!(h, acosh(x / r₁))
-                else
-                    push!(h, asin(x / r₁))
-                end
-            end
-        elseif real(s[1]) ≈ 0
-            push!(h, asinh(x / imag.(s[1])))
-        end
-    end
-
-    Δ = expand_derivatives(Differential(x)(p))
-    push!(h, log(0.5 * Δ + sqrt(p)))
-    return h
-end
-
 ###############################################################################
 
 enqueue_expr!(S, q, eq, x) = enqueue_expr!!(S, q, ops(eq)..., x)