try SymPy.jl

Author: anandijain

Project.toml

@@ -12,16 +12,20 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
+Conda = "1"
 DataDrivenDiffEq = "0.8"
 DataStructures = "0.18"
 PyCall = "1"
 SymbolicUtils = "0.19"
 Symbolics = "4"
+SymPy = "1"
 julia = "1.6"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "PyCall"]
+test = ["Conda", "PyCall", "SymPy", "Test"]

test/axiom.jl

@@ -1,10 +1,4 @@
-using SymbolicUtils
-using SymbolicUtils.Rewriters
-using Symbolics
-# using SymbolicNumericIntegration
 
-using PyCall
-sympy = pyimport("sympy")
 
 @variables x a b c d e p t m n z
 
@@ -15,16 +9,18 @@ function convert_axiom(name::AbstractString)
 
     D = Dict{Any,Int}()
 
-    for (lineno,line) in enumerate(readlines(fd))
+    for (lineno, line) in enumerate(readlines(fd))
         line = strip(line)
-        if length(line) < 3 || line[1:2] == "--" continue end
+        if length(line) < 3 || line[1:2] == "--"
+            continue
+        end
         ma = match(re, line)
 
         if ma != nothing
             line = ma.match
             println(line)
             line = replace(line, "%e" => MathConstants.e)
-            line  = replace(line, "%i" => im)
+            line = replace(line, "%i" => im)
 
             try
                 expr = Meta.parse(line)
@@ -41,7 +37,7 @@ function convert_axiom(name::AbstractString)
                 P[1] = substitute(P[1], D)
                 P[4] = substitute(P[4], D)
 
-                printstyled(length(L)+1, ": "; color=:red)
+                printstyled(length(L) + 1, ": "; color=:red)
                 println(P[1])
                 push!(L, P)
             catch err
@@ -83,16 +79,16 @@ function convert_axiom(sym)
             "Wester Problems.input"
         end
 
-        name = joinpath("test/0", name)
-        println(name)
-        return convert_axiom(name)
+    name = joinpath("test/0", name)
+    println(name)
+    return convert_axiom(name)
 end
 
 function test_point(complex_plane, radius)
     if complex_plane
-        return radius * sqrt(rand()) * cis(2π*rand())
+        return radius * sqrt(rand()) * cis(2π * rand())
     else
-        return Complex(radius * (2*rand() - 1))
+        return Complex(radius * (2 * rand() - 1))
     end
 end
 
@@ -116,36 +112,37 @@ function test_axiom(L, try_sympy=true; kwargs...)
     n_diff = 0
     n_catch = 0
 
-    for (i,p) in enumerate(L)
+    for (i, p) in enumerate(L)
         # try
-            printstyled(i, ": "; color=:yellow)
-            eq = p[1]
-            x = p[2]
-            ans = p[4]
-
-            printstyled(eq, '\n'; color=:green)
-            sol = integrate(eq, x; kwargs...)[1]
-
-            println(">>>>", sol)
-
-            if isequal(sol, 0)
-                printstyled("\tFailure\n"; color=:red)
-                n_fail += 1
-            elseif accept_integral(value(sol), value(ans), x)
-                printstyled("\tSuccess:\t", sol, '\n'; color=:cyan)
-                n_ok += 1
-            else
-                printstyled("\tDiscrepancy:\t", sol, '\n'; color=:yellow)
-                n_diff += 1
-            end
+        printstyled(i, ": "; color=:yellow)
+        eq = p[1]
+        x = p[2]
+        ans = p[4]
+
+        printstyled(eq, '\n'; color=:green)
+        sol = integrate(eq, x; kwargs...)[1]
+
+        println(">>>>", sol)
+
+        if isequal(sol, 0)
+            printstyled("\tFailure\n"; color=:red)
+            n_fail += 1
+        elseif accept_integral(value(sol), value(ans), x)
+            printstyled("\tSuccess:\t", sol, '\n'; color=:cyan)
+            n_ok += 1
+        else
+            printstyled("\tDiscrepancy:\t", sol, '\n'; color=:yellow)
+            n_diff += 1
+        end
 
-            if try_sympy
-                s = pythonize(eq)
-                py = sympy.integrate(s, sympy.Symbol(string(x)))
-                printstyled("\tSymPy      :\t", string(py)[10:end], '\n'; color=:magenta)
-            end
+        if try_sympy
+            s = Symbolics.symbolics_to_sympy(s)
+            s_x = Symbolics.symbolics_to_sympy(x)
+            py = SymPy.integrate(s, s_x)
+            printstyled("\tSymPy      :\t", string(py)[10:end], '\n'; color=:magenta)
+        end
 
-            printstyled("\tAnswer: \t", ans, '\n'; color=:blue)
+        printstyled("\tAnswer: \t", ans, '\n'; color=:blue)
         # catch e
         #    printstyled(i, ": ", e, '\n'; color=:red)
         #    n_catch += 1
@@ -158,22 +155,3 @@ function test_axiom(L, try_sympy=true; kwargs...)
     printstyled("Discrepancy = ", n_diff, '\n'; color=:yellow)
     printstyled("Exception   = ", n_catch, '\n'; color=:cyan)
 end
-
-#############################################################################
-
-pythonize(eq::SymbolicUtils.Add) = "(" * join(pythonize.(arguments(eq)), ")+(") * ")"
-pythonize(eq::SymbolicUtils.Mul) = "(" * join(pythonize.(arguments(eq)), ")*(") * ")"
-
-function pythonize(eq::SymbolicUtils.Pow)
-    a = pythonize.(arguments(eq))
-    "(" * a[1] * ")**(" * a[2] * ")"
-end
-
-pythonize(eq::SymbolicUtils.Term) = string(operation(eq)) * '(' * pythonize(arguments(eq)[1]) * ')'
-
-function pythonize(eq)
-    s = string(eq)
-    s = replace(s, "π" => "pi")
-    s = replace(s, "im" => "j")
-    s
-end

test/runtests.jl

@@ -1,6 +1,11 @@
 using SymbolicNumericIntegration
 using Symbolics
+
+using SymbolicUtils
+using SymbolicUtils.Rewriters
+
 using Test
+using PyCall, Conda, SymPy
 
 include("axiom.jl")
 
@@ -13,26 +18,26 @@ include("axiom.jl")
     based on http://integral-table.com/ with modifications
 """
 basic_integrals = [
-# Basic Forms
+    # Basic Forms
     1,
     x^2,
     4x^3,
-# Integrals of Rational Functions
-    1/x,
-    1/(2x + 5),
-    1/(x + 1)^2,
+    # Integrals of Rational Functions
+    1 / x,
+    1 / (2x + 5),
+    1 / (x + 1)^2,
     (x + 3)^3,
-    x*(x - 2)^4,
-    1/(1 + x^2),
-    1/(9 + x^2),
-    x/(4 + x^2),
-    x^2/(16 + x^2),
-    x^3/(1 + x^2),
-    1/(x^2 - 5x + 6),
-    1/(x^2 + x + 1),
-    x/(x + 4)^2,
-    x/(x^2 + x + 1),
-# Integrals with Roots
+    x * (x - 2)^4,
+    1 / (1 + x^2),
+    1 / (9 + x^2),
+    x / (4 + x^2),
+    x^2 / (16 + x^2),
+    x^3 / (1 + x^2),
+    1 / (x^2 - 5x + 6),
+    1 / (x^2 + x + 1),
+    x / (x + 4)^2,
+    x / (x^2 + x + 1),
+    # Integrals with Roots
     sqrt(x - 2),
     1 / sqrt(x - 1),
     1 / sqrt(x + 1),
@@ -45,8 +50,8 @@ basic_integrals = [
     sqrt(x / (4 - x)),
     sqrt(x / (4 + x)),
     x * sqrt(2x + 3),
-    sqrt(x*(x+2)),
-    sqrt(x^3*(x+3)),
+    sqrt(x * (x + 2)),
+    sqrt(x^3 * (x + 3)),
     sqrt(x^2 + 4),
     sqrt(x^2 - 4),
     sqrt(4 - x^2),
@@ -64,7 +69,7 @@ basic_integrals = [
     x * sqrt(x^2 - 5x + 6),
     1 / sqrt(x^2 - 5x + 6),
     1 / (4 + x^2)^1.5,
-# Integrals with Logarithms
+    # Integrals with Logarithms
     log(x),
     x * log(x),
     x^2 * log(x),
@@ -80,7 +85,7 @@ basic_integrals = [
     log(x)^3,
     x * log(x)^2,
     x^2 * log(x)^2,
-# Integrals with Exponentials
+    # Integrals with Exponentials
     exp(x),
     sqrt(x) * exp(x),
     x * exp(x),
@@ -91,7 +96,7 @@ basic_integrals = [
     x^3 * exp(2x),
     exp(x^2),
     x * exp(x^2),
-# Integrals with Trigonometric Functions
+    # Integrals with Trigonometric Functions
     sin(4x),
     sin(x)^2,
     sin(x)^3,
@@ -116,7 +121,7 @@ basic_integrals = [
     sec(x)^2 * tan(x),
     csc(x),
     sec(x) * csc(x),
-# Products of Trigonometric Functions and Monomials
+    # Products of Trigonometric Functions and Monomials
     x * cos(x),
     x * cos(3x),
     x^2 * cos(x),
@@ -132,14 +137,14 @@ basic_integrals = [
     x^3 * sin(x),
     x^4 * cos(2x),
     sin(x)^2 * cos(x)^3,
-# Products of Trigonometric Functions and Exponentials
+    # Products of Trigonometric Functions and Exponentials
     exp(x) * sin(x),
     exp(3x) * sin(2x),
     exp(x) * cos(x),
     exp(2x) * cos(7x),
     x * exp(x) * sin(x),
     x * exp(x) * cos(x),
-# Integrals of Hyperbolic Functions
+    # Integrals of Hyperbolic Functions
     cosh(x),
     exp(x) * cosh(x),
     sinh(3x),
@@ -152,15 +157,15 @@ basic_integrals = [
     sin(x) * sinh(x),
     sinh(x) * cosh(x),
     sinh(3x) * cosh(5x),
-# Misc
+    # Misc
     exp(x) / (1 + exp(x)),
     cos(exp(x)) * sin(exp(x)) * exp(x),
     cos(exp(x))^2 * sin(exp(x)) * exp(x),
-    1 / (x*log(x)),
+    1 / (x * log(x)),
     (log(x - 1) + (x - 1)^-1) * log(x),
     1 / (exp(x) - 1),
     1 / (exp(x) + 5),
-    sqrt(x)*log(x),
+    sqrt(x) * log(x),
     log(log(x)) / x,
     x^3 * exp(x^2),
     sin(log(x)),
@@ -169,24 +174,24 @@ basic_integrals = [
     1 / (exp(2x) - 1),
     exp(x) / (exp(2x) - 1),
     x / (exp(2x) - 1),
-# derivative-divide examples (Lamangna 7.10.2)
+    # derivative-divide examples (Lamangna 7.10.2)
     exp(x) * exp(exp(x)),
     exp(sqrt(x)) / sqrt(x),
-    log(log(x)) / (x*log(x)),
+    log(log(x)) / (x * log(x)),
     log(cos(x)) * tan(x),
-# rothstein-Trager examples (Lamangna 7.10.9)
+    # rothstein-Trager examples (Lamangna 7.10.9)
     1 / (x^3 - x),
     1 / (x^3 + 1),
     1 / (x^2 - 8),
     (x + 1) / (x^2 + 1),
     x / (x^4 - 4),
     x^3 / (x^4 + 1),
     1 / (x^4 + 1),
-# bypass = true
+    # bypass = true
     β,      # turn of bypass = true
-    exp(x)/x - exp(x)/x^2,
-    cos(x)/x - sin(x)/x^2,
-    1/log(x) - 1/log(x)^2,
+    exp(x) / x - exp(x) / x^2,
+    cos(x) / x - sin(x) / x^2,
+    1 / log(x) - 1 / log(x)^2,
 ]
 
 function test_integrals(; kw...)
@@ -203,7 +208,7 @@ function test_integrals(; kw...)
             printstyled(k, ": "; color=:blue)
             k += 1
             printstyled(eq, " =>\n"; color=:green)
-            solved, unsolved = integrate(eq; args...)
+            solved, unsolved = SymbolicNumericIntegration.integrate(eq; args...)
             printstyled('\t', solved; color=:cyan)
             if isequal(unsolved, 0)
                 println()
@@ -215,14 +220,16 @@ function test_integrals(; kw...)
     end
 
     n = length(misses)
-    if n > 0 println("**** missess (n=$n) *****") end
+    if n > 0
+        println("**** missess (n=$n) *****")
+    end
     for eq in misses
         printstyled(eq, '\n'; color=:red)
     end
     return n
 end
 
-@testset "integral" begin 
-    n = test_integrals(;symbolic=false, verbose=false, homotopy=true, num_steps=2, num_trials=10)
-    @test n==11
+@testset "integral" begin
+    n = test_integrals(; symbolic=false, verbose=false, homotopy=true, num_steps=2, num_trials=10)
+    @test n == 10 
 end