Skip to content

Commit 6133194

Browse files
shahriariravanianshahriariravanian
authored
Merge pull request #45 from shahriariravanian/main
Improvements in the integration rules and separating cache
2 parents 99004ee + 63dfad9 commit 6133194

File tree

10 files changed

+198
-74
lines changed

10 files changed

+198
-74
lines changed

‎Project.toml‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "SymbolicNumericIntegration"
22
uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
33
authors = ["Shahriar Iravanian <siravan@svtsim.com>"]
4-
version = "1.0.0"
4+
version = "1.0.1"
55

66
[deps]
77
DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"
@@ -27,3 +27,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2727

2828
[targets]
2929
test = ["PyCall", "SymPy", "Test"]
30+

‎src/SymbolicNumericIntegration.jl‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,17 +2,16 @@ module SymbolicNumericIntegration
22

33
using SymbolicUtils
44
using SymbolicUtils: istree, operation, arguments
5+
using Symbolics
56
using Symbolics: value, get_variables, expand_derivatives
67
using SymbolicUtils.Rewriters
7-
using Symbolics
88

99
using DataDrivenDiffEq
1010

1111
include("utils.jl")
12+
include("cache.jl")
1213
include("roots.jl")
1314

14-
export find_roots
15-
1615
include("rules.jl")
1716
include("candidates.jl")
1817
include("homotopy.jl")

‎src/cache.jl‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
const DEBUG_CACHE = true
2+
3+
mutable struct ExprCache
4+
eq::Any# the primary expression
5+
f::Any# compiled eq
6+
δeq::Any# the symbolic derivative of eq
7+
δf::Any# compiled δeq
8+
end
9+
10+
cache(eq::ExprCache) = eq
11+
cache(eq) = ExprCache(eq, nothing, nothing, nothing)
12+
13+
expr(c::ExprCache) = c.eq
14+
expr(c) = c
15+
16+
function deriv!(c::ExprCache, x)
17+
if c.δeq == nothing
18+
c.δeq = expand_derivatives(Differential(x)(expr(c)))
19+
end
20+
return c.δeq
21+
end
22+
23+
function deriv!(c, x)
24+
if DEBUG_CACHE
25+
error("ExprCache object expected")
26+
end
27+
return expand_derivatives(Differential(x)(c))
28+
end
29+
30+
function fun!(c::ExprCache, x)
31+
if c.f == nothing
32+
c.f = build_function(expr(c), x; expression = false)
33+
end
34+
return c.f
35+
end
36+
37+
function fun!(c, x)
38+
if DEBUG_CACHE
39+
error("ExprCache object expected")
40+
end
41+
return build_function(c, x; expression = false)
42+
end
43+
44+
function deriv_fun!(c::ExprCache, x)
45+
if c.δf == nothing
46+
c.δf = build_function(deriv!(c, x), x; expression = false)
47+
end
48+
return c.δf
49+
end
50+
51+
function deriv_fun!(c, x)
52+
if DEBUG_CACHE
53+
error("ExprCache object expected")
54+
end
55+
return build_function(deriv!(c, x), x; expression = false)
56+
end
57+
58+
Base.show(io::IO, c::ExprCache) = show(expr(c))

‎src/candidates.jl‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ using DataStructures
33
# this is the main heurisctic used to find the test fragments
44
function generate_basis(eq, x; homotopy = true)
55
# if homotopy return generate_homotopy(eq, x) end
6-
eq = expand(eq)
6+
eq = expand(expr(eq))
77
S = 0 #Set{Any}()
88
for t in terms(eq)
99
q = equivalent(t, x)
@@ -24,16 +24,17 @@ function generate_basis(eq, x; homotopy = true)
2424

2525
S += sum(c₁ * c₂ for c₁ in C₁ for c₂ in C₂)
2626
end
27-
return unique([one(x); [equivalent(t, x) for t in terms(S)]])
27+
return cache.(unique([one(x); [equivalent(t, x) for t in terms(S)]]))
2828
end
2929

3030
function expand_basis(basis, x)
31-
b = sum(basis)
32-
eq = (1 + x) * (b + Differential(x)(b))
33-
eq = expand(expand_derivatives(eq))
31+
b = sum(expr.(basis))
32+
δb = sum(deriv!.(basis, x))
33+
eq = (1 + x) * (b + δb)
34+
eq = expand(eq)
3435
S = Set{Any}()
3536
enqueue_expr!(S, eq, x)
36-
return [one(x); [s for s in S]]
37+
return cache.([one(x); [s for s in S]])
3738
end
3839

3940
function closure(eq, x; max_terms = 50)

‎src/homotopy.jl‎

Lines changed: 37 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -14,13 +14,13 @@ function next_variable!(f, eq)
1414
end
1515

1616
function transformer(eq::SymbolicUtils.Add, f)
17-
sum(transformer(t, f) for t in arguments(eq); init = 0)
17+
return sum(transformer(t, f) for t in arguments(eq); init = 0)
1818
end
1919
function transformer(eq::SymbolicUtils.Mul, f)
20-
prod(transformer(t, f) for t in arguments(eq); init = 1)
20+
return prod(transformer(t, f) for t in arguments(eq); init = 1)
2121
end
2222
function transformer(eq::SymbolicUtils.Div, f)
23-
transformer(arguments(eq)[1], f) * transformer(inv(arguments(eq)[2]), f)
23+
return transformer(arguments(eq)[1], f) * transformer(arguments(eq)[2]^-1, f)
2424
end
2525

2626
function transformer(eq::SymbolicUtils.Pow, f)
@@ -63,6 +63,7 @@ function substitute_x(eq, x, sub)
6363
end
6464

6565
function generate_homotopy(eq, x)
66+
eq = eq isa Num ? eq.val : eq
6667
q, sub = transform(eq, x)
6768
S = 0
6869

@@ -87,30 +88,44 @@ function ∂(x)
8788
return isequal(d, 0) ? 1 : d
8889
end
8990

90-
partial_int_rules = [@rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
91-
@rule 𝛷(^(cos(~x), ~k::is_neg)) => 𝛷(^(sec(~x), -~k))
92-
@rule 𝛷(^(tan(~x), ~k::is_neg)) => 𝛷(^(cot(~x), -~k))
93-
@rule 𝛷(^(csc(~x), ~k::is_neg)) => 𝛷(^(sin(~x), -~k))
94-
@rule 𝛷(^(sec(~x), ~k::is_neg)) => 𝛷(^(cos(~x), -~k))
95-
@rule 𝛷(^(cot(~x), ~k::is_neg)) => 𝛷(^(tan(~x), -~k))
96-
@rule 𝛷(^(sinh(~x), ~k::is_neg)) => 𝛷(^(csch(~x), -~k))
97-
@rule 𝛷(^(cosh(~x), ~k::is_neg)) => 𝛷(^(sech(~x), -~k))
98-
@rule 𝛷(^(tanh(~x), ~k::is_neg)) => 𝛷(^(coth(~x), -~k))
99-
@rule 𝛷(^(csch(~x), ~k::is_neg)) => 𝛷(^(sinh(~x), -~k))
100-
@rule 𝛷(^(sech(~x), ~k::is_neg)) => 𝛷(^(cosh(~x), -~k))
101-
@rule 𝛷(^(coth(~x), ~k::is_neg)) => 𝛷(^(tanh(~x), -~k))
102-
@rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
91+
partial_int_rules = [@rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
10392
@rule 𝛷(cos(~x)) => (sin(~x), ∂(~x))
10493
@rule 𝛷(tan(~x)) => (log(cos(~x)), ∂(~x))
105-
@rule 𝛷(csc(~x)) => (log(sin(~x)^-1 - cos(~x) * sin(~x)^-1), ∂(~x))
106-
@rule 𝛷(sec(~x)) => (log(cos(~x)^-1 + sin(~x) * cos(~x)^-1), ∂(~x))
94+
@rule 𝛷(csc(~x)) => (log(csc(~x) + cot(~x)), ∂(~x))
95+
@rule 𝛷(sec(~x)) => (log(sec(~x) + tan(~x)), ∂(~x))
10796
@rule 𝛷(cot(~x)) => (log(sin(~x)), ∂(~x))
10897
@rule 𝛷(sinh(~x)) => (cosh(~x), ∂(~x))
10998
@rule 𝛷(cosh(~x)) => (sinh(~x), ∂(~x))
11099
@rule 𝛷(tanh(~x)) => (log(cosh(~x)), ∂(~x))
111-
@rule 𝛷(csch(~x)) => (log(sinh(~x)^-1 + cosh(~x) * sinh(~x)^-1), ∂(~x))
112-
@rule 𝛷(sech(~x)) => (log(cosh(~x)^-1 + sinh(~x) * cosh(~x)^-1), ∂(~x))
100+
@rule 𝛷(csch(~x)) => (log(tanh(~x / 2)), ∂(~x))
101+
@rule 𝛷(sech(~x)) => (atan(sinh(~x)), ∂(~x))
113102
@rule 𝛷(coth(~x)) => (log(sinh(~x)), ∂(~x))
103+
@rule 𝛷(^(sin(~x), -1)) => (log(csc(~x) + cot(~x)), ∂(~x))
104+
@rule 𝛷(^(cos(~x), -1)) => (log(sec(~x) + tan(~x)), ∂(~x))
105+
@rule 𝛷(^(tan(~x), -1)) => (log(sin(~x)), ∂(~x))
106+
@rule 𝛷(^(csc(~x), -1)) => (cos(~x), ∂(~x))
107+
@rule 𝛷(^(sec(~x), -1)) => (sin(~x), ∂(~x))
108+
@rule 𝛷(^(cot(~x), -1)) => (log(cos(~x)), ∂(~x))
109+
@rule 𝛷(^(sinh(~x), -1)) => (log(tanh(~x / 2)), ∂(~x))
110+
@rule 𝛷(^(cosh(~x), -1)) => (atan(sinh(~x)), ∂(~x))
111+
@rule 𝛷(^(tanh(~x), -1)) => (log(sinh(~x)), ∂(~x))
112+
@rule 𝛷(^(csch(~x), -1)) => (cosh(~x), ∂(~x))
113+
@rule 𝛷(^(sech(~x), -1)) => (sinh(~x), ∂(~x))
114+
@rule 𝛷(^(coth(~x), -1)) => (log(cosh(~x)), ∂(~x))
115+
116+
# @rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
117+
# @rule 𝛷(^(cos(~x), ~k::is_neg)) => 𝛷(^(sec(~x), -~k))
118+
# @rule 𝛷(^(tan(~x), ~k::is_neg)) => 𝛷(^(cot(~x), -~k))
119+
# @rule 𝛷(^(csc(~x), ~k::is_neg)) => 𝛷(^(sin(~x), -~k))
120+
# @rule 𝛷(^(sec(~x), ~k::is_neg)) => 𝛷(^(cos(~x), -~k))
121+
# @rule 𝛷(^(cot(~x), ~k::is_neg)) => 𝛷(^(tan(~x), -~k))
122+
# @rule 𝛷(^(sinh(~x), ~k::is_neg)) => 𝛷(^(csch(~x), -~k))
123+
# @rule 𝛷(^(cosh(~x), ~k::is_neg)) => 𝛷(^(sech(~x), -~k))
124+
# @rule 𝛷(^(tanh(~x), ~k::is_neg)) => 𝛷(^(coth(~x), -~k))
125+
# @rule 𝛷(^(csch(~x), ~k::is_neg)) => 𝛷(^(sinh(~x), -~k))
126+
# @rule 𝛷(^(sech(~x), ~k::is_neg)) => 𝛷(^(cosh(~x), -~k))
127+
# @rule 𝛷(^(coth(~x), ~k::is_neg)) => 𝛷(^(tanh(~x), -~k))
128+
114129
@rule 𝛷(asin(~x)) => (~x * asin(~x) + sqrt(1 - ~x * ~x), ∂(~x))
115130
@rule 𝛷(acos(~x)) => (~x * acos(~x) + sqrt(1 - ~x * ~x), ∂(~x))
116131
@rule 𝛷(atan(~x)) => (~x * atan(~x) + log(~x * ~x + 1), ∂(~x))
@@ -132,6 +147,8 @@ partial_int_rules = [@rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
132147
@rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), 1);
133148
@rule 𝛷(^(sqrt(~x), -1)) => 𝛷(^(~x, -0.5))
134149
@rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
150+
@rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1),
151+
∂(~x))
135152
@rule 𝛷(1 / ~x) => 𝛷(^(~x, -1))
136153
@rule 𝛷(^(~x, ~k)) => (^(~x, ~k + 1), ∂(~x))
137154
@rule 𝛷(exp(~x)) => (exp(~x), ∂(~x))

0 commit comments

Comments
 (0)