symbolic integration updated

Author: shahriariravanian

src/candidates.jl

@@ -65,7 +65,7 @@ function closure(eq, x; max_terms = 50)
 end
 
 function candidate_pow_minus(p, k; abstol = 1e-8)
-    if isnan(poly_deg(p))
+    if !is_univar_poly(p)
         return [p^k, p^(k + 1), log(p)]
     end
 
@@ -104,10 +104,14 @@ function candidate_pow_minus(p, k; abstol = 1e-8)
 end
 
 function candidate_sqrt(p, k)
-    x = var(p)
-
     h = Any[p^k, p^(k + 1)]
-
+    
+    if !is_univar_poly(p)
+        return h
+    end    
+    
+    x = var(p)
+    
     if poly_deg(p) == 2
         r, s = find_roots(p, x)
         l = leading(p, x)

src/homotopy.jl

@@ -64,11 +64,14 @@ Symbolics.@register_symbolic Ei(z)
 Symbolics.@register_symbolic Si(z)
 Symbolics.@register_symbolic Ci(z)
 Symbolics.@register_symbolic Li(z)
+Symbolics.@register_symbolic erfi(z)
+
 
 Symbolics.derivative(::typeof(Ei), args::NTuple{1, Any}, ::Val{1}) = exp(args[1]) / args[1]
 Symbolics.derivative(::typeof(Si), args::NTuple{1, Any}, ::Val{1}) = sin(args[1]) / args[1]
 Symbolics.derivative(::typeof(Ci), args::NTuple{1, Any}, ::Val{1}) = cos(args[1]) / args[1]
 Symbolics.derivative(::typeof(Li), args::NTuple{1, Any}, ::Val{1}) = 1 / log(args[1])
+Symbolics.derivative(::typeof(erfi), args::NTuple{1, Any}, ::Val{1}) = 2 / sqrt(2) * exp(args[1]^2)
 
 @syms 𝑥 si(𝑥) ci(𝑥) ei(𝑥) li(𝑥)
 
@@ -77,8 +80,12 @@ Symbolics.derivative(::typeof(Li), args::NTuple{1, Any}, ::Val{1}) = 1 / log(arg
 guard_zero(x) = isequal(x, 0) ? one(x) : x
 
 function generate_homotopy(eq, x)
-    eq = eq isa Num ? eq.val : eq
-    x = x isa Num ? x.val : x
+    eq = value(eq)
+    x = value(x)
+    
+    if is_add(eq)
+        return unique(union([generate_homotopy(t,x) for t in args(eq)]...))
+    end
 
     p = transform(eq, x)
     q, sub, ks = rename_factors(p, (si => Si, ci => Ci, ei => Ei, li => Li))
@@ -107,7 +114,7 @@ partial_int_rules = [
                      # trigonometric functions
                      @rule 𝛷(sin(~x)) => (cos(~x) + si(~x), ~x)
                      @rule 𝛷(cos(~x)) => (sin(~x) + ci(~x), ~x)
-                     @rule 𝛷(tan(~x)) => (log(cos(~x)), ~x)
+                     @rule 𝛷(tan(~x)) => (log(cos(~x)) + ~x*tan(~x), ~x)
                      @rule 𝛷(csc(~x)) => (log(csc(~x) + cot(~x)) + log(sin(~x)), ~x)
                      @rule 𝛷(sec(~x)) => (log(sec(~x) + tan(~x)) + log(cos(~x)), ~x)
                      @rule 𝛷(cot(~x)) => (log(sin(~x)), ~x)

src/integral.jl

@@ -29,7 +29,7 @@ Keyword Arguments:
 - `num_trials` (default: `10`): the number of trials in each step (no changes to the basis)
 - `show_basis` (default: `false`): if true, the basis (list of candidate terms) is printed
 - `bypass` (default: `false`): if true do not integrate terms separately but consider all at once
-- `symbolic` (default: `false`): try symbolic integration first
+- `symbolic` (default: `false`): try symbolic integration first (will be forced if `eq` has constant parameters)
 - `max_basis` (default: `100`): the maximum number of candidate terms to consider
 - `verbose` (default: `false`): print a detailed report
 - `complex_plane` (default: `true`): generate random test points on the complex plane (if false, the points will be on real axis)
@@ -47,7 +47,7 @@ Output:
 function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 10,
                    radius = 1.0,
                    show_basis = false, opt = STLSQ(exp.(-10:1:0)), bypass = false,
-                   symbolic = true, max_basis = 100, verbose = false, complex_plane = true,
+                   symbolic = false, max_basis = 100, verbose = false, complex_plane = true,
                    homotopy = true, use_optim = false)
     eq = expand(eq)
 
@@ -165,6 +165,25 @@ function integrate_term(eq, x, l; kwargs...)
         result(l, "Successful", y)
         return y, 0, 0
     end
+    
+    params = const_params(eq, x)
+    
+    if !isempty(params) && !symbolic
+        @warn("The input expression has constant parameters: [$(join(params, ", "))], forcing `symbolic = true`")
+        symbolic = true
+    end
+    
+    if symbolic
+        y = integrate_symbolic(eq, x; abstol, radius)
+        if y == nothing
+            if !isempty(params)
+                @warn("Symbolic integration failed. Try changing constant parameters ([$(join(params, ", "))]) to numerical values.")
+                return 0, eq, Inf
+            end
+        else
+            return y, 0, 0
+        end
+    end
 
     eq = cache(eq)
     basis1 = generate_basis(eq, x, true)
@@ -190,19 +209,7 @@ function integrate_term(eq, x, l; kwargs...)
     for i in 1:num_steps
         if length(basis1) > max_basis
             break
-        end
-
-        if symbolic
-            y, ϵ = try_symbolic(Float64, expr(eq), x, expr.(basis1), deriv!.(basis1, x);
-                                kwargs...)
-
-            if !isequal(y, 0) && accept_solution(eq, x, y, radius) < abstol
-                result(l, "Successful symbolic", y)
-                return y, 0, 0
-            else
-                inform(l, "Failed symbolic")
-            end
-        end
+        end        
 
         for j in 1:num_trials
             basis = isodd(j) ? basis1 : basis2

src/sparse.jl

@@ -2,7 +2,7 @@
 function solve_sparse(eq, x, basis, radius; kwargs...)
     args = Dict(kwargs)
     abstol, opt, complex_plane = args[:abstol], args[:opt], args[:complex_plane]
-
+    
     A, X = init_basis_matrix(eq, x, basis, radius, complex_plane; abstol)
 
     # find a linearly independent subset of the basis
@@ -52,8 +52,14 @@ function init_basis_matrix(eq, x, basis, radius, complex_plane; abstol = 1e-6)
 
     # A is an nxn matrix holding the values of the fragments at n random points
     A = zeros(Complex{Float64}, (n, n))
-    X = zeros(Complex{Float64}, n)
-
+    X = zeros(Complex{Float64}, n)    
+    
+    S = subs_symbols(eq, x)
+    if !isempty(S)
+        eq = substitute(eq, S)
+        basis = [substitute(y, S) for y in basis]
+    end
+    
     eq_fun = fun!(eq, x)
     Δbasis_fun = deriv_fun!.(basis, x)
 
@@ -162,3 +168,54 @@ function find_dense(A, basis; abstol = 1e-6)
     end
     return nothing, Inf
 end
+
+
+###############################################################################
+
+
+function hints(eq, x, basis; radius=5.0, abstol=1e-6, opt=STLSQ(exp.(-10:1:0)), complex_plane=true)
+    A, X = init_basis_matrix(eq, x, basis, radius, complex_plane; abstol)
+
+    # find a linearly independent subset of the basis
+    l = find_independent_subset(A; abstol)
+    A, basis = A[l, l], basis[l]
+    
+    n, m = size(A)
+
+    try
+        b = ones((1, n))
+        solver = SparseLinearSolver(opt,
+                                    options = DataDrivenCommonOptions(verbose = false,
+                                                                      maxiters = 1000))
+        res, _... = solver(A', b)
+        q = DataDrivenSparse.coef(first(res))
+        err = abs(rms(A * q' .- 1)) 
+        if err < abstol
+            sel = abs.(q) .> abstol
+            h = [basis[i] for i in 1:length(basis) if sel[i]]
+        else 
+            h = []
+        end
+        return h, err
+    catch e
+        println("Error from hints", e)        
+    end
+    
+    return 0, Inf
+end
+
+
+function best_hints(eq, x, basis; radius=5.0, abstol=1e-6, opt=STLSQ(exp.(-10:1:0)), complex_plane=true, num_trials=10)        
+    H = []
+    L = Int[]
+    
+    for _ in 1:num_trials
+        h, err = hints(eq, x, basis; radius, abstol, opt, complex_plane)
+        push!(H, h)
+        push!(L, err < abstol ? length(h) : length(basis))
+    end
+    
+    _, idx = findmin(L)
+    return H[idx]
+end
+