Merge pull request #40 from shahriariravanian/main

various performance improvements based of profiling

Author: shahriariravanian

Project.toml

@@ -1,11 +1,12 @@
 name = "SymbolicNumericIntegration"
 uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 authors = ["Shahriar Iravanian <[email protected]>"]
-version = "0.8.6"
+version = "0.9.0"
 
 [deps]
 DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
@@ -15,9 +16,9 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 DataDrivenDiffEq = "0.8"
 DataStructures = "0.18"
 PyCall = "1"
+SymPy = "1"
 SymbolicUtils = "0.19"
 Symbolics = "4"
-SymPy = "1"
 julia = "1.6"
 
 [extras]

src/integral.jl

@@ -280,16 +280,17 @@ function init_basis_matrix!(T, A, X, x, eq, Δbasis, radius, complex_plane; abst
     k = 1
     i = 1
 
+    eq_fun = build_function(eq, x; expression = false)
+    Δbasis_fun = build_function.(Δbasis, x; expression = false)
+
     while k <= n
         try
-            x₀ = test_point(complex_plane, radius)
-            X[k] = x₀
-            d = Dict(x => x₀)
+            X[k] = test_point(complex_plane, radius)
+            b₀ = eq_fun(X[k])
 
-            b₀ = Complex{T}(substitute(eq, d))
             if is_proper(b₀)
                 for j in 1:n
-                    A[k, j] = Complex{T}(substitute(Δbasis[j], d)) / b₀
+                    A[k, j] = Δbasis_fun[j](X[k]) / b₀
                 end
                 if all(is_proper, A[k, :])
                     k += 1
@@ -304,17 +305,21 @@ end
 function modify_basis_matrix!(T, A, X, x, eq, ∂eq, Δbasis, radius; abstol = 1e-6)
     n = size(A, 1)
     k = 1
+
+    eq_fun = build_function(eq, x; expression = false)
+    ∂eq_fun = build_function(∂eq, x; expression = false)
+    Δbasis_fun = build_function.(Δbasis, x; expression = false)
+
     for k in 1:n
-        d = Dict(x => X[k])
         # One Newton iteration toward the poles
         # note the + sign instead of the usual - in Newton-Raphson's method. This is
         # because we are moving toward the poles and not zeros.
-        x₀ = X[k] + Complex{T}(substitute(eq, d)) / Complex{T}(substitute(∂eq, d))
+
+        x₀ = X[k] + eq_fun(X[k]) / ∂eq_fun(X[k])
         X[k] = x₀
-        d = Dict(x => x₀)
-        b₀ = Complex{T}(substitute(eq, d))
+        b₀ = eq_fun(X[k])
         for j in 1:n
-            A[k, j] = Complex{T}(substitute(Δbasis[j], d)) / b₀
+            A[k, j] = Δbasis_fun[j](X[k]) / b₀
         end
     end
 end

src/numeric_utils.jl

@@ -31,42 +31,28 @@ end
 """
     converts float to int or small rational numbers
 """
-function nice_parameters(p; abstol = 1e-3)
-    c = lcm(collect(1:10)...)
-    n = length(p)
-    q = Array{Any}(undef, n)
-    for i in 1:n
-        den = 1
-        while den < 10
-            if abs(round(p[i] * den) - p[i] * den) < abstol
-                a = round(Int, p[i] * den) // den
-                q[i] = (denominator(a) == 1 ? numerator(a) : a)
-                den = 10
-            else
-                q[i] = Float64(p[i])
-            end
-            den += 1
-        end
-    end
-    q
+
+function nice_parameter(u::Complex{T}; abstol = 1e-6) where {T <: Real}
+    α = nice_parameter(real(u); abstol)
+    β = nice_parameter(imag(u); abstol)
+    return β ≈ 0 ? α : Complex(α, β)
 end
 
-function nice_parameter(u::T; abstol = 1e-3, M = 10) where {T <: Real}
-    c = lcm(collect(1:M)...)
-    for den in 1:M
+nice_parameter(x; abstol = 1e-6) = small_rational(x; abstol)
+
+nice_parameters(p; abstol = 1e-6) = nice_parameter.(p; abstol)
+
+function small_rational(x::T; abstol = 1e-6, M = 20) where {T <: Real}
+    a = floor(Int, x)
+    r = x - a
+    for den in 2:M
         try
-            if abs(round(u * den) - u * den) < abstol
-                a = round(Int, u * den) // den
-                return (denominator(a) == 1 ? numerator(a) : a)
+            if abs(round(r * den) - r * den) < abstol
+                s = a + round(Int, r * den) // den
+                return (denominator(s) == 1 ? numerator(s) : s)
             end
         catch e
         end
     end
-    return u
-end
-
-function nice_parameter(u::Complex{T}; abstol = 1e-3, M = 10) where {T <: Real}
-    α = nice_parameter(real(u))
-    β = nice_parameter(imag(u))
-    return β ≈ 0 ? α : Complex(α, β)
+    return x
 end

src/utils.jl

@@ -1,7 +1,12 @@
 """
     isdependent returns true if eq is dependent on x
 """
-isdependent(eq, x) = !isequal(expand_derivatives(Differential(x)(eq)), 0)
+#isdependent(eq, x) = !isequal(expand_derivatives(Differential(x)(eq)), 0)
+
+function isdependent(eq, x)
+    vars = get_variables(eq)
+    return length(vars) == 1 && isequal(x, vars[1])
+end
 
 """
     is_number(x) returns true if x is a concrete numerical type

test/axiom.jl

@@ -1,4 +1,3 @@
-
 @variables x a b c d e p t m n z
 
 function convert_axiom(name::AbstractString)
@@ -104,7 +103,7 @@ function accept_integral(sol, ans, x; radius = 1.0, abstol = 1e-3, n = 5)
     return false
 end
 
-function test_axiom(L, try_sympy = true; kwargs...)
+function test_axiom(L; kwargs...)
     n_ok = 0
     n_fail = 0
     n_diff = 0
@@ -120,7 +119,7 @@ function test_axiom(L, try_sympy = true; kwargs...)
         printstyled(eq, '\n'; color = :green)
         sol = integrate(eq, x; kwargs...)[1]
 
-        println(">>>>", sol)
+        println(">>>>\t", sol)
 
         if isequal(sol, 0)
             printstyled("\tFailure\n"; color = :red)
@@ -133,13 +132,6 @@ function test_axiom(L, try_sympy = true; kwargs...)
             n_diff += 1
         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)
         # catch e
         #    printstyled(i, ": ", e, '\n'; color=:red)

test/runtests.jl

@@ -1,11 +1,11 @@
 using SymbolicNumericIntegration
+using SymbolicNumericIntegration: value
 using Symbolics
 
 using SymbolicUtils
 using SymbolicUtils.Rewriters
 
 using Test
-using PyCall, SymPy
 
 include("axiom.jl")