Skip to content

Commit b79569e

Browse files
shahriariravanianshahriariravanian
committed
homotopy rules simplified
1 parent a0b761a commit b79569e

File tree

4 files changed

+87
-99
lines changed

4 files changed

+87
-99
lines changed

‎src/candidates.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ function expand_basis(basis, x; Kmax=1000)
3232
b = sum(expr.(basis))
3333

3434
Kb = complexity(b) # Kolmogorov complexity
35-
if Kb > Kmax
35+
if is_proper(Kb) && Kb > Kmax
3636
return basis, false
3737
end
3838

‎src/homotopy.jl‎

Lines changed: 74 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,3 @@
1-
@syms 𝑥
2-
@syms u[20]
3-
41
transformer(eq) = transformer(ops(eq)...)
52

63
function transformer(::Mul, eq)
@@ -33,19 +30,22 @@ function transformer(::Any, eq)
3330
end
3431

3532
function transform(eq, x)
36-
eq = substitute(eq, Dict(x => 𝑥))
3733
p = transformer(eq)
38-
p = p[isdependent.(first.(p), 𝑥)]
39-
34+
p = p[isdependent.(first.(p), x)]
4035
return p
4136
end
4237

43-
function rename_factors(p)
44-
n = length(p)
38+
@syms u[20]
4539

40+
function rename_factors(p, ab)
41+
n = length(p)
4642
q = 1
47-
sub = Dict()
4843
ks = Int[]
44+
sub = Dict()
45+
46+
for (a,b) in ab
47+
sub[a] = b
48+
end
4949

5050
for (i,(y,k)) in enumerate(p)
5151
μ = u[i]
@@ -69,119 +69,114 @@ Symbolics.derivative(::typeof(Si), args::NTuple{1, Any}, ::Val{1}) = sin(args[1]
6969
Symbolics.derivative(::typeof(Ci), args::NTuple{1, Any}, ::Val{1}) = cos(args[1]) / args[1]
7070
Symbolics.derivative(::typeof(Li), args::NTuple{1, Any}, ::Val{1}) = 1 / log(args[1])
7171

72-
@syms si(𝑥) ci(𝑥) ei(𝑥) li(𝑥)
72+
@syms 𝑥 si(𝑥) ci(𝑥) ei(𝑥) li(𝑥)
7373

7474
##############################################################################
7575

76-
function substitute_x(eq, x, sub)
77-
eq = substitute(eq, sub)
78-
return substitute(eq, Dict(𝑥 => x))
79-
end
80-
8176
guard_zero(x) = isequal(x, 0) ? one(x) : x
8277

8378
function generate_homotopy(eq, x)
8479
eq = eq isa Num ? eq.val : eq
8580
x = x isa Num ? x.val : x
8681

8782
p = transform(eq, x)
88-
q, sub, ks = rename_factors(p)
83+
q, sub, ks = rename_factors(p, (si => Si, ci => Ci, ei => Ei, li => Li))
8984
S = 0
9085

91-
for i in 1:length(sub)
86+
for i in 1:length(ks)
9287
μ = u[i]
93-
h₁, ∂h₁ = apply_partial_int_rules(sub[μ])
94-
h₁ = substitute(h₁, Dict(si => Si, ci => Ci, ei => Ei, li => Li))
95-
h₁ = substitute_x(h₁, x, sub)
88+
h₁, ∂h₁ = apply_partial_int_rules(sub[μ], x)
89+
h₁ = substitute(h₁, sub)
9690

9791
for j = 1:ks[i]
98-
h₂ = substitute_x((q / μ^j) / ∂h₁, x, sub)
99-
S += expand((1 + h₁) * guard_zero(1 + h₂))
92+
h₂ = substitute((q / μ^j) / ∂h₁, sub)
93+
S += expand((ω + h₁) * (ω + h₂))
10094
end
10195
end
10296

97+
S = substitute(S, Dict(ω => 1))
98+
10399
ζ = [x^k for k=1:maximum(ks)+1]
104100

105101
unique([one(x); ζ; [equivalent(t, x) for t in terms(S)]])
106102
end
107103

108104
##############################################################################
109105

110-
function ∂(x)
111-
d = expand_derivatives(Differential(𝑥)(x))
112-
return isequal(d, 0) ? 1 : d
113-
end
114-
115106
@syms 𝛷(x)
116107

117108
partial_int_rules = [
118109
# trigonometric functions
119-
@rule 𝛷(sin(~x)) => (cos(~x) + si(~x), ∂(~x))
120-
@rule 𝛷(cos(~x)) => (sin(~x) + ci(~x), ∂(~x))
121-
@rule 𝛷(tan(~x)) => (log(cos(~x)), ∂(~x))
122-
@rule 𝛷(csc(~x)) => (log(csc(~x) + cot(~x)), ∂(~x))
123-
@rule 𝛷(sec(~x)) => (log(sec(~x) + tan(~x)), ∂(~x))
124-
@rule 𝛷(cot(~x)) => (log(sin(~x)), ∂(~x))
110+
@rule 𝛷(sin(~x)) => (cos(~x) + si(~x), ~x)
111+
@rule 𝛷(cos(~x)) => (sin(~x) + ci(~x), ~x)
112+
@rule 𝛷(tan(~x)) => (log(cos(~x)), ~x)
113+
@rule 𝛷(csc(~x)) => (log(csc(~x) + cot(~x)) + log(sin(~x)), ~x)
114+
@rule 𝛷(sec(~x)) => (log(sec(~x) + tan(~x)) + log(cos(~x)), ~x)
115+
@rule 𝛷(cot(~x)) => (log(sin(~x)), ~x)
125116
# hyperbolic functions
126-
@rule 𝛷(sinh(~x)) => (cosh(~x), ∂(~x))
127-
@rule 𝛷(cosh(~x)) => (sinh(~x), ∂(~x))
128-
@rule 𝛷(tanh(~x)) => (log(cosh(~x)), ∂(~x))
129-
@rule 𝛷(csch(~x)) => (log(tanh(~x / 2)), ∂(~x))
130-
@rule 𝛷(sech(~x)) => (atan(sinh(~x)), ∂(~x))
131-
@rule 𝛷(coth(~x)) => (log(sinh(~x)), ∂(~x))
117+
@rule 𝛷(sinh(~x)) => (cosh(~x), ~x)
118+
@rule 𝛷(cosh(~x)) => (sinh(~x), ~x)
119+
@rule 𝛷(tanh(~x)) => (log(cosh(~x)), ~x)
120+
@rule 𝛷(csch(~x)) => (log(tanh(~x / 2)), ~x)
121+
@rule 𝛷(sech(~x)) => (atan(sinh(~x)), ~x)
122+
@rule 𝛷(coth(~x)) => (log(sinh(~x)), ~x)
132123
# 1/trigonometric functions
133-
@rule 𝛷(1 / sin(~x)) => (log(csc(~x) + cot(~x)) + log(sin(~x)), ∂(~x))
134-
@rule 𝛷(1 / cos(~x)) => (log(sec(~x) + tan(~x)) + log(cos(~x)), ∂(~x))
135-
@rule 𝛷(1 / tan(~x)) => (log(sin(~x)) + log(tan(~x)), ∂(~x))
136-
@rule 𝛷(1 / csc(~x)) => (cos(~x) + log(csc(~x)), ∂(~x))
137-
@rule 𝛷(1 / sec(~x)) => (sin(~x) + log(sec(~x)), ∂(~x))
138-
@rule 𝛷(1 / cot(~x)) => (log(cos(~x)) + log(cot(~x)), ∂(~x))
124+
@rule 𝛷(1 / sin(~x)) => (log(csc(~x) + cot(~x)) + log(sin(~x)), ~x)
125+
@rule 𝛷(1 / cos(~x)) => (log(sec(~x) + tan(~x)) + log(cos(~x)), ~x)
126+
@rule 𝛷(1 / tan(~x)) => (log(sin(~x)) + log(tan(~x)), ~x)
127+
@rule 𝛷(1 / csc(~x)) => (cos(~x) + log(csc(~x)), ~x)
128+
@rule 𝛷(1 / sec(~x)) => (sin(~x) + log(sec(~x)), ~x)
129+
@rule 𝛷(1 / cot(~x)) => (log(cos(~x)) + log(cot(~x)), ~x)
139130
# 1/hyperbolic functions
140-
@rule 𝛷(1 / sinh(~x)) => (log(tanh(~x / 2)) + log(sinh(~x)), ∂(~x))
141-
@rule 𝛷(1 / cosh(~x)) => (atan(sinh(~x)) + log(cosh(~x)), ∂(~x))
142-
@rule 𝛷(1 / tanh(~x)) => (log(sinh(~x)) + log(tanh(~x)), ∂(~x))
143-
@rule 𝛷(1 / csch(~x)) => (cosh(~x) + log(csch(~x)), ∂(~x))
144-
@rule 𝛷(1 / sech(~x)) => (sinh(~x) + log(sech(~x)), ∂(~x))
145-
@rule 𝛷(1 / coth(~x)) => (log(cosh(~x)) + log(coth(~x)), ∂(~x))
131+
@rule 𝛷(1 / sinh(~x)) => (log(tanh(~x / 2)) + log(sinh(~x)), ~x)
132+
@rule 𝛷(1 / cosh(~x)) => (atan(sinh(~x)) + log(cosh(~x)), ~x)
133+
@rule 𝛷(1 / tanh(~x)) => (log(sinh(~x)) + log(tanh(~x)), ~x)
134+
@rule 𝛷(1 / csch(~x)) => (cosh(~x) + log(csch(~x)), ~x)
135+
@rule 𝛷(1 / sech(~x)) => (sinh(~x) + log(sech(~x)), ~x)
136+
@rule 𝛷(1 / coth(~x)) => (log(cosh(~x)) + log(coth(~x)), ~x)
146137
# inverse trigonometric functions
147-
@rule 𝛷(asin(~x)) => (~x * asin(~x) + sqrt(1 - ~x * ~x), ∂(~x))
148-
@rule 𝛷(acos(~x)) => (~x * acos(~x) + sqrt(1 - ~x * ~x), ∂(~x))
149-
@rule 𝛷(atan(~x)) => (~x * atan(~x) + log(~x * ~x + 1), ∂(~x))
150-
@rule 𝛷(acsc(~x)) => (~x * acsc(~x) + atanh(1 - ^(~x, -2)), ∂(~x))
151-
@rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ∂(~x))
152-
@rule 𝛷(acot(~x)) => (~x * acot(~x) + log(~x * ~x + 1), ∂(~x))
138+
@rule 𝛷(asin(~x)) => (~x * asin(~x) + sqrt(1 - ~x * ~x), ~x)
139+
@rule 𝛷(acos(~x)) => (~x * acos(~x) + sqrt(1 - ~x * ~x), ~x)
140+
@rule 𝛷(atan(~x)) => (~x * atan(~x) + log(~x * ~x + 1), ~x)
141+
@rule 𝛷(acsc(~x)) => (~x * acsc(~x) + atanh(1 - ^(~x, -2)), ~x)
142+
@rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ~x)
143+
@rule 𝛷(acot(~x)) => (~x * acot(~x) + log(~x * ~x + 1), ~x)
153144
# inverse hyperbolic functions
154-
@rule 𝛷(asinh(~x)) => (~x * asinh(~x) + sqrt(~x * ~x + 1), ∂(~x))
155-
@rule 𝛷(acosh(~x)) => (~x * acosh(~x) + sqrt(~x * ~x - 1), ∂(~x))
156-
@rule 𝛷(atanh(~x)) => (~x * atanh(~x) + log(~x + 1), ∂(~x))
157-
@rule 𝛷(acsch(~x)) => (acsch(~x), ∂(~x))
158-
@rule 𝛷(asech(~x)) => (asech(~x), ∂(~x))
159-
@rule 𝛷(acoth(~x)) => (~x * acot(~x) + log(~x + 1), ∂(~x))
145+
@rule 𝛷(asinh(~x)) => (~x * asinh(~x) + sqrt(~x * ~x + 1), ~x)
146+
@rule 𝛷(acosh(~x)) => (~x * acosh(~x) + sqrt(~x * ~x - 1), ~x)
147+
@rule 𝛷(atanh(~x)) => (~x * atanh(~x) + log(~x + 1), ~x)
148+
@rule 𝛷(acsch(~x)) => (acsch(~x), ~x)
149+
@rule 𝛷(asech(~x)) => (asech(~x), ~x)
150+
@rule 𝛷(acoth(~x)) => (~x * acot(~x) + log(~x + 1), ~x)
160151
# logarithmic and exponential functions
161152
@rule 𝛷(log(~x)) => (~x + ~x * log(~x) +
162153
sum(candidate_pow_minus(~x, -1); init = one(~x)),
163-
∂(~x))
164-
@rule 𝛷(1 / log(~x)) => (log(log(~x)) + li(~x), ∂(~x))
165-
@rule 𝛷(exp(~x)) => (exp(~x) + ei(~x), ∂(~x))
166-
@rule 𝛷(^(exp(~x), ~k::is_neg)) => (^(exp(-~x), -~k), ∂(~x))
154+
~x)
155+
@rule 𝛷(1 / log(~x)) => (log(log(~x)) + li(~x), ~x)
156+
@rule 𝛷(exp(~x)) => (exp(~x) + ei(~x), ~x)
157+
@rule 𝛷(^(exp(~x), ~k::is_neg)) => (^(exp(-~x), -~k), ~x)
167158
# square-root functions
168159
@rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k);
169-
init = one(~x)), 1);
170-
@rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), ∂(~x));
171-
@rule 𝛷(1 / sqrt(~x)) => (sum(candidate_sqrt(~x, -0.5); init = one(~x)), ∂(~x));
160+
init = one(~x)), ~x);
161+
@rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), ~x);
162+
@rule 𝛷(1 / sqrt(~x)) => (sum(candidate_sqrt(~x, -0.5); init = one(~x)), ~x);
172163
# rational functions
173164
@rule 𝛷(1 / ^(~x::is_poly, ~k::is_pos_int)) => (sum(candidate_pow_minus(~x, -~k);
174-
init = one(~x)), 1)
165+
init = one(~x)), ~x)
175166
@rule 𝛷(1 / ~x::is_poly) => (sum(candidate_pow_minus(~x, -1);
176-
init = one(~x)), 1)
177-
@rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
167+
init = one(~x)), ~x)
168+
@rule 𝛷(^(~x, -1)) => (log(~x), ~x)
178169
@rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1),
179-
∂(~x))
180-
@rule 𝛷(1 / ~x) => (log(~x), ∂(~x))
181-
@rule 𝛷(^(~x, ~k::is_pos_int)) => (sum(^(~x, i+1) for i=1:~k+1), ∂(~x))
170+
~x)
171+
@rule 𝛷(1 / ~x) => (log(~x), ~x)
172+
@rule 𝛷(^(~x, ~k::is_pos_int)) => (sum(^(~x, i+1) for i=1:~k+1), ~x)
182173
@rule 𝛷(1) => (𝑥, 1)
183-
@rule 𝛷(~x) => ((~x + ^(~x, 2)), ∂(~x))]
174+
@rule 𝛷(~x) => ((~x + ^(~x, 2)), ~x)]
184175

185-
function apply_partial_int_rules(eq)
186-
expand(Fixpoint(Prewalk(Chain(partial_int_rules))))(𝛷(value(eq)))
176+
function apply_partial_int_rules(eq, x)
177+
y, dy = expand(Fixpoint(Prewalk(Chain(partial_int_rules))))(𝛷(value(eq)))
178+
D = Differential(x)
179+
return y, guard_zero(expand_derivatives(D(dy)))
187180
end
181+
182+

‎src/integral.jl‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,6 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 1
3939
symbolic = true, max_basis = 100, verbose = false, complex_plane = true,
4040
homotopy = true, use_optim = false)
4141
eq = expand(eq)
42-
eq = apply_div_rule(eq)
4342

4443
if x == nothing
4544
x = var(eq)
@@ -256,9 +255,7 @@ end
256255

257256
# integrate_basis is used for debugging and should not be called in the course of normal execution
258257
function integrate_basis(eq, x = var(eq); abstol = 1e-6, radius = 1.0, complex_plane = true)
259-
eq = expand(eq)
260-
eq = apply_div_rule(eq)
261-
eq = cache(eq)
258+
eq = cache(expand(eq))
262259
basis = generate_basis(eq, x, false)
263260
n = length(basis)
264261
A, X = init_basis_matrix(eq, x, basis, radius, complex_plane; abstol)

‎src/rules.jl‎

Lines changed: 11 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -244,38 +244,34 @@ factor_rational(eq) = Prewalk(PassThrough(Chain(rational_rules)))(Ω(value(eq)))
244244

245245
###############################################################################
246246

247-
div_rule = @rule ~x / ~y => ~x * ^(~y, -1)
248-
249-
apply_div_rule(eq) = Prewalk(PassThrough(div_rule))(value(eq))
250-
251-
###############################################################################
252-
253247
inv_rules = [@rule Ω(1 / ~x) => ~x
254248
@rule Ω(~x / ~y) => Ω(~x) * ~y
255249
@rule Ω(^(~x, -1)) => ~x
256250
@rule Ω(^(~x, ~k)) => ^(Ω(~x), ~k)
257-
@rule Ω(~x) => ^(~x, -1)]
251+
@rule Ω(~x) => 1 / ~x]
258252

259253
inverse(eq) = Prewalk(PassThrough(Chain(inv_rules)))(Ω(value(eq)))
260254

261255
###############################################################################
262256

257+
@syms ω
258+
263259
is_sym(x) = first(ops(x)) isa Sym
264260

265261
h_rules = [
266-
@rule +(~~xs) => 1 + sum(~~xs)
267-
@rule *(~~xs) => 1 + sum(~~xs)
268-
@rule ~x / ~y => 1 + ~x + ~y
269-
@rule ^(~x, ~y) => 1 + ~x + ~y
270-
@rule log(~x) => 1 + ~x
271-
@rule (~f)(~x) => 1 + ~x
272-
@rule ~x::is_sym => 1
262+
@rule +(~~xs) => ω + sum(~~xs)
263+
@rule *(~~xs) => ω + sum(~~xs)
264+
@rule ~x / ~y => ω + ~x + ~y
265+
@rule ^(~x, ~y) => ω + ~x + ~y
266+
@rule (~f)(~x) => ω + ~x
267+
@rule ~x::is_sym => ω
273268
]
274269

275270
# complexity returns a measure of the complexity of an equation
276271
# it is roughly similar ro kolmogorov complexity
277272
function complexity(eq)
278273
_, eq = ops(eq)
279-
return Prewalk(PassThrough(Chain(h_rules)))(eq)
274+
h = Prewalk(PassThrough(Chain(h_rules)))(eq)
275+
return substitute(h, Dict(ω => 1))
280276
end
281277

0 commit comments

Comments
 (0)