fixing homotopy

Author: shahriariravanian

src/candidates.jl

@@ -2,10 +2,6 @@ using DataStructures
 
 # this is the main heurisctic used to find the test fragments
 function generate_basis(eq, x; homotopy=false)
-    if homotopy
-         return generate_homotopy(eq, x)
-    end
-
     eq = expand(eq)
     S = Set{Any}()
     for t in terms(eq)

src/homotopy.jl

@@ -22,17 +22,23 @@ transformer(eq::SymbolicUtils.Div, f) = transformer(arguments(eq)[1],f) * transf
 function transformer(eq::SymbolicUtils.Pow, f)
     y, k = arguments(eq)
 
-    r = nice_parameter(k)
-    if r isa Rational || isinteger(r)
-        a, b = numerator(r), denominator(r)
+    # if is_poly(y)
+    #     return next_variable!(f, y)^k
+    # end
+
+    # r = nice_parameter(k)
+    # if r isa Rational || isinteger(r)
+    if isinteger(k)
+        a, b = k, 1
+        # a, b = numerator(r), denominator(r)
         if k < 0
             y = inv(y)
         end
         f.deg = max(f.deg, abs(a))
         μ = next_variable!(f, b == 1 ?  y : y ^(1/b))
         return μ ^ abs(a)
     else
-        return next_variable!(f, y ^ k)
+        return next_variable!(f, y^k)
     end
 end
 
@@ -48,7 +54,11 @@ end
 function transform(eq, x)
     eq = substitute(eq, Dict(x => 𝑥))
     f = Transform(1, Dict(), 1, false)
-    return transformer(eq, f), f.sub, f.deg
+    p = transformer(eq, f)
+    if !any(is_poly, values(f.sub))
+        p *= next_variable!(f, 1)
+    end
+    return p, f.sub, f.deg
 end
 
 function homotopy_integrand(eq, x)
@@ -68,30 +78,30 @@ function homotopy_integrand(eq, x)
 end
 
 function expand_integrand(I, x, deg)
-    E = sum((Differential(x)^i)(I) for i=1:deg-1; init=I) #* (1+x)
-    S = Set{Any}()
-    enqueue_expr_ex!(S, expand(expand_derivatives(E)), x)
-    return [one(x); [s for s in S]]
-
+    # E = sum((Differential(x)^i)(I) for i=1:deg-1; init=I) #* (1+x)
     # S = Set{Any}()
-    # # T = Set{Any}()
-    # Q₁ = Queue{Any}()
-    #
-    # enqueue_expr_ex!(S, Q₁, expand(I), x)
-    #
-    # D = Differential(x)
-    #
-    # for i = 1:deg
-    #     Q₂ = Queue{Any}()
-    #     while !isempty(Q₁) # && length(S) < max_terms
-    #         y = dequeue!(Q₁)
-    #         E = expand(expand_derivatives(D(y)))
-    #         enqueue_expr_ex!(S, Q₂, E, x)
-    #     end
-    #     Q₁ = Q₂
-    # end
-    #
+    # enqueue_expr_ex!(S, expand(expand_derivatives(E)), x)
     # return [one(x); [s for s in S]]
+
+    S = Set{Any}()
+    # T = Set{Any}()
+    Q₁ = Queue{Any}()
+
+    enqueue_expr_ex!(S, Q₁, expand(I), x)
+
+    D = Differential(x)
+
+    for i = 1:deg
+        Q₂ = Queue{Any}()
+        while !isempty(Q₁) # && length(S) < max_terms
+            y = dequeue!(Q₁)
+            E = expand(expand_derivatives(D(y)))
+            enqueue_expr_ex!(S, Q₂, E, x)
+        end
+        Q₁ = Q₂
+    end
+
+    return [one(x); [s for s in S]]
 end
 
 function expand_integrand(I, x, deg)
@@ -113,7 +123,8 @@ function substitute_x(eq, x, sub)
 end
 
 function generate_homotopy2(eq, x)
-    q, sub, deg = transform(eq, x)
+    q, sub = transform(eq, x)
+    d = degree(q)
     n = length(sub)
 
     S = Set{Any}()
@@ -126,17 +137,17 @@ function generate_homotopy2(eq, x)
         h₁ = substitute_x(h₁, x, sub)
         h₂ = substitute_x(h₂, x, sub)
 
-        H = sum((Differential(x)^i)(h₁ * h₂) for i=1:deg-1; init=h₁*h₂)
-        I = expand(expand_derivatives(H))
+        H = sum((Differential(x)^i)(h₂) for i=1:d-1; init=(1 + h₂))
+        I = expand(expand_derivatives((1 + h₁) * H))
         enqueue_expr_ex!(S, I, x)
     end
 
-    H = sum((Differential(x)^i)(eq) for i=1:deg-1; init=eq)
-    I = expand((1+x) * expand_derivatives(H))
+    # H = sum((Differential(x)^i)(eq) for i=1:deg-1; init=eq)
+    # I = expand((1+x) * expand_derivatives(H))
     # enqueue_expr_ex!(S, I, x)
 
     # return [one(x); [s for s in S]]
-    return [[x^i for i=0:deg]; [s for s in S]]
+    return [one(x); [s for s in S]]
 end
 
 ##############################################################################

src/integral.jl

@@ -22,7 +22,7 @@ Base.signbit(z::Complex{T}) where T<:Number = signbit(real(z))
     a pair of expressions, solved is the solved integral and unsolved is the residual unsolved
     portion of the input
 """
-function integrate(eq, x=nothing; abstol=1e-6, num_steps=3, num_trials=5, radius=1.0,
+function integrate(eq, x=nothing; abstol=1e-6, num_steps=2, num_trials=5, radius=1.0,
                    show_basis=false, opt = STLSQ(exp.(-10:1:0)), bypass=false,
                    symbolic=true, bypart=true, max_basis=100,
                    verbose=false, complex_plane=true, prune_basis=false, homotopy=false)
@@ -52,6 +52,8 @@ function integrate(eq, x=nothing; abstol=1e-6, num_steps=3, num_trials=5, radius
     #     eq = q
     # end
 
+    homotopy = homotopy && !bypass
+
     s₁, u₁, ϵ = integrate_sum(eq, x, l; bypass, abstol, num_trials, num_steps,
                               radius, show_basis, opt, symbolic,
                               max_basis, verbose, complex_plane, prune_basis, homotopy)

src/rules.jl

@@ -178,8 +178,6 @@ function decompose_rational(eq)
     end
     F = unique(F)
 
-    println(F)
-
     n = length(F)
     A = zeros(Complex, (n,n))
     b = zeros(Complex, n)

test/runtests.jl

@@ -190,19 +190,20 @@ basic_integrals = [
 ]
 
 function test_integrals(; kw...)
-    bypass = false
+    args = Dict(kw)
     misses = []
     k = 1
 
     for eq in basic_integrals
         if isequal(eq, β)
             printstyled("**** bypass on ****\n"; color=:red)
             bypass = true
+            args[:bypass] = true
         else
             printstyled(k, ": "; color=:blue)
             k += 1
             printstyled(eq, " =>\n"; color=:green)
-            solved, unsolved = integrate(eq; kw...)
+            solved, unsolved = integrate(eq; args...)
             printstyled('\t', solved; color=:cyan)
             if isequal(unsolved, 0)
                 println()