Skip to content

perf: hoist loop-invariant substitution out of JacobianSpectrum #146

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
oameye merged 1 commit into from
Jul 26, 2026
Merged

perf: hoist loop-invariant substitution out of JacobianSpectrum #146

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "HarmonicSteadyState"
uuid = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
version = "0.5.0"
version = "0.5.1"
authors = ["Orjan Ameye <orjan.ameye@hotmail.com>", "Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]

[deps]
Expand Down
24 changes: 17 additions & 7 deletions src/LinearResponse/Lorentzian_spectrum.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,18 +105,28 @@ function JacobianSpectrum(
solutions = get_variable_solutions(res; branch=branch, index=index)
λs, vs = eigen(res.jacobian(solutions))

# The harmonic of a uv pair does not depend on the eigenvalue, so it is substituted
# once per pair instead of once per (eigenvalue, pair). The symbolic substitution and
# the `Dict` it needs dominate this loop, so hoisting them matters.
substitutions = Dict(solution_dict)
uv_pairs = _get_uv_pairs(hvars)
a_idxs = _get_as(hvars)
uv_ωnums = [
real(
SymbolicUtils.unwrap_const(
Symbolics.unwrap(Symbolics.substitute(hvars[pair][1].ω, substitutions))
),
) for pair in uv_pairs
]

for (j, λ) in enumerate(λs)
eigvec = vs[:, j] # the eigenvector

# 2 peaks for each pair of uv variables
for pair in _get_uv_pairs(hvars)
for (p, pair) in enumerate(uv_pairs)
u, v = hvars[pair]
eigvec_2d = eigvec[pair] # fetch the relevant part of the Jacobian eigenvector
ωnum = real(
SymbolicUtils.unwrap_const(
Symbolics.unwrap(Symbolics.substitute(u.ω, Dict(solution_dict)))
),
)
ωnum = uv_ωnums[p]
# ^ the harmonic (numerical now) associated to this harmonic variable

# eigvec_2d is associated to a natural variable -> this variable gets Lorentzian peaks
Expand All @@ -128,7 +138,7 @@ function JacobianSpectrum(
end

# 1 peak for a-type variable
for a_idx in _get_as(hvars)
for a_idx in a_idxs
a = hvars[a_idx]
eigvec_1d = eigvec[a_idx]
peak = 2 * norm(eigvec_1d) * Lorentzian(; ω0=abs(imag(λ)), Γ=real(λ))
Expand Down
Loading