Merge pull request #50 from shahriariravanian/main

upgraded to work with Symbolics 5

Author: shahriariravanian

Project.toml

@@ -1,10 +1,11 @@
 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"
+DataDrivenSparse = "5b588203-7d8b-4fab-a537-c31a7f73f46b"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
@@ -14,12 +15,12 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-DataDrivenDiffEq = "0.8"
+DataDrivenDiffEq = "1.0"
 DataStructures = "0.18"
 Optim = "1"
+SymbolicUtils = "1"
+Symbolics = "5"
 SpecialFunctions = "2"
-SymbolicUtils = "0.19"
-Symbolics = "4"
 julia = "1.6"
 
 [extras]

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/SymbolicNumericIntegration.jl

@@ -5,8 +5,9 @@ using SymbolicUtils: istree, operation, arguments
 using Symbolics
 using Symbolics: value, get_variables, expand_derivatives
 using SymbolicUtils.Rewriters
+using SymbolicUtils: exprtype, BasicSymbolic
 
-using DataDrivenDiffEq
+using DataDrivenDiffEq, DataDrivenSparse
 
 include("utils.jl")
 include("special.jl")

src/candidates.jl

@@ -1,10 +1,10 @@
 using DataStructures
 
 # this is the main heurisctic used to find the test fragments
-function generate_basis(eq, x, try_kernel = true; homotopy = true)
-    if homotopy && !try_kernel
+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
@@ -13,33 +13,39 @@ function generate_basis(eq, x, try_kernel = true; homotopy = true)
     for t in terms(eq)
         q = equivalent(t, x)
         f = kernel(q)
+        p = q / f
 
-        C₁ = closure(f, x)
-        p = q * inverse(f)
-
-        if homotopy
-            if isdependent(p, x)
-                C₂ = generate_homotopy(p, x)
-            else
-                C₂ = 1
-            end
+        if isdependent(p, x)
+            C₂ = generate_homotopy(p, x)
         else
-            C₂ = find_candidates(p, x)
+            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)
+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
+
     δ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]])
+    return cache.([s for s in S]), true
 end
 
 function closure(eq, x; max_terms = 50)
@@ -55,23 +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]])
-end
-
-function find_candidates(eq, x)
-    eq = apply_d_rules(eq)
-    D = Differential(x)
-
-    S = Set{Any}()
-    q = Queue{Any}()
-    enqueue_expr!(S, q, eq, x)
-
-    for y in q
-        ∂y = expand_derivatives(D(y))
-        enqueue_expr!(S, ∂y, x)
-    end
-
-    return unique([one(x); [s 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)
@@ -142,27 +132,31 @@ end
 
 ###############################################################################
 
-function enqueue_expr!(S, q, eq::SymbolicUtils.Add, x)
+enqueue_expr!(S, q, eq, x) = enqueue_expr!!(S, q, ops(eq)..., x)
+
+function enqueue_expr!!(S, q, ::Add, eq, x)
     for t in arguments(eq)
         enqueue_expr!(S, q, t, x)
     end
 end
 
-function enqueue_expr!(S, q, eq, x)
+function enqueue_expr!!(S, q, ::Any, eq, x)
     y = eq / coef(eq, x)
     if y ∉ S && isdependent(y, x)
         enqueue!(q, y)
         push!(S, y)
     end
 end
 
-function enqueue_expr!(S, eq::SymbolicUtils.Add, x)
+enqueue_expr!(S, eq, x) = enqueue_expr!!(S, ops(eq)..., x)
+
+function enqueue_expr!!(S, ::Add, eq, x)
     for t in arguments(eq)
         enqueue_expr!(S, t, x)
     end
 end
 
-function enqueue_expr!(S, eq, x)
+function enqueue_expr!!(S, ::Any, eq, x)
     y = eq / coef(eq, x)
     if y ∉ S && isdependent(y, x)
         push!(S, y)