upgraded to use Symbolics 5

Author: shahriariravanian

Project.toml

@@ -1,7 +1,7 @@
 name = "SymbolicNumericIntegration"
 uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 authors = ["Shahriar Iravanian <[email protected]>"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"

docs/Project.toml

@@ -4,4 +4,4 @@ SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 
 [compat]
 Documenter = "0.27"
-SymbolicNumericIntegration = "0.9, 1"
+SymbolicNumericIntegration = "1"

docs/make.jl

@@ -1,7 +1,7 @@
 using Documenter, SymbolicNumericIntegration
 
-cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
-cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
+# cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
+# cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
 
 include("pages.jl")
 
@@ -11,7 +11,7 @@ makedocs(sitename = "SymbolicNumericIntegration.jl",
          clean = true, doctest = false, linkcheck = true,
          strict = [
              :doctest,
-             :linkcheck,
+             # :linkcheck,
              :parse_error,
              :example_block,
              # Other available options are

src/candidates.jl

@@ -4,7 +4,7 @@ using DataStructures
 function generate_basis(eq, x, try_kernel = true)
     if !try_kernel
         S = sum(generate_homotopy(expr(eq), x))
-        return cache.(unique([one(x); [equivalent(t, x) for t in terms(S)]]))
+        return cache.(unique([equivalent(t, x) for t in terms(S)]))
     end
 
     S = 0
@@ -16,32 +16,36 @@ function generate_basis(eq, x, try_kernel = true)
         p = q / f
 
         if isdependent(p, x)
-			C₂ = generate_homotopy(p, x)
+            C₂ = generate_homotopy(p, x)
         else
-        	C₂ = 1
+            C₂ = 1
         end
-        
+
         C₁ = closure(f, x)
-        
+
         S += sum(c₁ * c₂ for c₁ in C₁ for c₂ in C₂)
     end
-    return cache.(unique([one(x); [equivalent(t, x) for t in terms(S)]]))
+    return cache.(unique([equivalent(t, x) for t in terms(S)]))
 end
 
-function expand_basis(basis, x; Kmax=1000)
-    b = sum(expr.(basis))    
+function expand_basis(basis, x; Kmax = 1000)
+    if isempty(basis)
+        return basis, false
+    end
+
+    b = sum(expr.(basis))
+
+    Kb = complexity(b)# Kolmogorov complexity
+    if is_proper(Kb) && Kb > Kmax
+        return basis, false
+    end
 
-    Kb = complexity(b)		# Kolmogorov complexity
-	if is_proper(Kb) && Kb > Kmax
-		return basis, false
-	end
-	
-    δb = sum(deriv!.(basis, x))	    
+    δb = sum(deriv!.(basis, x))
     eq = (1 + x) * (b + δb)
     eq = expand(eq)
     S = Set{Any}()
     enqueue_expr!(S, eq, x)
-    return cache.([one(x); [s for s in S]]), true
+    return cache.([s for s in S]), true
 end
 
 function closure(eq, x; max_terms = 50)
@@ -57,7 +61,7 @@ function closure(eq, x; max_terms = 50)
         y = dequeue!(q)
         enqueue_expr!(S, q, expand_derivatives(D(y)), x)
     end
-    unique([one(x); [s for s in S]; [s * x for s in S]])
+    unique([[s for s in S]; [s * x for s in S]])
 end
 
 function candidate_pow_minus(p, k; abstol = 1e-8)

src/homotopy.jl

@@ -5,56 +5,57 @@ function transformer(::Mul, eq)
 end
 
 function transformer(::Div, eq)
-	a = transformer(arguments(eq)[1])
-	b = transformer(arguments(eq)[2])
-	b = [(1/q, k) for (q, k) in b]
+    a = transformer(arguments(eq)[1])
+    b = transformer(arguments(eq)[2])
+    b = [(1 / q, k) for (q, k) in b]
     return [a; b]
 end
 
 function transformer(::Pow, eq)
     y, k = arguments(eq)
     if is_number(k)
-    	r = nice_parameter(k)
-    	if denominator(r) == 1
-	    	return [(y, k)]
-	    else
-	    	return [(y^(1/denominator(r)), numerator(r))]
-	    end
-    else
-    	return [(eq, 1)]
-	end
+        r = nice_parameter(k)
+        if r isa Integer || r isa Rational
+            if denominator(r) == 1
+                return [(y, k)]
+            else
+                return [(y^(1 / denominator(r)), numerator(r))]
+            end
+        end
+    end
+    return [(eq, 1)]
 end
 
 function transformer(::Any, eq)
-    return [(eq, 1)]    
+    return [(eq, 1)]
 end
 
 function transform(eq, x)
     p = transformer(eq)
-    p = p[isdependent.(first.(p), x)]    
+    p = p[isdependent.(first.(p), x)]
     return p
 end
 
 @syms u[20]
 
-function rename_factors(p, ab)
-	n = length(p)
-	q = 1
-	ks = Int[]
-	sub = Dict()
-
-	for (a,b) in ab
-		sub[a] = b
-	end
-	
-	for (i,(y,k)) in enumerate(p)
-		μ = u[i]
-		q *= μ ^ k
-		sub[μ] = y
-		push!(ks, k)
-	end
-	
-	return q, sub, ks
+function rename_factors(p, ab = ())
+    n = length(p)
+    q = 1
+    ks = Int[]
+    sub = Dict()
+
+    for (a, b) in ab
+        sub[a] = b
+    end
+
+    for (i, (y, k)) in enumerate(p)
+        μ = u[i]
+        q *= μ^k
+        sub[μ] = y
+        push!(ks, k)
+    end
+
+    return q, sub, ks
 end
 
 ##############################################################################
@@ -79,29 +80,26 @@ function generate_homotopy(eq, x)
     eq = eq isa Num ? eq.val : eq
     x = x isa Num ? x.val : x
 
-	p = transform(eq, x)
+    p = transform(eq, x)
     q, sub, ks = rename_factors(p, (si => Si, ci => Ci, ei => Ei, li => Li))
     S = 0
 
     for i in 1:length(ks)
-		μ = u[i]
-		h₁, ∂h₁ = apply_partial_int_rules(sub[μ], x)
-	    h₁ = substitute(h₁, sub)
-		
-    	for j = 1:ks[i]
-		    h₂ = substitute((q / μ^j) / ∂h₁, sub)
-		    S += expand((ω + h₁) * (ω + h₂))
-		end
-    end    
-    
+        μ = u[i]
+        h₁, ∂h₁ = apply_partial_int_rules(sub[μ], x)
+        h₁ = substitute(h₁, sub)
+
+        for j in 1:ks[i]
+            h₂ = substitute((q / μ^j) / ∂h₁, sub)
+            S += expand((ω + h₁) * (ω + h₂))
+        end
+    end
+
     S = substitute(S, Dict(ω => 1))
-    
-    ζ = [x^k for k=1:maximum(ks)+1]
-    
-    unique([one(x); ζ; [equivalent(t, x) for t in terms(S)]])
+    unique([x; [equivalent(t, x) for t in terms(S)]])
 end
 
-##############################################################################
+########################## Main Integration Rules ##################################
 
 @syms 𝛷(x)
 
@@ -159,17 +157,22 @@ partial_int_rules = [
                      @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k);
                                                              init = one(~x)), ~x);
                      @rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), ~x);
-                     @rule 𝛷(1 / sqrt(~x)) => (sum(candidate_sqrt(~x, -0.5); init = one(~x)), ~x);
+                     @rule 𝛷(1 / sqrt(~x)) => (sum(candidate_sqrt(~x, -0.5);
+                                                   init = one(~x)), ~x);
                      # rational functions                                                              
-                     @rule 𝛷(1 / ^(~x::is_poly, ~k::is_pos_int)) => (sum(candidate_pow_minus(~x, -~k);
-                                                                 init = one(~x)), ~x)
+                     @rule 𝛷(1 / ^(~x::is_poly, ~k::is_pos_int)) => (sum(candidate_pow_minus(~x,
+                                                                                             -~k);
+                                                                         init = one(~x)),
+                                                                     ~x)
                      @rule 𝛷(1 / ~x::is_poly) => (sum(candidate_pow_minus(~x, -1);
-                                                                 init = one(~x)), ~x)
+                                                      init = one(~x)), ~x);
                      @rule 𝛷(^(~x, -1)) => (log(~x), ~x)
                      @rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1),
                                                         ~x)
                      @rule 𝛷(1 / ~x) => (log(~x), ~x)
-                     @rule 𝛷(^(~x, ~k::is_pos_int)) => (sum(^(~x, i+1) for i=1:~k+1), ~x)
+                     @rule 𝛷(^(~x, ~k::is_pos_int)) => (sum(^(~x, i + 1)
+                                                            for i in 1:(~k + 1)),
+                                                        ~x)
                      @rule 𝛷(1) => (𝑥, 1)
                      @rule 𝛷(~x) => ((~x + ^(~x, 2)), ~x)]
 
@@ -178,5 +181,3 @@ function apply_partial_int_rules(eq, x)
     D = Differential(x)
     return y, guard_zero(expand_derivatives(D(dy)))
 end
-    
-

src/integral.jl

@@ -135,9 +135,9 @@ function integrate_term(eq, x, l; kwargs...)
     args = Dict(kwargs)
     abstol, num_steps, num_trials, show_basis, symbolic, verbose, max_basis,
     radius = args[:abstol], args[:num_steps],
-                       args[:num_trials], args[:show_basis], args[:symbolic],
-                       args[:verbose],
-                       args[:max_basis], args[:radius]
+             args[:num_trials], args[:show_basis], args[:symbolic],
+             args[:verbose],
+             args[:max_basis], args[:radius]
 
     attempt(l, "Integrating term", eq)
 
@@ -203,11 +203,11 @@ function integrate_term(eq, x, l; kwargs...)
         inform(l, "Failed numeric")
 
         if i < num_steps
-            basis1, ok1 = expand_basis(basis1, x)
-            basis2, ok2 = expand_basis(basis2, x)
-            
-            if !ok1 && ~ok2            	
-            	break
+            basis1, ok1 = expand_basis(prune_basis(eq, x, basis1, radius; kwargs...), x)
+            basis2, ok2 = expand_basis(prune_basis(eq, x, basis2, radius; kwargs...), x)
+
+            if !ok1 && ~ok2
+                break
             end
 
             if show_basis
@@ -230,8 +230,7 @@ end
 ###############################################################################
 
 """
-    the core of the randomized parameter-fitting algorithm
-
+	1try_integrate` is the main dispatch point to call different sparse solvers	
     `try_integrate` tries to find a linear combination of the basis, whose
     derivative is equal to eq
 
@@ -242,7 +241,10 @@ end
 function try_integrate(eq, x, basis, radius; kwargs...)
     args = Dict(kwargs)
     use_optim = args[:use_optim]
-    basis = basis[2:end]    # remove 1 from the beginning
+
+    if isempty(basis)
+        return 0, Inf
+    end
 
     if use_optim
         return solve_optim(eq, x, basis, radius; kwargs...)
@@ -253,12 +255,13 @@ end
 
 #################################################################################
 
-# integrate_basis is used for debugging and should not be called in the course of normal execution
-function integrate_basis(eq, x = var(eq); abstol = 1e-6, radius = 1.0, complex_plane = true)	
+"""
+	integrate_basis is used for debugging and should not be called in the course of normal execution
+"""
+function integrate_basis(eq, x = var(eq); abstol = 1e-6, radius = 1.0, complex_plane = true)
     eq = cache(expand(eq))
     basis = generate_basis(eq, x, false)
     n = length(basis)
-	A, X = init_basis_matrix(eq, x, basis, radius, complex_plane; abstol)
-	return basis, A, X
+    A, X = init_basis_matrix(eq, x, basis, radius, complex_plane; abstol)
+    return basis, A, X
 end
-

src/optim.jl

@@ -1,33 +1,18 @@
 using Optim
 
-function solve_optim(T, eq, x, basis, radius, rounds = 2; kwargs...)
+function solve_optim(eq, x, basis, radius; kwargs...)
     args = Dict(kwargs)
     abstol, complex_plane, verbose = args[:abstol], args[:complex_plane], args[:verbose]
 
     n = length(basis)
     λ = 1e-9
 
-    A = zeros(Complex{T}, (2n, n))
-    X = zeros(Complex{T}, 2n)
-    init_basis_matrix!(T, A, X, x, eq, basis, radius, complex_plane; abstol)
-    # modify_basis_matrix!(T, A, X, x, eq, basis, radius; abstol)
+    A, X = init_basis_matrix(eq, x, basis, radius, complex_plane; abstol)
 
     l = find_independent_subset(A; abstol)
     A, basis, n = A[:, l], basis[l], sum(l)
     p = rank_basis(A, basis)
 
-    # q, ϵ = find_minimizer(A, λ)
-
-    # if ϵ > abstol
-    # 	return 0, ϵ
-    # end
-
-    # qa = q
-    # μ = maximum(abs.(qa))
-
-    # for ρ in exp10.(-1:-1:-8)    
-    # 	l = abs.(qa) .> ρ * μ	    
-
     qm = zeros(n)
     ϵm = 1e6
     lm = qm .> 0