homotopy rules expanded

Author: shahriariravanian

src/candidates.jl

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

src/homotopy.jl

@@ -138,26 +138,38 @@ function generate_homotopy2(eq, x)
         h₂ = substitute_x(h₂, x, sub)
 
         H = sum((Differential(x)^i)(h₂) for i=1:d-1; init=(1 + h₂))
-        I = expand(expand_derivatives((1 + h₁) * H))
+        I = expand(expand_derivatives((1 + h₁) * H))        
+        # H = sum((Differential(x)^i)(h₁*h₂) for i=1:d-1; init=(h₁ + h₂ + h₁*h₂))
+        # I = expand(expand_derivatives(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))
-    # enqueue_expr_ex!(S, I, x)
-
-    # return [one(x); [s for s in S]]
     return [one(x); [s for s in S]]
 end
 
 ##############################################################################
 
 ∂(x) = expand_derivatives(Differential(𝑥)(x))
+cabs(x) = sqrt(x * conj(x))
 
 H_rules = [
     # @rule 𝛷(+(~~xs)) => sum(map(𝛷, ~~xs))
     # @rule 𝛷(*(~~xs)) => prod(map(𝛷, ~~xs))
 
+    @rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
+    @rule 𝛷(^(cos(~x), ~k::is_neg)) => 𝛷(^(sec(~x), -~k))
+    @rule 𝛷(^(tan(~x), ~k::is_neg)) => 𝛷(^(cot(~x), -~k))
+    @rule 𝛷(^(csc(~x), ~k::is_neg)) => 𝛷(^(sin(~x), -~k))
+    @rule 𝛷(^(sec(~x), ~k::is_neg)) => 𝛷(^(cos(~x), -~k))
+    @rule 𝛷(^(cot(~x), ~k::is_neg)) => 𝛷(^(tan(~x), -~k))
+
+    @rule 𝛷(^(sinh(~x), ~k::is_neg)) => 𝛷(^(csch(~x), -~k))
+    @rule 𝛷(^(cosh(~x), ~k::is_neg)) => 𝛷(^(sech(~x), -~k))
+    @rule 𝛷(^(tanh(~x), ~k::is_neg)) => 𝛷(^(coth(~x), -~k))
+    @rule 𝛷(^(csch(~x), ~k::is_neg)) => 𝛷(^(sinh(~x), -~k))
+    @rule 𝛷(^(sech(~x), ~k::is_neg)) => 𝛷(^(cosh(~x), -~k))
+    @rule 𝛷(^(coth(~x), ~k::is_neg)) => 𝛷(^(tanh(~x), -~k))
+
     @rule 𝛷(sin(~x)) => cos(~x) * ∂(~x)^-1
     @rule 𝛷(cos(~x)) => sin(~x) * ∂(~x)^-1
     @rule 𝛷(tan(~x)) => log(cos(~x)) * ∂(~x)^-1
@@ -175,8 +187,8 @@ H_rules = [
     @rule 𝛷(asin(~x)) => (~x*asin(~x) + sqrt(1 - ~x*~x)) * ∂(~x)^-1
     @rule 𝛷(acos(~x)) => (~x*acos(~x) + sqrt(1 - ~x*~x)) * ∂(~x)^-1
     @rule 𝛷(atan(~x)) => (~x*atan(~x) + log(~x*~x + 1)) * ∂(~x)^-1
-    @rule 𝛷(acsc(~x)) => acsc(~x) * ∂(~x)^-1
-    @rule 𝛷(asec(~x)) => asec(~x) * ∂(~x)^-1
+    @rule 𝛷(acsc(~x)) => (~x*acsc(~x) + acosh(cabs(~x))) * ∂(~x)^-1
+    @rule 𝛷(asec(~x)) => (~x*asec(~x) + acosh(cabs(~x))) * ∂(~x)^-1
     @rule 𝛷(acot(~x)) => (~x*acot(~x) + log(~x*~x + 1)) * ∂(~x)^-1
 
     @rule 𝛷(asinh(~x)) => (~x*asinh(~x) + sqrt(~x*~x + 1)) * ∂(~x)^-1

src/integral.jl

@@ -351,8 +351,7 @@ function try_integrate(T, eq, x, basis, Δbasis, radius; kwargs...)
 end
 
 function init_basis_matrix!(T, A, X, x, eq, Δbasis, radius, complex_plane; abstol=1e-6)
-    n = size(A, 1)
-    # X = zeros(Complex{T}, n)
+    n = size(A, 1)    
     k = 1
     i = 1
 
@@ -372,7 +371,7 @@ function init_basis_matrix!(T, A, X, x, eq, Δbasis, radius, complex_plane; abst
                 end
             end
         catch e
-            println(e)
+            println("Error from init_basis_matrix!: ", e)
         end
     end
 end
@@ -408,7 +407,7 @@ function sparse_fit(T, A, x, basis, Δbasis, opt; abstol=1e-6)
         q = nice_parameter.(q₀)
         return sum(q[i]*basis[i] for i = 1:length(basis) if q[i] != 0; init=zero(x)), abs(ϵ)
     catch e
-        println(e)
+        println("Error from spart_fit", e)
         return nothing, Inf
     end
 end

src/rules.jl

@@ -170,7 +170,7 @@ function decompose_rational(eq)
     r = nice_parameter.(r)
     s = nice_parameter.(s)
 
-    F = [(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s] ∪
+    F = Any[(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s] ∪
         [x*(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s]
     for i in eachindex(r)
         μ = sum(r[1:i] .== r[i])

test/axiom.jl

@@ -5,34 +5,52 @@ using SymbolicUtils.Rewriters
 using PyCall
 sympy = pyimport("sympy")
 
-@syms x t 𝐞
+function convert_axiom(name::AbstractString)
+    @syms 𝐞 a b c d e f g h i j k l m n o p q r s t u v w x y z
 
-r1 = @rule ^(𝐞, ~k) => exp(~k)
-r2 = @rule 𝐞 => MathConstants.e
+    axion_rules = [
+        @rule ^(𝐞, ~k) => exp(~k)
+        @rule 𝐞 => MathConstants.e
+    ]
 
-function convert_axiom(name::AbstractString)
-    f = open(name, "r")
-    r = r"\[.*\]"
+    fd = open(name, "r")
+    re = r"\[.*\]"
     L = []
 
-    for l in readlines(f)
-        l = strip(l)
-        m = match(r, l)
-        if m != nothing
-            l = replace(m.match, "%e" => 𝐞)
-            l = replace(m.match, "%i" => im)
+    D = Dict{Any,Int}()
+
+    for (lineno,line) in enumerate(readlines(fd))
+        line = strip(line)
+        ma = match(re, line)
+        if ma != nothing
+            line = replace(ma.match, "%e" => 𝐞)
+            line  = replace(ma.match, "%i" => im)
+
             try
-                expr = Meta.parse(l)
-                p = eval(expr)
-                p[1] = Prewalk(PassThrough(Chain([r1,r2])))(p[1])
-                p[4] = Prewalk(PassThrough(Chain([r1,r2])))(p[4])
-                push!(L, p)
-            catch e
+                expr = Meta.parse(line)
+                P = eval(expr)
+
+                print(lineno, ": ", P[1], " => ")
+                P[1] = Prewalk(PassThrough(Chain(axion_rules)))(P[1])
+                P[4] = Prewalk(PassThrough(Chain(axion_rules)))(P[4])
+
+                for ν in get_variables(P[1])
+                    if !isequal(ν, x) && !haskey(D, ν)
+                        D[ν] = length(D) + 1
+                    end
+                end
+
+                P[1] = substitute(P[1], D)
+                P[4] = substitute(P[4], D)
+
+                println(P[1])
+                push!(L, P)
+            catch err
             end
         end
     end
 
-    close(f)
+    close(fd)
     return L
 end