homotomy expansion fixed

Author: shahriariravanian

Project.toml

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

src/candidates.jl

@@ -4,9 +4,9 @@ using DataStructures
 function generate_basis(eq, x; homotopy=false)
     # if homotopy return generate_homotopy2(eq, x) end
     eq = expand(eq)
-    S = Set{Any}()
+    S = 0 #Set{Any}()
     for t in terms(eq)
-        q = t / coef(t, x)
+        q = equivalent(t, x)
         f = kernel(q)
 
         C₁ = closure(f, x)
@@ -22,30 +22,18 @@ function generate_basis(eq, x; homotopy=false)
             C₂ = find_candidates(p, x)
         end
 
-        for c₁ in C₁
-            enqueue_expr!(S, c₁, x)
-        end
-
-        for c₂ in C₂
-            enqueue_expr!(S, c₂, x)
-        end
-
-        for c₁ in C₁
-            for c₂ in C₂
-                enqueue_expr!(S, expand(c₁*c₂), x)
-            end
-        end
+        S += sum(c₁*c₂ for c₁ in C₁ for c₂ in C₂)
     end
-    return unique([one(x); [s for s in S]])
+    return unique([one(x); [equivalent(t,x) for t in terms(S)]])
 end
 
 function expand_basis(basis, x; homotopy=false)
-    if homotopy
-        I, deg = homotopy_integrand(sum(basis), x)
-        return expand_integrand(I, x, deg)
-    else
-        return unique([basis; basis*x])
-    end
+    b = sum(basis)
+    eq = (1 + x) * (b + Differential(x)(b))
+    eq = expand(expand_derivatives(eq))
+    S = Set{Any}()
+    enqueue_expr!(S, eq, x)
+    return [one(x); [s for s in S]]
 end
 
 function closure(eq, x; max_terms=50)
@@ -84,7 +72,12 @@ function candidate_pow_minus(p, k; abstol=1e-8)
     end
 
     x = var(p)
-    r, s = find_roots(p, x)
+    d = poly_deg(p)
+
+    for j = 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)

src/homotopy.jl

@@ -20,9 +20,12 @@ transformer(eq::SymbolicUtils.Div, f) = transformer(arguments(eq)[1],f) * transf
 function transformer(eq::SymbolicUtils.Pow, f)
     y, k = arguments(eq)
 
-    if isinteger(k)
-        μ = next_variable!(f, k < 0 ? inv(y) : y)
-        return μ ^ abs(k)
+    if is_pos_int(k)
+        μ = next_variable!(f, y)
+        return μ ^ k
+    elseif is_neg_int(k)
+        μ = next_variable!(f, inv(y))
+        return μ ^ -k
     else
         return next_variable!(f, y^k)
     end
@@ -54,32 +57,29 @@ function substitute_x(eq, x, sub)
 end
 
 function generate_homotopy(eq, x)
-    q, sub = transform(eq, x)
-    d = degree(q)
-    n = length(sub)
+    q, sub = transform(eq, x)    
+    S = 0
 
-    S = Set{Any}()
-
-    for i = 1:n
+    for i = 1:length(sub)
         μ = u[i]
         h₁, ∂h₁ = apply_partial_int_rules(sub[μ])
         h₂ = expand_derivatives(Differential(μ)(q))
 
         h₁ = substitute_x(h₁, x, sub)
         h₂ = substitute_x(h₂ * ∂h₁^-1, x, sub)
 
-        H = sum((Differential(x)^i)(h₂) for i=1:d-1; init=(1 + h₂))
-        I = expand(expand_derivatives((1 + h₁) * H))
-        enqueue_expr!(S, I, x)
+        S += expand((1 + h₁) * (1 + h₂))
     end
 
-    return [one(x); [s for s in S]]
+    unique([one(x); [equivalent(t,x) for t in terms(S)]])
 end
 
 ##############################################################################
 
-∂(x) = expand_derivatives(Differential(𝑥)(x))
-cabs(x) = sqrt(x * conj(x))
+function ∂(x)
+    d = expand_derivatives(Differential(𝑥)(x))
+    return isequal(d, 0) ? 1 : d
+end
 
 partial_int_rules = [
     @rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
@@ -131,7 +131,6 @@ partial_int_rules = [
     @rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x,0.5); init=one(~x)), 1)
     @rule 𝛷(^(sqrt(~x),-1)) => 𝛷(^(~x,-0.5))
 
-
     @rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
     @rule 𝛷(1 / ~x) => 𝛷(^(~x, -1))
     @rule 𝛷(^(~x, ~k)) => (^(~x, ~k+1), ∂(~x))

src/integral.jl

@@ -93,7 +93,7 @@ function integrate_sum(eq, x, l; bypass=false, kwargs...)
         eq = apply_q_rules(apply_integration_rules(unsolved))
 
         if !isequal(eq, unsolved)
-            # eq = expand(eq)
+            eq = expand(eq)
             unsolved = 0
             ϵ₀ = 0
             ts = bypass ? [eq] : terms(eq)
@@ -169,7 +169,7 @@ function integrate_term(eq, x, l; kwargs...)
         inform(l, "Generating basis (|β| = $(length(basis)))", basis)
     end
 
-    if prune_basis || length(basis) > max_basis
+    if prune_basis # || length(basis) > max_basis
         basis, ok = prune(basis, eq, x)
         if ok && show_basis
             inform(l, "Prunning the basis (|β| = $(length(basis)))", basis)

src/roots.jl

@@ -33,7 +33,7 @@ function find_roots(T, p, x; abstol=1e-8, num_roots=0)
     s = Complex{T}[]
 
     while !isequal(p, 0)
-        q = expand(p / x)
+        q = expand(p * x^-1)
         if !is_poly(q) break end
         push!(r, 0)
         p = q
@@ -45,7 +45,7 @@ function find_roots(T, p, x; abstol=1e-8, num_roots=0)
     ∂p = expand_derivatives(Differential(x)(p))
 
     while length(zs) < n
-        z = solve_newton(Complex{T}, p, ∂p, x, cis(2π*rand()), zs; abstol)        
+        z = solve_newton(Complex{T}, p, ∂p, x, cis(2π*rand()), zs; abstol)
         if z != nothing
             if abs(imag(z)) < abstol
                 push!(zs, Complex(real(z)))

src/rules.jl

@@ -73,6 +73,7 @@ is_linear_poly(eq) = poly_deg(eq) == 1
 s_rules = [
     @rule Ω(+(~~xs)) => sum(map(Ω, ~~xs))
     @rule Ω(*(~~xs)) => prod(map(Ω, ~~xs))
+    @rule Ω(~x / ~y) => Ω(~x) * Ω(^(~y,-1))
     @rule Ω(^(~x::is_linear_poly, ~k::is_pos_int)) => ^(~x, ~k)
     @rule Ω(~x::is_var) => ~x
     @rule Ω(~x::is_linear_poly) => ~x
@@ -95,7 +96,9 @@ kernel(eq) = Prewalk(Chain(s_rules))(Ω(value(eq)))
 
 coef(eq::SymbolicUtils.Mul, x) = prod(t for t in arguments(eq) if !isdependent(t,x); init=1)
 coef(eq::SymbolicUtils.Add, x) = minimum(abs(coef(t,x)) for t in arguments(eq))
-coef(eq, x) = 1
+coef(eq, x) = is_number(eq) ? eq : 1
+
+equivalent(eq, x) = eq / coef(eq, x)
 
 ########################## Transformation Rules ###############################
 
@@ -173,7 +176,14 @@ end
 function decompose_rational(eq)
     if poly_deg(eq) == 1 return inverse(eq) end
     x = var(eq)
-    r, s = find_roots(eq, x)
+
+    d = poly_deg(eq)
+
+    for j = 1:10  # will try 10 times to find the roots
+        r, s = find_roots(eq, x)
+        if length(r) + length(s) >= d break end
+    end
+
     s = s[1:2:end]
     r = nice_parameter.(r)
     s = nice_parameter.(s)

test/axiom.jl

@@ -6,7 +6,7 @@ using Symbolics
 using PyCall
 sympy = pyimport("sympy")
 
-@syms 𝐞 a b c d e p
+@syms x 𝐞 a b c d e p t m n z
 
 axion_rules = [
     @rule ^(𝐞, ~k) => exp(~k)
@@ -26,8 +26,10 @@ function convert_axiom(name::AbstractString)
         ma = match(re, line)
 
         if ma != nothing
-            line = replace(ma.match, "%e" => 𝐞)
-            line  = replace(ma.match, "%i" => im)
+            line = ma.match
+            println(line)
+            line = replace(line, "%e" => 𝐞)
+            line  = replace(line, "%i" => im)
 
             try
                 expr = Meta.parse(line)
@@ -49,6 +51,8 @@ function convert_axiom(name::AbstractString)
                 println(P[1])
                 push!(L, P)
             catch err
+                printstyled(lineno, ": ", err, '\n'; color=:red)
+                printstyled('\t', line, '\n'; color=:red)
             end
         end
     end
@@ -117,8 +121,9 @@ function test_axiom(L, try_sympy=true; kwargs...)
     n_diff = 0
     n_catch = 0
 
-    for p in L
+    for (i,p) in enumerate(L)
         try
+            printstyled(i, ": "; color=:yellow)
             eq = p[1]
             x = p[2]
             ans = p[4]
@@ -145,7 +150,7 @@ function test_axiom(L, try_sympy=true; kwargs...)
 
             printstyled("\tAnswer: \t", ans, '\n'; color=:blue)
         catch e
-            printstyled(e, '\n'; color=:red)
+            printstyled(i, ": ", e, '\n'; color=:red)
             n_catch += 1
         end
     end

test/results.md

@@ -223,6 +223,6 @@ function test_integrals(; kw...)
 end
 
 @testset "integral" begin test_integrals(;
-        symbolic=false, bypart=false, prune_basis=false, verbose=false,
-        homotopy=true, num_steps=2, num_trials=10)
+    symbolic=false, bypart=false, prune_basis=false, verbose=false,
+    homotopy=true, num_steps=2, num_trials=10)
 end