ver 1.3; improved integration coverage; symbolic code refactored

Author: shahriariravanian

README.md

@@ -21,19 +21,18 @@ the documentation, which contains the unreleased features.
 ## Example
 
 ```julia
-julia> using SymbolicNumericIntegration
 julia> using Symbolics
+using SymbolicNumericIntegration
 
 julia> @variables x a b
 
 julia> integrate(3x^3 + 2x - 5)
-(x^2 + (3//4)*(x^4) - (5//1)*x, 0, 0)
 
 julia> integrate(exp(a * x), x; symbolic = true)
-(exp(a*x) / a, 0, 0)
+(x^2 + (3//4)*(x^4) - (5//1)*x, 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))
+(exp(a*x) / a, 0, 0)
 ```
 
 # Citation

docs/make.jl

@@ -9,7 +9,7 @@ makedocs(sitename = "SymbolicNumericIntegration.jl",
     authors = "Shahriar Iravanian",
     modules = [SymbolicNumericIntegration],
     clean = true, doctest = false, linkcheck = true,
-    warnonly = [:missing_docs],
+    # warnonly = [:missing_docs],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/SymbolicNumericIntegration/stable/"),
     pages = pages)

docs/src/index.md

@@ -22,52 +22,52 @@ Pkg.add("SymbolicNumericIntegration")
 Examples:
 
 ```julia
-julia> using SymbolicNumericIntegration
-julia> using Symbolics
+using Symbolics
+using SymbolicNumericIntegration
 
-julia> @variables x a b
+@variables x a b
 
 # if `detailed = true` (default), the output is a tuple of (solution, unsolved portion, err)
 
-julia> integrate(3x^3 + 2x - 5)
+integrate(3x^3 + 2x - 5)
 (x^2 + (3//4)*(x^4) - (5x), 0, 0)
 
-julia> integrate((5 + 2x)^-1)
+integrate((5 + 2x)^-1)
 ((1//2)*log((5//2) + x), 0, 0.0)
 
 # `detailed = false` simplifies the output to just the resulting integral
 
-julia> integrate(x^2 / (16 + x^2); detailed = false)
+integrate(x^2 / (16 + x^2); detailed = false)
 x + 4atan((-1//4)*x)
 
-julia> integrate(x^2 * log(x); detailed = false)
+integrate(x^2 * log(x); detailed = false)
 (1//3)*(x^3)*log(x) - (1//9)*(x^3)
 
-julia> integrate(sec(x) * tan(x); detailed = false)
+integrate(sec(x) * tan(x); detailed = false)
 sec(x)
 
 # 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(sin(a * x), x; detailed = false, symbolic = true)
+integrate(sin(a * x), x; detailed = false, symbolic = true)
 (-cos(a*x)) / a
 
-julia> integrate(x^2 * cos(a * x), x; detailed = false, symbolic = true)
+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(log(log(a * x)) / x, x; detailed = false, symbolic = true)
+integrate(log(log(a * x)) / x, x; detailed = false, symbolic = true)
 log(a*x)*log(log(a*x)) - log(a*x)
 
 # multiple symbolic constants
 
-julia> integrate(cosh(a * x) * exp(b * x), x; detailed = false, symbolic = true)
+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))
 
 # definite integration, passing a tuple of (x, lower bound, higher bound) in the 
 # second argument
 
-julia> integrate(x * sin(a * x), (x, 0, 1); symbolic = true, detailed = false)
+integrate(x * sin(a * x), (x, 0, 1); symbolic = true, detailed = false)
 (sin(a) - a*cos(a)) / (a^2)
 ```
 
@@ -81,13 +81,13 @@ SymbolicNumericIntegration.jl exports some special integral functions (defined o
 For examples:
 
 ```julia
-julia> integrate(exp(x + 1) / (x + 1))
+integrate(exp(x + 1) / (x + 1))
 (SymbolicNumericIntegration.Ei(1 + x), 0, 0.0)
 
-julia> integrate(x * cos(a*x^2 - 1) / (a*x^2 - 1), x; detailed=false, symbolic=true)
+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))), x; detailed=false, symbolic=true)
+integrate(1 / (x * log(log(x))), x; detailed = false, symbolic = true)
 SymbolicNumericIntegration.Li(log(x))
 ```
 
@@ -213,4 +213,4 @@ file and the
 [project]($link_project)
 file.
 """)
-```
+```

src/SymbolicNumericIntegration.jl

@@ -11,8 +11,8 @@ using DataDrivenDiffEq, DataDrivenSparse
 
 struct NumericalPlan
     abstol::Float64
-    radius::Float64    
-    complex_plane::Bool    
+    radius::Float64
+    complex_plane::Bool
     opt::DataDrivenDiffEq.AbstractDataDrivenAlgorithm
 end
 
@@ -36,7 +36,7 @@ include("sparse.jl")
 # include("optim.jl")
 include("integral.jl")
 
-export integrate, generate_basis
+export integrate, generate_basis, best_hints
 
 include("symbolic.jl")
 

src/candidates.jl

@@ -1,6 +1,6 @@
 using DataStructures
 
-# this is the main heurisctic used to find the test fragments
+# this is the old heurisctic used to find the test fragments
 function generate_basis(eq, x, try_kernel = true)
     if !try_kernel
         S = sum(generate_homotopy(expr(eq), x))

src/homotopy.jl

@@ -80,6 +80,9 @@ end
 
 guard_zero(x) = isequal(x, 0) ? one(x) : x
 
+# The core of ansatz generation.
+# The name is not accurate and should be changed (not really homotopy, 
+# just inspired from!)
 function generate_homotopy(eq, x)
     eq = value(eq)
     x = value(x)
@@ -95,7 +98,7 @@ function generate_homotopy(eq, x)
     for i in 1:length(ks)
         μ = u[i]
         y, dy = apply_partial_int_rules(sub[μ], x)
-        
+
         y = substitute(y, sub)
         ∂y = guard_zero(diff(dy, x))
 
@@ -157,8 +160,9 @@ partial_int_rules = [
     @rule 𝛷(~x, asech(~u)) => (asech(~u), ~u)
     @rule 𝛷(~x, acoth(~u)) => (~u * acot(~u) + log(~u + 1), ~u)
 # logarithmic and exponential functions
-    @rule 𝛷(~x, log(~u)) => (~u + ~u * log(~u) + sum(pow_minus_rule(~u, ~x, -1); init = one(~u)),
-    ~u);
+    @rule 𝛷(~x, log(~u)) => (~u + ~u * log(~u) +
+                             sum(pow_minus_rule(~u, ~x, -1); init = one(~u)),
+    ~u)
     @rule 𝛷(~x, 1 / log(~u)) => (log(log(~u)) + li(~u), ~u)
     @rule 𝛷(~x, exp(~u)) => (exp(~u) + ei(~u) + erfi_(~x), ~u)
     @rule 𝛷(~x, ^(exp(~u), ~k::is_neg)) => (^(exp(-~u), -~k), ~u)
@@ -167,10 +171,13 @@ partial_int_rules = [
     @rule 𝛷(~x, sqrt(~u)) => (sum(sqrt_rule(~u, ~x, 0.5); init = one(~u)), ~u);
     @rule 𝛷(~x, 1 / sqrt(~u)) => (sum(sqrt_rule(~u, ~x, -0.5); init = one(~u)), ~u);
 # rational functions                                                              
-    @rule 𝛷(~x, 1 / ^(~u::is_univar_poly, ~k::is_pos_int)) => (sum(pow_minus_rule(~u, ~x, -~k);
+    @rule 𝛷(~x, 1 / ^(~u::is_univar_poly, ~k::is_pos_int)) => (sum(pow_minus_rule(~u,
+            ~x,
+            -~k);
         init = one(~u)),
+    ~u)
+    @rule 𝛷(~x, 1 / ~u::is_univar_poly) => (sum(pow_minus_rule(~u, ~x, -1); init = one(~u)),
     ~u);
-    @rule 𝛷(~x, 1 / ~u::is_univar_poly) => (sum(pow_minus_rule(~u, ~x, -1); init = one(~u)), ~u);
     @rule 𝛷(~x, ^(~u, -1)) => (log(~u) + ~u * log(~u), ~u)
     @rule 𝛷(~x, ^(~u, ~k::is_neg_int)) => (sum(^(~u, i) for i in (~k + 1):-1), ~u)
     @rule 𝛷(~x, 1 / ~u) => (log(~u), ~u)
@@ -187,7 +194,7 @@ end
 
 function pow_minus_rule(p, x, k; abstol = 1e-8)
     if !is_univar_poly(p)
-        return [p^k, p^(k + 1), log(p), p*log(p)]
+        return [p^k, p^(k + 1), log(p), p * log(p)]
     end
 
     # x = var(p)
@@ -203,7 +210,7 @@ function pow_minus_rule(p, x, k; abstol = 1e-8)
     r = nice_parameter.(r)
     s = nice_parameter.(s)
 
-    # ∫ 1 / ((x-z₁)(x-z₂)) dx = ... + c₁ * log(x-z₁) + c₂ * log(x-z₂)
+    # applying ∫ 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]
@@ -225,9 +232,9 @@ end
 
 function sqrt_rule(p, x, k)
     h = Any[p^k, p^(k + 1)]
-    
+
     Δ = diff(p, x)
-    push!(h, log(Δ/2 + sqrt(p)))
+    push!(h, log(Δ / 2 + sqrt(p)))
 
     if !is_univar_poly(p)
         return h
@@ -254,4 +261,3 @@ function sqrt_rule(p, x, k)
 
     return h
 end
-