Skip to content

Commit 46315bc

Browse files
authored
minimal symbolic clean up (#268)
1 parent 0d5d110 commit 46315bc

19 files changed

+646
-579
lines changed

ext/SteadyStateDiffEqExt.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ using HarmonicBalance: HarmonicBalance, steady_state_sweep
55

66
using SteadyStateDiffEq: solve, NonlinearProblem, SteadyStateProblem, DynamicSS, remake
77
using LinearAlgebra: norm, eigvals
8-
using SteadyStateDiffEq.SciMLBase.SciMLStructures: isscimlstructure, Tunable, replace
8+
using SteadyStateDiffEq.SciMLBase.SciMLStructures: Tunable, replace
99

1010
function HarmonicBalance.steady_state_sweep(
1111
prob::SteadyStateProblem, alg::DynamicSS; varied::Pair, kwargs...

ext/TimeEvolution/TimeEvolution.jl

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -28,10 +28,9 @@ include("ODEProblem.jl")
2828
include("hysteresis_sweep.jl")
2929
include("plotting.jl")
3030

31-
export FFT
3231
export ParameterSweep
33-
export transform_solutions, plot, plot!, is_stable
32+
export transform_solutions, plot, plot!
3433
export ODEProblem, solve
35-
export plot_1D_solutions_branch, follow_branch
34+
export follow_branch
3635

3736
end

src/HarmonicBalance.jl

Lines changed: 47 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -24,8 +24,44 @@ function set_imaginary_tolerance(x::Float64)
2424
@eval(IM_TOL::Float64 = $x)
2525
end
2626

27-
include("Symbolics_customised.jl")
28-
include("Symbolics_utils.jl")
27+
using SymbolicUtils:
28+
SymbolicUtils,
29+
Postwalk,
30+
Sym,
31+
BasicSymbolic,
32+
isterm,
33+
ispow,
34+
isadd,
35+
isdiv,
36+
ismul,
37+
add_with_div,
38+
frac_maketerm,
39+
@compactified,
40+
issym
41+
42+
using Symbolics:
43+
Symbolics,
44+
Num,
45+
unwrap,
46+
wrap,
47+
get_variables,
48+
simplify,
49+
expand_derivatives,
50+
build_function,
51+
Equation,
52+
Differential,
53+
@variables,
54+
arguments,
55+
simplify_fractions,
56+
substitute,
57+
term,
58+
expand,
59+
operation
60+
61+
include("Symbolics/Symbolics_utils.jl")
62+
include("Symbolics/exponentials.jl")
63+
include("Symbolics/fourier.jl")
64+
include("Symbolics/drop_powers.jl")
2965

3066
include("modules/extention_functions.jl")
3167
include("utils.jl")
@@ -74,14 +110,14 @@ export get_krylov_equations
74110
include("modules/FFTWExt.jl")
75111
using .FFTWExt
76112

77-
@setup_workload begin
78-
# Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
79-
# precompile file and potentially make loading faster.
80-
@compile_workload begin
81-
# all calls in this block will be precompiled, regardless of whether
82-
# they belong to your package or not (on Julia 1.8 and higher)
83-
include("precompilation.jl")
84-
end
85-
end
113+
# @setup_workload begin
114+
# # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
115+
# # precompile file and potentially make loading faster.
116+
# @compile_workload begin
117+
# # all calls in this block will be precompiled, regardless of whether
118+
# # they belong to your package or not (on Julia 1.8 and higher)
119+
# include("precompilation.jl")
120+
# end
121+
# end
86122

87123
end # module

src/HarmonicEquation.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,8 @@ function slow_flow!(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degre
6363
drop = [d(var, fast_time, degree) => 0 for var in get_variables(eom)]
6464

6565
eom.equations = substitute_all(substitute_all(eom.equations, drop), replace)
66-
return eom.variables = substitute_all(eom.variables, replace)
66+
eom.variables = substitute_all(eom.variables, replace)
67+
return nothing
6768
end
6869

6970
"""

src/HarmonicVariable.jl

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,3 +51,33 @@ Symbolics.get_variables(var::HarmonicVariable)::Num = Num(first(get_variables(va
5151

5252
Base.isequal(v1::HarmonicVariable, v2::HarmonicVariable)::Bool =
5353
isequal(v1.symbol, v2.symbol)
54+
55+
"The derivative of f w.r.t. x of degree deg"
56+
function d(f::Num, x::Num, deg=1)::Num
57+
return isequal(deg, 0) ? f : (Differential(x)^deg)(f)
58+
end
59+
d(funcs::Vector{Num}, x::Num, deg=1) = Num[d(f, x, deg) for f in funcs]
60+
61+
"Declare a variable in the the currently active namespace"
62+
function declare_variable(name::String)
63+
var_sym = Symbol(name)
64+
@eval($(var_sym) = first(Symbolics.@variables $var_sym))
65+
return eval(var_sym)
66+
end
67+
68+
"Declare a variable that is a function of another variable in the the current namespace"
69+
function declare_variable(name::String, independent_variable::Num)
70+
# independent_variable = declare_variable(independent_variable) convert string into Num
71+
var_sym = Symbol(name)
72+
new_var = Symbolics.@variables $var_sym(independent_variable)
73+
@eval($(var_sym) = first($new_var)) # store the variable under "name" in this namespace
74+
return eval(var_sym)
75+
end
76+
77+
"Return the name of a variable (excluding independent variables)"
78+
function var_name(x::Num)
79+
var = Symbolics._toexpr(x)
80+
return var isa Expr ? String(var.args[1]) : String(var)
81+
end
82+
# var_name(x::Term) = String(Symbolics._toexpr(x).args[1])
83+
var_name(x::Sym) = String(x.name)

src/Symbolics/Symbolics_utils.jl

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
2+
expand_all(x::Num) = Num(expand_all(x.val))
3+
_apply_termwise(f, x::Num) = wrap(_apply_termwise(f, unwrap(x)))
4+
5+
"Expands using SymbolicUtils.expand and expand_exp_power (changes exp(x)^n to exp(x*n)"
6+
function expand_all(x)
7+
result = Postwalk(expand_exp_power)(SymbolicUtils.expand(x))
8+
return isnothing(result) ? x : result
9+
end
10+
expand_all(x::Complex{Num}) = expand_all(x.re) + im * expand_all(x.im)
11+
12+
"Apply a function f on every member of a sum or a product"
13+
function _apply_termwise(f, x::BasicSymbolic)
14+
@compactified x::BasicSymbolic begin
15+
Add => sum([f(arg) for arg in arguments(x)])
16+
Mul => prod([f(arg) for arg in arguments(x)])
17+
Div => _apply_termwise(f, x.num) / _apply_termwise(f, x.den)
18+
_ => f(x)
19+
end
20+
end
21+
22+
simplify_complex(x::Complex) = isequal(x.im, 0) ? x.re : x.re + im * x.im
23+
simplify_complex(x) = x
24+
function simplify_complex(x::BasicSymbolic)
25+
@compactified x::BasicSymbolic begin
26+
Add => _apply_termwise(simplify_complex, x)
27+
Mul => _apply_termwise(simplify_complex, x)
28+
Div => _apply_termwise(simplify_complex, x)
29+
_ => x
30+
end
31+
end
32+
33+
"""
34+
$(TYPEDSIGNATURES)
35+
36+
Perform substitutions in `rules` on `x`.
37+
`include_derivatives=true` also includes all derivatives of the variables of the keys of `rules`.
38+
"""
39+
subtype = Union{Num,Equation,BasicSymbolic}
40+
function substitute_all(x::subtype, rules::Dict; include_derivatives=true)
41+
if include_derivatives
42+
rules = merge(
43+
rules,
44+
Dict([Differential(var) => Differential(rules[var]) for var in keys(rules)]),
45+
)
46+
end
47+
return substitute(x, rules)
48+
end
49+
"Variable substitution - dictionary"
50+
function substitute_all(dict::Dict, rules::Dict)::Dict
51+
new_keys = substitute_all.(keys(dict), rules)
52+
new_values = substitute_all.(values(dict), rules)
53+
return Dict(zip(new_keys, new_values))
54+
end
55+
Collections = Union{Dict,Pair,Vector,OrderedDict}
56+
substitute_all(v::AbstractArray, rules) = [substitute_all(x, rules) for x in v]
57+
substitute_all(x::subtype, rules::Collections) = substitute_all(x, Dict(rules))
58+
# Collections = Union{Dict,OrderedDict}
59+
# function substitute_all(x, rules::Collections; include_derivatives=true)
60+
# if include_derivatives
61+
# rules = merge(
62+
# rules,
63+
# Dict([Differential(var) => Differential(rules[var]) for var in keys(rules)]),
64+
# )
65+
# end
66+
# return substitute(x, rules)
67+
# end
68+
# "Variable substitution - dictionary"
69+
# function substitute_all(dict::Dict, rules::Dict)::Dict
70+
# new_keys = substitute_all.(keys(dict), rules)
71+
# new_values = substitute_all.(values(dict), rules)
72+
# return Dict(zip(new_keys, new_values))
73+
# end
74+
# substitute_all(v::AbstractArray, rules::Collections) = [substitute_all(x, rules) for x in v]
75+
76+
get_independent(x::Num, t::Num) = get_independent(x.val, t)
77+
function get_independent(x::Complex{Num}, t::Num)
78+
return get_independent(x.re, t) + im * get_independent(x.im, t)
79+
end
80+
get_independent(v::Vector{Num}, t::Num) = [get_independent(el, t) for el in v]
81+
get_independent(x, t::Num) = x
82+
83+
function get_independent(x::BasicSymbolic, t::Num)
84+
@compactified x::BasicSymbolic begin
85+
Add => sum([get_independent(arg, t) for arg in arguments(x)])
86+
Mul => prod([get_independent(arg, t) for arg in arguments(x)])
87+
Div => !is_function(x.den, t) ? get_independent(x.num, t) / x.den : 0
88+
Pow => !is_function(x.base, t) && !is_function(x.exp, t) ? x : 0
89+
Term => !is_function(x, t) ? x : 0
90+
Sym => !is_function(x, t) ? x : 0
91+
_ => x
92+
end
93+
end
94+
95+
"Return all the terms contained in `x`"
96+
get_all_terms(x::Num) = unique(_get_all_terms(Symbolics.expand(x).val))
97+
function get_all_terms(x::Equation)
98+
return unique(cat(get_all_terms(Num(x.lhs)), get_all_terms(Num(x.rhs)); dims=1))
99+
end
100+
function _get_all_terms(x::BasicSymbolic)
101+
@compactified x::BasicSymbolic begin
102+
Add => vcat([_get_all_terms(term) for term in SymbolicUtils.arguments(x)]...)
103+
Mul => Num.(SymbolicUtils.arguments(x))
104+
Div => Num.([_get_all_terms(x.num)..., _get_all_terms(x.den)...])
105+
_ => Num(x)
106+
end
107+
end
108+
_get_all_terms(x) = Num(x)
109+
110+
function is_harmonic(x::Num, t::Num)::Bool
111+
all_terms = get_all_terms(x)
112+
t_terms = setdiff(all_terms, get_independent(all_terms, t))
113+
isempty(t_terms) && return true
114+
trigs = is_trig.(t_terms)
115+
116+
if !prod(trigs)
117+
return false
118+
else
119+
powers = [max_power(first(term.val.arguments), t) for term in t_terms[trigs]]
120+
return all(isone, powers)
121+
end
122+
end
123+
124+
is_harmonic(x::Equation, t::Num) = is_harmonic(x.lhs, t) && is_harmonic(x.rhs, t)
125+
is_harmonic(x, t) = is_harmonic(Num(x), Num(t))
126+
127+
"Return true if `f` is a function of `var`."
128+
is_function(f, var) = any(isequal.(get_variables(f), var))

src/Symbolics/drop_powers.jl

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
"""
2+
$(SIGNATURES)
3+
Remove parts of `expr` where the combined power of `vars` is => `deg`.
4+
5+
# Example
6+
```julia-repl
7+
julia> @variables x,y;
8+
julia>drop_powers((x+y)^2, x, 2)
9+
y^2 + 2*x*y
10+
julia>drop_powers((x+y)^2, [x,y], 2)
11+
0
12+
julia>drop_powers((x+y)^2 + (x+y)^3, [x,y], 3)
13+
x^2 + y^2 + 2*x*y
14+
```
15+
"""
16+
function drop_powers(expr::Num, vars::Vector{Num}, deg::Int)
17+
Symbolics.@variables ϵ
18+
subs_expr = deepcopy(expr)
19+
rules = Dict([var => ϵ * var for var in unique(vars)])
20+
subs_expr = Symbolics.expand(substitute_all(subs_expr, rules))
21+
max_deg = max_power(subs_expr, ϵ)
22+
removal = Dict([ϵ^d => Num(0) for d in deg:max_deg])
23+
res = substitute_all(substitute_all(subs_expr, removal), Dict=> Num(1)))
24+
return Symbolics.expand(res)
25+
end
26+
27+
function drop_powers(expr::Vector{Num}, var::Vector{Num}, deg::Int)
28+
return [drop_powers(x, var, deg) for x in expr]
29+
end
30+
31+
# calls the above for various types of the first argument
32+
function drop_powers(eq::Equation, var::Vector{Num}, deg::Int)
33+
return drop_powers(eq.lhs, var, deg) .~ drop_powers(eq.lhs, var, deg)
34+
end
35+
function drop_powers(eqs::Vector{Equation}, var::Vector{Num}, deg::Int)
36+
return [
37+
Equation(drop_powers(eq.lhs, var, deg), drop_powers(eq.rhs, var, deg)) for eq in eqs
38+
]
39+
end
40+
drop_powers(expr, var::Num, deg::Int) = drop_powers(expr, [var], deg)
41+
drop_powers(x, vars, deg::Int) = drop_powers(Num(x), vars, deg)
42+
43+
"Return the highest power of `y` occuring in the term `x`."
44+
function max_power(x::Num, y::Num)
45+
terms = get_all_terms(x)
46+
powers = power_of.(terms, y)
47+
return maximum(powers)
48+
end
49+
50+
max_power(x::Vector{Num}, y::Num) = maximum(max_power.(x, y))
51+
max_power(x::Complex, y::Num) = maximum(max_power.([x.re, x.im], y))
52+
max_power(x, t) = max_power(Num(x), Num(t))
53+
54+
"Return the power of `y` in the term `x`"
55+
function power_of(x::Num, y::Num)
56+
issym(y.val) ? nothing : error("power of " * string(y) * " is ambiguous")
57+
return power_of(x.val, y.val)
58+
end
59+
60+
function power_of(x::BasicSymbolic, y::BasicSymbolic)
61+
if ispow(x) && issym(y)
62+
return isequal(x.base, y) ? x.exp : 0
63+
elseif issym(x) && issym(y)
64+
return isequal(x, y) ? 1 : 0
65+
else
66+
return 0
67+
end
68+
end
69+
70+
power_of(x, y) = 0

src/Symbolics/exponentials.jl

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
expand_exp_power(expr::Num) = expand_exp_power(expr.val)
2+
simplify_exp_products(x::Num) = simplify_exp_products(x.val)
3+
4+
"Returns true if expr is an exponential"
5+
isexp(expr) = isterm(expr) && expr.f == exp
6+
7+
"Expand powers of exponential such that exp(x)^n => exp(x*n) "
8+
function expand_exp_power(expr::BasicSymbolic)
9+
@compactified expr::BasicSymbolic begin
10+
Add => sum([expand_exp_power(arg) for arg in arguments(expr)])
11+
Mul => prod([expand_exp_power(arg) for arg in arguments(expr)])
12+
_ => ispow(expr) && isexp(expr.base) ? exp(expr.base.arguments[1] * expr.exp) : expr
13+
end
14+
end
15+
expand_exp_power(expr) = expr
16+
17+
"Simplify products of exponentials such that exp(a)*exp(b) => exp(a+b)
18+
This is included in SymbolicUtils as of 17.0 but the method here avoid other simplify calls"
19+
function simplify_exp_products(expr::BasicSymbolic)
20+
@compactified expr::BasicSymbolic begin
21+
Add => _apply_termwise(simplify_exp_products, expr)
22+
Div => _apply_termwise(simplify_exp_products, expr)
23+
Mul => simplify_exp_products_mul(expr)
24+
_ => expr
25+
end
26+
end
27+
function simplify_exp_products(x::Complex{Num})
28+
return Complex{Num}(simplify_exp_products(x.re.val), simplify_exp_products(x.im.val))
29+
end
30+
function simplify_exp_products_mul(expr)
31+
ind = findall(x -> isexp(x), arguments(expr))
32+
rest_ind = setdiff(1:length(arguments(expr)), ind)
33+
rest = isempty(rest_ind) ? 1 : prod(arguments(expr)[rest_ind])
34+
total = isempty(ind) ? 0 : sum(getindex.(arguments.(arguments(expr)[ind]), 1))
35+
if SymbolicUtils.is_literal_number(total)
36+
(total == 0 && return rest)
37+
else
38+
return rest * exp(total)
39+
end
40+
end
41+
simplify_exp_products(x) = x

0 commit comments

Comments
 (0)