performance improvements (especially in init_basis_matrix)

Author: shahriariravanian

src/integral.jl

@@ -279,17 +279,18 @@ function init_basis_matrix!(T, A, X, x, eq, Δbasis, radius, complex_plane; abst
     n = size(A, 1)
     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₀)
-
-            b₀ = Complex{T}(substitute(eq, d))
+        try            
+            X[k] = test_point(complex_plane, radius)
+            b₀ = eq_fun(X[k])
+            
             if is_proper(b₀)
-                for j in 1:n
-                    A[k, j] = Complex{T}(substitute(Δbasis[j], d)) / b₀
+                for j in 1:n            
+                    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,33 @@ 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-3, M = 10) where {T <: Real}
+    α = nice_parameter(real(u))
+    β = nice_parameter(imag(u))
+    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-7) = small_rational(x; abstol)
+
+nice_parameters(p; abstol = 1e-3) = nice_parameter.(p)
+
+function small_rational(x::T; abstol=1e-6, M=20) where {T <: Real}
+    # c = lcm(collect(1:M)...)
+   	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
+    return x
 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(α, β)
-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