From 172de6f1b5f77878a15d3ade11b6b25f694b0df5 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 2 Jun 2026 12:10:02 +0000 Subject: [PATCH 1/7] Add populate_cache for batched multi-RHS warming of virtual matrices Introduce `populate_cache` for VirtualPTDF, VirtualLODF, and VirtualMODF so callers can pre-warm row caches over an iterable of components/contingencies using the solver backends' batched multi-RHS `solve_sparse!` instead of one single-RHS solve per row. This amortizes the ABA factorization solves that dominate optimization-problem builds on large systems. - VirtualPTDF / VirtualLODF: gather requested arc rows, build one sparse RHS of the corresponding BA columns, solve in a single batched call, then scatter (PTDF) or apply the LODF map (A*X .* inv_PTDF_A_diag) across all columns at once. Accepts integer indices, arc bus-pair tuples, or branch-name strings. - VirtualMODF: batch the pre-contingency BA-column solves over the union of all monitored arcs and every contingency's modified arcs. The monitored-arc solve is contingency-independent, so it is computed once and reused; Woodbury factors per contingency are assembled from the precomputed solves with no further linear solves. Monitored set is supplied by the user. - RowCache: add set_persistent_row!/pin_row! so populated rows are pinned and never evicted by later lazy lookups, plus warn_if_over_capacity. - Backend dispatch routes to KLU's solve_sparse! or AccelerateWrapper's, with a column-by-column fallback for other factorization backends. - Export populate_cache; add tests covering PTDF/LODF/MODF equivalence with the lazy getindex path, pinning, identifier resolution, and error handling. https://claude.ai/code/session_011TFwa4Hsdse4XYEeuHmYro --- src/PowerNetworkMatrices.jl | 2 + src/populate_cache.jl | 390 ++++++++++++++++++++++++++++++++++++ src/row_cache.jl | 57 ++++++ test/test_populate_cache.jl | 136 +++++++++++++ 4 files changed, 585 insertions(+) create mode 100644 src/populate_cache.jl create mode 100644 test/test_populate_cache.jl diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl index f2b4151fb..fbc760e45 100644 --- a/src/PowerNetworkMatrices.jl +++ b/src/PowerNetworkMatrices.jl @@ -41,6 +41,7 @@ export factorize export get_post_modification_ptdf_row export get_system_uuid export is_factorized +export populate_cache export depth_first_search export find_subnetworks @@ -147,6 +148,7 @@ include("woodbury_kernel.jl") include("virtual_ptdf_modification.jl") include("modf_reduction_consistency.jl") include("virtual_modf_calculations.jl") +include("populate_cache.jl") include("system_utils.jl") include("serialization.jl") diff --git a/src/populate_cache.jl b/src/populate_cache.jl new file mode 100644 index 000000000..2f0d1d570 --- /dev/null +++ b/src/populate_cache.jl @@ -0,0 +1,390 @@ +# Bulk cache population for the lazy virtual network matrices. +# +# `getindex` on a `VirtualPTDF` / `VirtualLODF` / `VirtualMODF` computes one row +# at a time, issuing a single right-hand-side (RHS) linear solve per row through +# the ABA factorization. When the optimization-problem build queries many rows +# (every monitored branch, every contingency), those solves dominate. +# +# `populate_cache` amortizes that cost: it gathers the requested rows, builds a +# single sparse RHS matrix of the corresponding `BA` columns, and solves them +# all in one batched call via the backend's `solve_sparse!` (KLU or Apple +# Accelerate). A batched solve reuses the factorization's working set across +# columns and issues one libklu/libSparse call per block of 64 instead of one +# per row, which is markedly faster on large systems. The resulting rows are +# pinned in the row cache so later `getindex` calls are pure cache hits. + +# --- Backend dispatch: batched (multi-RHS) sparse solve -------------------- + +# KLU's `solve_sparse!` is imported unqualified into PNM; Apple Accelerate's is +# reached through the `AccelerateWrapper` module (mirrors `ptdf_calculations.jl`). +_solve_multi_rhs!( + K::KLULinSolveCache{Float64}, + B::SparseArrays.SparseMatrixCSC{Float64, Int}, + out::Matrix{Float64}, +) = solve_sparse!(K, B, out) + +_solve_multi_rhs!( + K::AAFactorCache, + B::SparseArrays.SparseMatrixCSC{Float64, Int}, + out::Matrix{Float64}, +) = AccelerateWrapper.solve_sparse!(K, B, out) + +# Generic fallback for any other factorization backend: solve column-by-column. +function _solve_multi_rhs!(K, B::SparseArrays.SparseMatrixCSC{Float64, Int}, out::Matrix{Float64}) + n = size(B, 1) + col = zeros(n) + @inbounds for j in axes(B, 2) + fill!(col, 0.0) + for p in SparseArrays.nzrange(B, j) + col[SparseArrays.rowvals(B)[p]] = SparseArrays.nonzeros(B)[p] + end + _solve_factorization(K, col) + copyto!(view(out, :, j), col) + end + return out +end + +""" + _solve_arc_columns(mat, arc_rows) -> Matrix{Float64} + +Solve `ABA · X = BA[valid_ix, arc_rows]` for all requested arcs at once, +returning the `(n_valid × length(arc_rows))` dense solution. The caller must +hold `mat.solver_lock` (the batched solve mutates the factorization's scratch). +""" +function _solve_arc_columns(mat, arc_rows::Vector{Int}) + valid_ix = mat.valid_ix + B = mat.BA[valid_ix, arc_rows] + out = Matrix{Float64}(undef, length(valid_ix), length(arc_rows)) + _solve_multi_rhs!(mat.K, B, out) + return out +end + +# --- Component resolution -------------------------------------------------- + +# Resolve a user-supplied component to an integer arc/row index. Accepts the +# same identifiers as `getindex`: an integer index, an arc bus-pair tuple, or a +# branch name. +_resolve_arc_index(::PowerNetworkMatrix, c::Integer) = Int(c) +_resolve_arc_index(mat::PowerNetworkMatrix, c::Tuple{Int, Int}) = get_arc_lookup(mat)[c] +function _resolve_arc_index(mat::PowerNetworkMatrix, c::AbstractString) + _, arc = get_branch_multiplier(mat, String(c)) + return get_arc_lookup(mat)[arc] +end + +# --- VirtualPTDF ----------------------------------------------------------- + +""" + populate_cache(vptdf::VirtualPTDF, components) -> Nothing + +Precompute and pin the PTDF rows for an iterable of `components`, using a single +batched multi-RHS solve instead of one solve per row. `components` may mix +integer arc indices, arc bus-pair tuples `(from, to)`, and branch-name strings. + +Populated rows are added to the cache's persistent set, so subsequent +`vptdf[component, :]` queries are cache hits and are never evicted by later lazy +lookups. Use this before an optimization-problem build to amortize the linear +solves over the monitored branch set. + +$(TYPEDSIGNATURES) +""" +function populate_cache(vptdf::VirtualPTDF, components) + rows = unique(Int[_resolve_arc_index(vptdf, c) for c in components]) + isempty(rows) && return nothing + + buscount = size(vptdf, 1) + ref_bus_positions = get_ref_bus_position(vptdf) + if !isempty(vptdf.dist_slack) && length(ref_bus_positions) != 1 + error( + "Distributed slack is not supported for systems with multiple reference buses.", + ) + end + use_dist_slack = length(vptdf.dist_slack) == buscount + if !use_dist_slack && !isempty(vptdf.dist_slack) + error("Distributed bus specification doesn't match the number of buses.") + end + tol = get_tol(vptdf) + + # Only solve rows not already resident; existing rows are pinned below. + new_rows = @lock vptdf.cache_lock Int[r for r in rows if !haskey(vptdf.cache, r)] + + if !isempty(new_rows) + sol = @lock vptdf.solver_lock _solve_arc_columns(vptdf, new_rows) + valid_ix = vptdf.valid_ix + @lock vptdf.cache_lock begin + full = zeros(buscount) + for (j, r) in enumerate(new_rows) + haskey(vptdf.cache, r) && continue # lost a race; keep the winner + fill!(full, 0.0) + @inbounds for i in eachindex(valid_ix) + full[valid_ix[i]] = sol[i, j] + end + if use_dist_slack + full .-= dot(full, vptdf.dist_slack_normalized) + end + stored = tol > eps() ? sparsify(full, tol) : copy(full) + set_persistent_row!(vptdf.cache, r, stored) + end + end + end + + @lock vptdf.cache_lock begin + for r in rows + pin_row!(vptdf.cache, r) + end + warn_if_over_capacity(vptdf.cache) + end + return nothing +end + +# --- VirtualLODF ----------------------------------------------------------- + +""" + populate_cache(vlodf::VirtualLODF, components) -> Nothing + +Precompute and pin the LODF rows for an iterable of `components` (integer arc +indices, arc bus-pair tuples, or branch-name strings) using a single batched +multi-RHS solve. The post-contingency scaling `(A · B⁻¹ · BA) .* inv_PTDF_A_diag` +is applied to all requested rows at once via one sparse-dense product. + +Populated rows are pinned in the cache so later `vlodf[component, :]` queries +are cache hits. + +$(TYPEDSIGNATURES) +""" +function populate_cache(vlodf::VirtualLODF, components) + rows = unique(Int[_resolve_arc_index(vlodf, c) for c in components]) + isempty(rows) && return nothing + tol = get_tol(vlodf) + n_bus = length(vlodf.temp_data[1]) + + new_rows = @lock vlodf.cache_lock Int[r for r in rows if !haskey(vlodf.cache, r)] + + if !isempty(new_rows) + sol = @lock vlodf.solver_lock _solve_arc_columns(vlodf, new_rows) + valid_ix = vlodf.valid_ix + # Scatter each solved column back to full-bus space, then apply the + # LODF map to every column at once: L = (A · Tmp) .* inv_PTDF_A_diag. + tmp = zeros(n_bus, length(new_rows)) + @inbounds for j in eachindex(new_rows), i in eachindex(valid_ix) + tmp[valid_ix[i], j] = sol[i, j] + end + lodf_cols = vlodf.A * tmp # (n_arcs × length(new_rows)) + lodf_cols .*= vlodf.inv_PTDF_A_diag # broadcast per-arc scaling down columns + @inbounds for (j, r) in enumerate(new_rows) + lodf_cols[r, j] = -1.0 # self-element convention + end + @lock vlodf.cache_lock begin + for (j, r) in enumerate(new_rows) + haskey(vlodf.cache, r) && continue + row = lodf_cols[:, j] + stored = tol > eps() ? sparsify(row, tol) : row + set_persistent_row!(vlodf.cache, r, stored) + end + end + end + + @lock vlodf.cache_lock begin + for r in rows + pin_row!(vlodf.cache, r) + end + warn_if_over_capacity(vlodf.cache) + end + return nothing +end + +# --- VirtualMODF ----------------------------------------------------------- + +# Resolve a contingency identifier to the `NetworkModification` used as the +# cache key. Accepts a modification, a `ContingencySpec`, a registered +# `PSY.Outage`, or a registered outage UUID. +_resolve_modification(::VirtualMODF, mod::NetworkModification) = mod +_resolve_modification(::VirtualMODF, ctg::ContingencySpec) = ctg.modification +function _resolve_modification(vmodf::VirtualMODF, outage::PSY.Outage) + return _resolve_modification(vmodf, IS.get_uuid(outage)) +end +function _resolve_modification(vmodf::VirtualMODF, uuid::Base.UUID) + haskey(vmodf.contingency_cache, uuid) || error( + "Contingency (UUID=$uuid) is not registered. Construct the VirtualMODF " * + "with the system containing this outage, or pass the NetworkModification " * + "/ ContingencySpec directly.", + ) + return vmodf.contingency_cache[uuid].modification +end + +""" + _woodbury_factors_from_base(base_full, BA, arc_sus, modifications, n_bus) -> WoodburyFactors + +Assemble Woodbury factors from precomputed pre-contingency solves. `base_full` +maps each modified arc index to `B⁻¹ · BA[:, arc]` scattered to full-bus space. +Identical math to `_compute_woodbury_factors_impl` but with the per-arc libklu +solves replaced by dictionary lookups (the batched solve already paid for them). +""" +function _woodbury_factors_from_base( + base_full::Dict{Int, Vector{Float64}}, + BA::SparseArrays.SparseMatrixCSC{Float64, Int}, + arc_sus::Vector{Float64}, + modifications::Tuple{Vararg{ArcModification}}, + n_bus::Int, +)::WoodburyFactors + M = length(modifications) + arc_indices = Vector{Int}(undef, M) + delta_b_vec = Vector{Float64}(undef, M) + for (j, mod) in enumerate(modifications) + arc_indices[j] = mod.arc_index + delta_b_vec[j] = mod.delta_b + end + + # Z[:, j] = B⁻¹ ν_j = (B⁻¹ BA[:, e_j]) / b_{e_j} + Z = Matrix{Float64}(undef, n_bus, M) + for j in 1:M + b_e = arc_sus[arc_indices[j]] + col = base_full[arc_indices[j]] + @inbounds for i in 1:n_bus + Z[i, j] = col[i] / b_e + end + end + + # K_mat[i, j] = ν_i⊤ B⁻¹ ν_j, iterating the sparse BA columns of arc e_i. + ba_nzv = SparseArrays.nonzeros(BA) + ba_rv = SparseArrays.rowvals(BA) + K_mat = zeros(M, M) + for i in 1:M + e_i = arc_indices[i] + b_i = arc_sus[e_i] + for j in 1:M + val = 0.0 + @inbounds for nz_idx in SparseArrays.nzrange(BA, e_i) + val += (ba_nzv[nz_idx] / b_i) * Z[ba_rv[nz_idx], j] + end + K_mat[i, j] = val + end + end + + W_mat = LinearAlgebra.diagm(1.0 ./ delta_b_vec) + K_mat + W_inv, is_island = _invert_woodbury_W(W_mat, Val(M)) + return WoodburyFactors(Z, W_inv, arc_indices, delta_b_vec, is_island) +end + +""" + _woodbury_correction_from_base(base_full, BA, arc_sus, monitored_idx, wf, n_bus) -> Vector{Float64} + +Post-modification PTDF row for `monitored_idx` from precomputed pre-contingency +solves. Mirrors `_apply_woodbury_correction_impl`, reusing `base_full[monitored_idx]` +(contingency-independent) instead of solving. +""" +function _woodbury_correction_from_base( + base_full::Dict{Int, Vector{Float64}}, + BA::SparseArrays.SparseMatrixCSC{Float64, Int}, + arc_sus::Vector{Float64}, + monitored_idx::Int, + wf::WoodburyFactors, + n_bus::Int, +)::Vector{Float64} + M = length(wf.arc_indices) + + b_mon = arc_sus[monitored_idx] + for (j, idx) in enumerate(wf.arc_indices) + idx == monitored_idx && (b_mon += wf.delta_b[j]) + end + abs(b_mon) < eps() && return zeros(n_bus) + + b_mon_pre = arc_sus[monitored_idx] + temp = base_full[monitored_idx] ./ b_mon_pre # z_m (fresh vector; base_full untouched) + + ba_nzv = SparseArrays.nonzeros(BA) + ba_rv = SparseArrays.rowvals(BA) + zm_Z = zeros(M) + @inbounds for nz_idx in SparseArrays.nzrange(BA, monitored_idx) + coeff = ba_nzv[nz_idx] / b_mon_pre + row = ba_rv[nz_idx] + for j in 1:M + zm_Z[j] += coeff * wf.Z[row, j] + end + end + + correction_coeff = wf.W_inv * zm_Z + LinearAlgebra.mul!(temp, wf.Z, correction_coeff, -1.0, 1.0) + temp .*= b_mon + return temp +end + +""" + populate_cache(vmodf::VirtualMODF, contingencies; monitored) -> Nothing + +Precompute and pin the post-contingency PTDF rows for an iterable of +`contingencies`, each evaluated over the user-supplied `monitored` arc set, using +batched multi-RHS solves. + +`contingencies` may mix `NetworkModification`, `ContingencySpec`, registered +`PSY.Outage`, and registered outage UUIDs. `monitored` is an iterable of arc +identifiers (integer indices, bus-pair tuples, or branch names). + +The acceleration comes from a single batched solve over the union of all arcs +referenced — every contingency's modified arcs plus every monitored arc. Because +the pre-contingency solve for a monitored arc is contingency-independent, it is +computed once and reused across all contingencies; the Woodbury factors per +contingency are then assembled from these precomputed solves with no further +linear solves. This replaces the `O(n_contingencies × (M + n_monitored))` +one-at-a-time solves of repeated `getindex` with one batched solve over the +distinct arcs. + +Results are written into the per-contingency row caches (and Woodbury cache) and +pinned, so subsequent `vmodf[monitored, contingency]` queries are cache hits. + +$(TYPEDSIGNATURES) +""" +function populate_cache(vmodf::VirtualMODF, contingencies; monitored) + mods = unique(NetworkModification[_resolve_modification(vmodf, c) for c in contingencies]) + mon_idx = unique(Int[_resolve_arc_index(vmodf, m) for m in monitored]) + (isempty(mods) || isempty(mon_idx)) && return nothing + + n_bus = length(vmodf.temp_data[1]) + tol = get_tol(vmodf) + BA = vmodf.BA + arc_sus = vmodf.arc_susceptances + + @lock vmodf.solver_lock begin + # Union of all arcs needing a pre-contingency solve: monitored arcs and + # every contingency's modified arcs. + arc_set = Set{Int}(mon_idx) + for mod in mods + for am in mod.arc_modifications + push!(arc_set, am.arc_index) + end + end + all_arcs = collect(arc_set) + + # One batched solve for B⁻¹ BA[:, arc] over the distinct arcs. + sol = _solve_arc_columns(vmodf, all_arcs) + valid_ix = vmodf.valid_ix + base_full = Dict{Int, Vector{Float64}}() + sizehint!(base_full, length(all_arcs)) + for (j, arc) in enumerate(all_arcs) + full = zeros(n_bus) + @inbounds for i in eachindex(valid_ix) + full[valid_ix[i]] = sol[i, j] + end + base_full[arc] = full + end + + for mod in mods + wf = get!(vmodf.woodbury_cache, mod) do + _woodbury_factors_from_base(base_full, BA, arc_sus, mod.arc_modifications, n_bus) + end + rc = get!(vmodf.row_caches, mod) do + RowCache(vmodf.max_cache_size_bytes, Set{Int}(), n_bus * sizeof(Float64)) + end + for m in mon_idx + if haskey(rc, m) + pin_row!(rc, m) + continue + end + row = _woodbury_correction_from_base(base_full, BA, arc_sus, m, wf, n_bus) + stored = tol > eps() ? sparsify(row, tol) : row + set_persistent_row!(rc, m, stored) + end + warn_if_over_capacity(rc) + end + end + return nothing +end diff --git a/src/row_cache.jl b/src/row_cache.jl index 2d09ca086..31423721c 100644 --- a/src/row_cache.jl +++ b/src/row_cache.jl @@ -128,6 +128,63 @@ function Base.getindex( return cache.temp_cache[key] end +""" +Stores `val` for `key` and pins `key` so it is never evicted by LRU. + +Used by `populate_cache` to bulk-fill rows computed via multi-RHS solves and +guarantee they stay warm for later queries. Unlike `setindex!`, this never +triggers `purge_one!`: pinned rows are added to `persistent_cache_keys`, and a +cache populated beyond `max_num_keys` simply grows (callers should warn via +[`warn_if_over_capacity`](@ref)). + +# Arguments +- `cache::RowCache`: + cache where the row vector is stored and pinned. +- `key::Int`: + row number (enumerated branch index) for the row vector. +- `val`: + the row vector (dense `Vector{Float64}` or sparsified `SparseVector{Float64}`). +""" +function set_persistent_row!( + cache::RowCache{T}, + key::Int, + val::T, +) where {T <: Union{Vector{Float64}, SparseArrays.SparseVector{Float64}}} + if !haskey(cache.temp_cache, key) + push!(cache.access_order, key) + end + cache.temp_cache[key] = val + push!(cache.persistent_cache_keys, key) + return +end + +""" +Pin an already-stored `key` so a row populated lazily is also protected from +eviction. No-op for keys absent from `temp_cache`. +""" +function pin_row!(cache::RowCache, key::Int) + haskey(cache.temp_cache, key) && push!(cache.persistent_cache_keys, key) + return +end + +""" + warn_if_over_capacity(cache::RowCache) + +Emit a single warning when the cache holds more rows than `max_num_keys`. This +happens when `populate_cache` pins more rows than the configured +`max_cache_size` can hold; the pinned rows remain resident (never evicted), so +the only risk is higher memory use. +""" +function warn_if_over_capacity(cache::RowCache) + if length(cache.temp_cache) > cache.max_num_keys + @warn "populate_cache pinned $(length(cache.temp_cache)) rows, exceeding " * + "the cache capacity (max_num_keys = $(cache.max_num_keys)). Pinned " * + "rows are not evicted; increase `max_cache_size` to avoid memory " * + "pressure." maxlog = 1 + end + return +end + """ Shows the number of rows stored in cache """ diff --git a/test/test_populate_cache.jl b/test/test_populate_cache.jl new file mode 100644 index 000000000..e9eba387c --- /dev/null +++ b/test/test_populate_cache.jl @@ -0,0 +1,136 @@ +@testset "populate_cache VirtualPTDF: batched solve matches lazy getindex" for solver in + ("KLU", "AppleAccelerate") + if !PowerNetworkMatrices._has_apple_accelerate_backend() && solver == "AppleAccelerate" + @info "Skipped AppleAccelerate populate_cache tests (backend unavailable)" + continue + end + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + v_lazy = VirtualPTDF(sys; linear_solver = solver) + v_pop = VirtualPTDF(sys; linear_solver = solver) + + n_arcs = length(v_pop.axes[1]) + arc_idxs = collect(1:min(6, n_arcs)) + + @test populate_cache(v_pop, arc_idxs) === nothing + + for i in arc_idxs + # Row is resident and pinned after population + @test haskey(v_pop.cache.temp_cache, i) + @test i ∈ v_pop.cache.persistent_cache_keys + # Batched multi-RHS solve must agree with one-at-a-time lazy compute + @test isapprox(v_pop[i, :], v_lazy[i, :]; atol = 1e-10) + end +end + +@testset "populate_cache VirtualPTDF: tuple and string identifiers" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + v_lazy = VirtualPTDF(sys) + v_pop = VirtualPTDF(sys) + + arc_tuples = collect(keys(v_pop.lookup[1]))[1:4] + populate_cache(v_pop, arc_tuples) + for a in arc_tuples + idx = v_pop.lookup[1][a] + @test idx ∈ v_pop.cache.persistent_cache_keys + @test isapprox(v_pop[a, :], v_lazy[a, :]; atol = 1e-10) + end +end + +@testset "populate_cache VirtualPTDF: pins previously-lazy rows" begin + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + v = VirtualPTDF(sys) + # Warm a row lazily (not pinned yet) + _ = v[1, :] + @test haskey(v.cache.temp_cache, 1) + @test 1 ∉ v.cache.persistent_cache_keys + # populate_cache must pin it without recomputing incorrectly + populate_cache(v, [1, 2]) + @test 1 ∈ v.cache.persistent_cache_keys + @test 2 ∈ v.cache.persistent_cache_keys +end + +@testset "populate_cache VirtualLODF: batched solve matches lazy getindex" for solver in + ("KLU", "AppleAccelerate") + if !PowerNetworkMatrices._has_apple_accelerate_backend() && solver == "AppleAccelerate" + @info "Skipped AppleAccelerate populate_cache tests (backend unavailable)" + continue + end + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + v_lazy = VirtualLODF(sys; linear_solver = solver) + v_pop = VirtualLODF(sys; linear_solver = solver) + + n_arcs = length(v_pop.axes[1]) + arc_idxs = collect(1:min(6, n_arcs)) + + populate_cache(v_pop, arc_idxs) + for i in arc_idxs + @test haskey(v_pop.cache.temp_cache, i) + @test i ∈ v_pop.cache.persistent_cache_keys + @test isapprox(v_pop[i, :], v_lazy[i, :]; atol = 1e-10) + # Self-element convention preserved by the batched path + @test v_pop[i, i] == -1.0 + end +end + +@testset "populate_cache VirtualMODF: batched contingencies match lazy getindex" begin + sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") + v_lazy = VirtualMODF(sys5) + v_pop = VirtualMODF(sys5) + + n_arcs = size(v_pop, 1) + monitored = collect(1:min(4, n_arcs)) + + ctgs = ContingencySpec[] + for e in 1:3 + b_e = v_pop.arc_susceptances[e] + uuid = Base.UUID(UInt128(7000 + e)) + ctg = ContingencySpec( + uuid, + NetworkModification("populate_ctg_$e", [ArcModification(e, -b_e)]), + ) + v_pop.contingency_cache[uuid] = ctg + v_lazy.contingency_cache[uuid] = ctg + push!(ctgs, ctg) + end + + @test populate_cache(v_pop, ctgs; monitored = monitored) === nothing + + for ctg in ctgs + mod = ctg.modification + @test haskey(v_pop.woodbury_cache, mod) + @test haskey(v_pop.row_caches, mod) + rc = v_pop.row_caches[mod] + for m in monitored + @test haskey(rc, m) + @test m ∈ rc.persistent_cache_keys + @test isapprox(v_pop[m, ctg], v_lazy[m, ctg]; atol = 1e-8) + end + end +end + +@testset "populate_cache VirtualMODF: tuple-monitored and outage resolution errors" begin + sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") + vmodf = VirtualMODF(sys5) + + e = 1 + b_e = vmodf.arc_susceptances[e] + uuid = Base.UUID(UInt128(7100)) + ctg = ContingencySpec( + uuid, + NetworkModification("populate_tuple_ctg", [ArcModification(e, -b_e)]), + ) + vmodf.contingency_cache[uuid] = ctg + + # Monitor by arc bus-pair tuple + mon_tuple = vmodf.axes[1][2] + populate_cache(vmodf, [ctg]; monitored = [mon_tuple]) + m_idx = vmodf.lookup[1][mon_tuple] + @test haskey(vmodf.row_caches[ctg.modification], m_idx) + + # Unregistered UUID must error clearly + @test_throws ErrorException populate_cache( + vmodf, + [Base.UUID(UInt128(99999))]; + monitored = [1], + ) +end From 002b3d082d32344022685b06f69dab1ecf36c5cb Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 2 Jun 2026 15:16:17 +0000 Subject: [PATCH 2/7] Adapt populate_cache to psy6 shared-core backend; use getters not getproperty Merge of origin/psy6 moved the virtual matrices' shared state into VirtualFactorCore. Rework populate_cache so every wrapper field is read through its proper getter (get_core, get_cache, get_cache_lock, get_dist_slack, get_dist_slack_normalized, get_inv_PTDF_A_diag, get_row_caches, get_woodbury_cache, get_max_cache_size_bytes, get_contingency_cache, get_arc_lookup, get_tol, get_ref_bus_position) instead of relying on the wrappers' getproperty forwarding. Core container fields are accessed directly on the VirtualFactorCore, matching _compute_ptdf_row/_compute_lodf_row and the Woodbury kernel. The batched multi-RHS solve and the from-precomputed-solve Woodbury assembly are unchanged and remain byte-for-byte consistent with the post-merge woodbury_kernel.jl. Tests read cache/core/contingency state through getters too. --- src/populate_cache.jl | 115 +++++++++++++++++++++--------------- test/test_populate_cache.jl | 44 +++++++------- 2 files changed, 88 insertions(+), 71 deletions(-) diff --git a/src/populate_cache.jl b/src/populate_cache.jl index 2f0d1d570..4b80b3988 100644 --- a/src/populate_cache.jl +++ b/src/populate_cache.jl @@ -2,8 +2,9 @@ # # `getindex` on a `VirtualPTDF` / `VirtualLODF` / `VirtualMODF` computes one row # at a time, issuing a single right-hand-side (RHS) linear solve per row through -# the ABA factorization. When the optimization-problem build queries many rows -# (every monitored branch, every contingency), those solves dominate. +# the shared `VirtualFactorCore` ABA factorization. When the optimization-problem +# build queries many rows (every monitored branch, every contingency), those +# solves dominate. # # `populate_cache` amortizes that cost: it gathers the requested rows, builds a # single sparse RHS matrix of the corresponding `BA` columns, and solves them @@ -45,17 +46,17 @@ function _solve_multi_rhs!(K, B::SparseArrays.SparseMatrixCSC{Float64, Int}, out end """ - _solve_arc_columns(mat, arc_rows) -> Matrix{Float64} + _solve_arc_columns(core, arc_rows) -> Matrix{Float64} Solve `ABA · X = BA[valid_ix, arc_rows]` for all requested arcs at once, returning the `(n_valid × length(arc_rows))` dense solution. The caller must -hold `mat.solver_lock` (the batched solve mutates the factorization's scratch). +hold `core.solver_lock` (the batched solve mutates the factorization's scratch). """ -function _solve_arc_columns(mat, arc_rows::Vector{Int}) - valid_ix = mat.valid_ix - B = mat.BA[valid_ix, arc_rows] +function _solve_arc_columns(core::VirtualFactorCore, arc_rows::Vector{Int}) + valid_ix = core.valid_ix + B = core.BA[valid_ix, arc_rows] out = Matrix{Float64}(undef, length(valid_ix), length(arc_rows)) - _solve_multi_rhs!(mat.K, B, out) + _solve_multi_rhs!(core.K, B, out) return out end @@ -91,47 +92,53 @@ function populate_cache(vptdf::VirtualPTDF, components) rows = unique(Int[_resolve_arc_index(vptdf, c) for c in components]) isempty(rows) && return nothing - buscount = size(vptdf, 1) - ref_bus_positions = get_ref_bus_position(vptdf) - if !isempty(vptdf.dist_slack) && length(ref_bus_positions) != 1 + core = get_core(vptdf) + cache = get_cache(vptdf) + cache_lock = get_cache_lock(vptdf) + dist_slack = get_dist_slack(vptdf) + dist_slack_normalized = get_dist_slack_normalized(vptdf) + + buscount = size(core.BA, 1) + ref_bus_positions = get_ref_bus_position(core) + if !isempty(dist_slack) && length(ref_bus_positions) != 1 error( "Distributed slack is not supported for systems with multiple reference buses.", ) end - use_dist_slack = length(vptdf.dist_slack) == buscount - if !use_dist_slack && !isempty(vptdf.dist_slack) + use_dist_slack = length(dist_slack) == buscount + if !use_dist_slack && !isempty(dist_slack) error("Distributed bus specification doesn't match the number of buses.") end - tol = get_tol(vptdf) + tol = get_tol(core) # Only solve rows not already resident; existing rows are pinned below. - new_rows = @lock vptdf.cache_lock Int[r for r in rows if !haskey(vptdf.cache, r)] + new_rows = @lock cache_lock Int[r for r in rows if !haskey(cache, r)] if !isempty(new_rows) - sol = @lock vptdf.solver_lock _solve_arc_columns(vptdf, new_rows) - valid_ix = vptdf.valid_ix - @lock vptdf.cache_lock begin + sol = @lock core.solver_lock _solve_arc_columns(core, new_rows) + valid_ix = core.valid_ix + @lock cache_lock begin full = zeros(buscount) for (j, r) in enumerate(new_rows) - haskey(vptdf.cache, r) && continue # lost a race; keep the winner + haskey(cache, r) && continue # lost a race; keep the winner fill!(full, 0.0) @inbounds for i in eachindex(valid_ix) full[valid_ix[i]] = sol[i, j] end if use_dist_slack - full .-= dot(full, vptdf.dist_slack_normalized) + full .-= dot(full, dist_slack_normalized) end stored = tol > eps() ? sparsify(full, tol) : copy(full) - set_persistent_row!(vptdf.cache, r, stored) + set_persistent_row!(cache, r, stored) end end end - @lock vptdf.cache_lock begin + @lock cache_lock begin for r in rows - pin_row!(vptdf.cache, r) + pin_row!(cache, r) end - warn_if_over_capacity(vptdf.cache) + warn_if_over_capacity(cache) end return nothing end @@ -154,40 +161,45 @@ $(TYPEDSIGNATURES) function populate_cache(vlodf::VirtualLODF, components) rows = unique(Int[_resolve_arc_index(vlodf, c) for c in components]) isempty(rows) && return nothing - tol = get_tol(vlodf) - n_bus = length(vlodf.temp_data[1]) - new_rows = @lock vlodf.cache_lock Int[r for r in rows if !haskey(vlodf.cache, r)] + core = get_core(vlodf) + cache = get_cache(vlodf) + cache_lock = get_cache_lock(vlodf) + inv_PTDF_A_diag = get_inv_PTDF_A_diag(vlodf) + tol = get_tol(core) + n_bus = length(core.temp_data[1]) + + new_rows = @lock cache_lock Int[r for r in rows if !haskey(cache, r)] if !isempty(new_rows) - sol = @lock vlodf.solver_lock _solve_arc_columns(vlodf, new_rows) - valid_ix = vlodf.valid_ix + sol = @lock core.solver_lock _solve_arc_columns(core, new_rows) + valid_ix = core.valid_ix # Scatter each solved column back to full-bus space, then apply the # LODF map to every column at once: L = (A · Tmp) .* inv_PTDF_A_diag. tmp = zeros(n_bus, length(new_rows)) @inbounds for j in eachindex(new_rows), i in eachindex(valid_ix) tmp[valid_ix[i], j] = sol[i, j] end - lodf_cols = vlodf.A * tmp # (n_arcs × length(new_rows)) - lodf_cols .*= vlodf.inv_PTDF_A_diag # broadcast per-arc scaling down columns + lodf_cols = core.A * tmp # (n_arcs × length(new_rows)) + lodf_cols .*= inv_PTDF_A_diag # broadcast per-arc scaling down columns @inbounds for (j, r) in enumerate(new_rows) lodf_cols[r, j] = -1.0 # self-element convention end - @lock vlodf.cache_lock begin + @lock cache_lock begin for (j, r) in enumerate(new_rows) - haskey(vlodf.cache, r) && continue + haskey(cache, r) && continue row = lodf_cols[:, j] stored = tol > eps() ? sparsify(row, tol) : row - set_persistent_row!(vlodf.cache, r, stored) + set_persistent_row!(cache, r, stored) end end end - @lock vlodf.cache_lock begin + @lock cache_lock begin for r in rows - pin_row!(vlodf.cache, r) + pin_row!(cache, r) end - warn_if_over_capacity(vlodf.cache) + warn_if_over_capacity(cache) end return nothing end @@ -203,12 +215,13 @@ function _resolve_modification(vmodf::VirtualMODF, outage::PSY.Outage) return _resolve_modification(vmodf, IS.get_uuid(outage)) end function _resolve_modification(vmodf::VirtualMODF, uuid::Base.UUID) - haskey(vmodf.contingency_cache, uuid) || error( + contingency_cache = get_contingency_cache(vmodf) + haskey(contingency_cache, uuid) || error( "Contingency (UUID=$uuid) is not registered. Construct the VirtualMODF " * "with the system containing this outage, or pass the NetworkModification " * "/ ContingencySpec directly.", ) - return vmodf.contingency_cache[uuid].modification + return contingency_cache[uuid].modification end """ @@ -338,12 +351,16 @@ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) mon_idx = unique(Int[_resolve_arc_index(vmodf, m) for m in monitored]) (isempty(mods) || isempty(mon_idx)) && return nothing - n_bus = length(vmodf.temp_data[1]) - tol = get_tol(vmodf) - BA = vmodf.BA - arc_sus = vmodf.arc_susceptances + core = get_core(vmodf) + row_caches = get_row_caches(vmodf) + woodbury_cache = get_woodbury_cache(vmodf) + max_bytes = get_max_cache_size_bytes(vmodf) + n_bus = length(core.temp_data[1]) + tol = get_tol(core) + BA = core.BA + arc_sus = core.arc_susceptances - @lock vmodf.solver_lock begin + @lock core.solver_lock begin # Union of all arcs needing a pre-contingency solve: monitored arcs and # every contingency's modified arcs. arc_set = Set{Int}(mon_idx) @@ -355,8 +372,8 @@ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) all_arcs = collect(arc_set) # One batched solve for B⁻¹ BA[:, arc] over the distinct arcs. - sol = _solve_arc_columns(vmodf, all_arcs) - valid_ix = vmodf.valid_ix + sol = _solve_arc_columns(core, all_arcs) + valid_ix = core.valid_ix base_full = Dict{Int, Vector{Float64}}() sizehint!(base_full, length(all_arcs)) for (j, arc) in enumerate(all_arcs) @@ -368,11 +385,11 @@ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) end for mod in mods - wf = get!(vmodf.woodbury_cache, mod) do + wf = get!(woodbury_cache, mod) do _woodbury_factors_from_base(base_full, BA, arc_sus, mod.arc_modifications, n_bus) end - rc = get!(vmodf.row_caches, mod) do - RowCache(vmodf.max_cache_size_bytes, Set{Int}(), n_bus * sizeof(Float64)) + rc = get!(row_caches, mod) do + RowCache(max_bytes, Set{Int}(), n_bus * sizeof(Float64)) end for m in mon_idx if haskey(rc, m) diff --git a/test/test_populate_cache.jl b/test/test_populate_cache.jl index e9eba387c..630bf581d 100644 --- a/test/test_populate_cache.jl +++ b/test/test_populate_cache.jl @@ -8,15 +8,15 @@ v_lazy = VirtualPTDF(sys; linear_solver = solver) v_pop = VirtualPTDF(sys; linear_solver = solver) - n_arcs = length(v_pop.axes[1]) + n_arcs = length(PNM.get_arc_axis(v_pop)) arc_idxs = collect(1:min(6, n_arcs)) @test populate_cache(v_pop, arc_idxs) === nothing for i in arc_idxs # Row is resident and pinned after population - @test haskey(v_pop.cache.temp_cache, i) - @test i ∈ v_pop.cache.persistent_cache_keys + @test haskey(PNM.get_cache(v_pop).temp_cache, i) + @test i ∈ PNM.get_cache(v_pop).persistent_cache_keys # Batched multi-RHS solve must agree with one-at-a-time lazy compute @test isapprox(v_pop[i, :], v_lazy[i, :]; atol = 1e-10) end @@ -27,11 +27,11 @@ end v_lazy = VirtualPTDF(sys) v_pop = VirtualPTDF(sys) - arc_tuples = collect(keys(v_pop.lookup[1]))[1:4] + arc_tuples = collect(keys(PNM.get_arc_lookup(v_pop)))[1:4] populate_cache(v_pop, arc_tuples) for a in arc_tuples - idx = v_pop.lookup[1][a] - @test idx ∈ v_pop.cache.persistent_cache_keys + idx = PNM.get_arc_lookup(v_pop)[a] + @test idx ∈ PNM.get_cache(v_pop).persistent_cache_keys @test isapprox(v_pop[a, :], v_lazy[a, :]; atol = 1e-10) end end @@ -41,12 +41,12 @@ end v = VirtualPTDF(sys) # Warm a row lazily (not pinned yet) _ = v[1, :] - @test haskey(v.cache.temp_cache, 1) - @test 1 ∉ v.cache.persistent_cache_keys + @test haskey(PNM.get_cache(v).temp_cache, 1) + @test 1 ∉ PNM.get_cache(v).persistent_cache_keys # populate_cache must pin it without recomputing incorrectly populate_cache(v, [1, 2]) - @test 1 ∈ v.cache.persistent_cache_keys - @test 2 ∈ v.cache.persistent_cache_keys + @test 1 ∈ PNM.get_cache(v).persistent_cache_keys + @test 2 ∈ PNM.get_cache(v).persistent_cache_keys end @testset "populate_cache VirtualLODF: batched solve matches lazy getindex" for solver in @@ -59,13 +59,13 @@ end v_lazy = VirtualLODF(sys; linear_solver = solver) v_pop = VirtualLODF(sys; linear_solver = solver) - n_arcs = length(v_pop.axes[1]) + n_arcs = length(PNM.get_arc_axis(v_pop)) arc_idxs = collect(1:min(6, n_arcs)) populate_cache(v_pop, arc_idxs) for i in arc_idxs - @test haskey(v_pop.cache.temp_cache, i) - @test i ∈ v_pop.cache.persistent_cache_keys + @test haskey(PNM.get_cache(v_pop).temp_cache, i) + @test i ∈ PNM.get_cache(v_pop).persistent_cache_keys @test isapprox(v_pop[i, :], v_lazy[i, :]; atol = 1e-10) # Self-element convention preserved by the batched path @test v_pop[i, i] == -1.0 @@ -88,8 +88,8 @@ end uuid, NetworkModification("populate_ctg_$e", [ArcModification(e, -b_e)]), ) - v_pop.contingency_cache[uuid] = ctg - v_lazy.contingency_cache[uuid] = ctg + PNM.get_contingency_cache(v_pop)[uuid] = ctg + PNM.get_contingency_cache(v_lazy)[uuid] = ctg push!(ctgs, ctg) end @@ -97,9 +97,9 @@ end for ctg in ctgs mod = ctg.modification - @test haskey(v_pop.woodbury_cache, mod) - @test haskey(v_pop.row_caches, mod) - rc = v_pop.row_caches[mod] + @test haskey(PNM.get_woodbury_cache(v_pop), mod) + @test haskey(PNM.get_row_caches(v_pop), mod) + rc = PNM.get_row_caches(v_pop)[mod] for m in monitored @test haskey(rc, m) @test m ∈ rc.persistent_cache_keys @@ -119,13 +119,13 @@ end uuid, NetworkModification("populate_tuple_ctg", [ArcModification(e, -b_e)]), ) - vmodf.contingency_cache[uuid] = ctg + PNM.get_contingency_cache(vmodf)[uuid] = ctg # Monitor by arc bus-pair tuple - mon_tuple = vmodf.axes[1][2] + mon_tuple = PNM.get_arc_axis(vmodf)[2] populate_cache(vmodf, [ctg]; monitored = [mon_tuple]) - m_idx = vmodf.lookup[1][mon_tuple] - @test haskey(vmodf.row_caches[ctg.modification], m_idx) + m_idx = PNM.get_arc_lookup(vmodf)[mon_tuple] + @test haskey(PNM.get_row_caches(vmodf)[ctg.modification], m_idx) # Unregistered UUID must error clearly @test_throws ErrorException populate_cache( From a0deb0b64faf39b3b4fc2f0ef53a4e6967154701 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 2 Jun 2026 15:21:25 +0000 Subject: [PATCH 3/7] Fix populate_cache tests for psy6 indexing; apply formatter wrapping The redesigned backend routes array getindex through to_index, which resolves the row against the arc-tuple lookup, so vptdf[i, :] with an integer row no longer works (integer+integer and arc-tuple+colon do). Rework the PTDF/LODF tests to index rows by arc tuple, and cover integer-keyed population by inspecting the row cache directly. VirtualMODF keeps its integer getindex(::Int, contingency) method, so those tests are unchanged in spirit. Also apply the JuliaFormatter wrapping flagged by reviewdog on src/populate_cache.jl (generic _solve_multi_rhs!, the mods comprehension, and the _woodbury_factors_from_base call) and restructure the solver-looped testsets to a begin/for form that stays within the line-length margin. --- src/populate_cache.jl | 17 +++++- test/test_populate_cache.jl | 116 ++++++++++++++++++++---------------- 2 files changed, 77 insertions(+), 56 deletions(-) diff --git a/src/populate_cache.jl b/src/populate_cache.jl index 4b80b3988..47be2e464 100644 --- a/src/populate_cache.jl +++ b/src/populate_cache.jl @@ -31,7 +31,11 @@ _solve_multi_rhs!( ) = AccelerateWrapper.solve_sparse!(K, B, out) # Generic fallback for any other factorization backend: solve column-by-column. -function _solve_multi_rhs!(K, B::SparseArrays.SparseMatrixCSC{Float64, Int}, out::Matrix{Float64}) +function _solve_multi_rhs!( + K, + B::SparseArrays.SparseMatrixCSC{Float64, Int}, + out::Matrix{Float64}, +) n = size(B, 1) col = zeros(n) @inbounds for j in axes(B, 2) @@ -347,7 +351,8 @@ pinned, so subsequent `vmodf[monitored, contingency]` queries are cache hits. $(TYPEDSIGNATURES) """ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) - mods = unique(NetworkModification[_resolve_modification(vmodf, c) for c in contingencies]) + mods = + unique(NetworkModification[_resolve_modification(vmodf, c) for c in contingencies]) mon_idx = unique(Int[_resolve_arc_index(vmodf, m) for m in monitored]) (isempty(mods) || isempty(mon_idx)) && return nothing @@ -386,7 +391,13 @@ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) for mod in mods wf = get!(woodbury_cache, mod) do - _woodbury_factors_from_base(base_full, BA, arc_sus, mod.arc_modifications, n_bus) + _woodbury_factors_from_base( + base_full, + BA, + arc_sus, + mod.arc_modifications, + n_bus, + ) end rc = get!(row_caches, mod) do RowCache(max_bytes, Set{Int}(), n_bus * sizeof(Float64)) diff --git a/test/test_populate_cache.jl b/test/test_populate_cache.jl index 630bf581d..17e985b91 100644 --- a/test/test_populate_cache.jl +++ b/test/test_populate_cache.jl @@ -1,74 +1,84 @@ -@testset "populate_cache VirtualPTDF: batched solve matches lazy getindex" for solver in - ("KLU", "AppleAccelerate") - if !PowerNetworkMatrices._has_apple_accelerate_backend() && solver == "AppleAccelerate" - @info "Skipped AppleAccelerate populate_cache tests (backend unavailable)" - continue - end - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") - v_lazy = VirtualPTDF(sys; linear_solver = solver) - v_pop = VirtualPTDF(sys; linear_solver = solver) - - n_arcs = length(PNM.get_arc_axis(v_pop)) - arc_idxs = collect(1:min(6, n_arcs)) - - @test populate_cache(v_pop, arc_idxs) === nothing - - for i in arc_idxs - # Row is resident and pinned after population - @test haskey(PNM.get_cache(v_pop).temp_cache, i) - @test i ∈ PNM.get_cache(v_pop).persistent_cache_keys - # Batched multi-RHS solve must agree with one-at-a-time lazy compute - @test isapprox(v_pop[i, :], v_lazy[i, :]; atol = 1e-10) +@testset "populate_cache VirtualPTDF: batched solve matches lazy getindex" begin + for solver in ("KLU", "AppleAccelerate") + if !PowerNetworkMatrices._has_apple_accelerate_backend() && + solver == "AppleAccelerate" + @info "Skipped AppleAccelerate populate_cache tests (backend unavailable)" + continue + end + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + v_lazy = VirtualPTDF(sys; linear_solver = solver) + v_pop = VirtualPTDF(sys; linear_solver = solver) + + arc_lookup = PNM.get_arc_lookup(v_pop) + arc_tuples = PNM.get_arc_axis(v_pop)[1:min(6, length(arc_lookup))] + + @test populate_cache(v_pop, arc_tuples) === nothing + + for a in arc_tuples + idx = arc_lookup[a] + # Row is resident and pinned after population + @test haskey(PNM.get_cache(v_pop).temp_cache, idx) + @test idx ∈ PNM.get_cache(v_pop).persistent_cache_keys + # Batched multi-RHS solve must agree with one-at-a-time lazy compute + @test isapprox(v_pop[a, :], v_lazy[a, :]; atol = 1e-10) + end end end -@testset "populate_cache VirtualPTDF: tuple and string identifiers" begin +@testset "populate_cache VirtualPTDF: integer indices populate the cache" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") v_lazy = VirtualPTDF(sys) v_pop = VirtualPTDF(sys) - arc_tuples = collect(keys(PNM.get_arc_lookup(v_pop)))[1:4] - populate_cache(v_pop, arc_tuples) - for a in arc_tuples - idx = PNM.get_arc_lookup(v_pop)[a] + arc_lookup = PNM.get_arc_lookup(v_pop) + arcs = PNM.get_arc_axis(v_pop)[1:4] + idxs = [arc_lookup[a] for a in arcs] + + populate_cache(v_pop, idxs) + for a in arcs + idx = arc_lookup[a] @test idx ∈ PNM.get_cache(v_pop).persistent_cache_keys - @test isapprox(v_pop[a, :], v_lazy[a, :]; atol = 1e-10) + @test isapprox(PNM.get_cache(v_pop).temp_cache[idx], v_lazy[a, :]; atol = 1e-10) end end @testset "populate_cache VirtualPTDF: pins previously-lazy rows" begin sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") v = VirtualPTDF(sys) + arc = PNM.get_arc_axis(v)[1] + idx = PNM.get_arc_lookup(v)[arc] # Warm a row lazily (not pinned yet) - _ = v[1, :] - @test haskey(PNM.get_cache(v).temp_cache, 1) - @test 1 ∉ PNM.get_cache(v).persistent_cache_keys + _ = v[arc, :] + @test haskey(PNM.get_cache(v).temp_cache, idx) + @test idx ∉ PNM.get_cache(v).persistent_cache_keys # populate_cache must pin it without recomputing incorrectly - populate_cache(v, [1, 2]) - @test 1 ∈ PNM.get_cache(v).persistent_cache_keys - @test 2 ∈ PNM.get_cache(v).persistent_cache_keys + populate_cache(v, [arc]) + @test idx ∈ PNM.get_cache(v).persistent_cache_keys end -@testset "populate_cache VirtualLODF: batched solve matches lazy getindex" for solver in - ("KLU", "AppleAccelerate") - if !PowerNetworkMatrices._has_apple_accelerate_backend() && solver == "AppleAccelerate" - @info "Skipped AppleAccelerate populate_cache tests (backend unavailable)" - continue - end - sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") - v_lazy = VirtualLODF(sys; linear_solver = solver) - v_pop = VirtualLODF(sys; linear_solver = solver) - - n_arcs = length(PNM.get_arc_axis(v_pop)) - arc_idxs = collect(1:min(6, n_arcs)) - - populate_cache(v_pop, arc_idxs) - for i in arc_idxs - @test haskey(PNM.get_cache(v_pop).temp_cache, i) - @test i ∈ PNM.get_cache(v_pop).persistent_cache_keys - @test isapprox(v_pop[i, :], v_lazy[i, :]; atol = 1e-10) - # Self-element convention preserved by the batched path - @test v_pop[i, i] == -1.0 +@testset "populate_cache VirtualLODF: batched solve matches lazy getindex" begin + for solver in ("KLU", "AppleAccelerate") + if !PowerNetworkMatrices._has_apple_accelerate_backend() && + solver == "AppleAccelerate" + @info "Skipped AppleAccelerate populate_cache tests (backend unavailable)" + continue + end + sys = PSB.build_system(PSB.PSITestSystems, "c_sys14") + v_lazy = VirtualLODF(sys; linear_solver = solver) + v_pop = VirtualLODF(sys; linear_solver = solver) + + arc_lookup = PNM.get_arc_lookup(v_pop) + arc_tuples = PNM.get_arc_axis(v_pop)[1:min(6, length(arc_lookup))] + + populate_cache(v_pop, arc_tuples) + for a in arc_tuples + idx = arc_lookup[a] + @test haskey(PNM.get_cache(v_pop).temp_cache, idx) + @test idx ∈ PNM.get_cache(v_pop).persistent_cache_keys + @test isapprox(v_pop[a, :], v_lazy[a, :]; atol = 1e-10) + # Self-element convention preserved by the batched path + @test v_pop[a, a] == -1.0 + end end end From 6632a37c3761507087dd0f132962e766357b958f Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 2 Jun 2026 15:24:28 +0000 Subject: [PATCH 4/7] Address Copilot review: monitored-arc error, lock scope, backend fallback - VirtualMODF populate_cache resolves monitored arc tuples through _monitored_arc_index, so an arc eliminated by network reduction yields the descriptive VirtualMODF error instead of a raw KeyError. - VirtualPTDF/VirtualLODF populate_cache build each stored row (scatter, dist-slack, sparsify) outside cache_lock and take the lock only for the double-check + insert, matching the cached_row_lookup pattern and shortening the critical section. - The generic _solve_multi_rhs! fallback now checks `applicable` and raises a clear, actionable error when a backend implements neither solve_sparse! nor _solve_factorization, instead of surfacing a raw MethodError. - Add shared RowCacheValue alias for the cached-row value type. --- src/populate_cache.jl | 56 +++++++++++++++++++++++++++++++------------ src/row_cache.jl | 4 ++++ 2 files changed, 45 insertions(+), 15 deletions(-) diff --git a/src/populate_cache.jl b/src/populate_cache.jl index 47be2e464..1df8ae3d8 100644 --- a/src/populate_cache.jl +++ b/src/populate_cache.jl @@ -30,7 +30,11 @@ _solve_multi_rhs!( out::Matrix{Float64}, ) = AccelerateWrapper.solve_sparse!(K, B, out) -# Generic fallback for any other factorization backend: solve column-by-column. +# Generic fallback for any other factorization backend: solve column-by-column +# through the single-RHS `_solve_factorization` seam (the documented extension +# point for new solver backends). Errors clearly when the backend implements +# neither a batched `solve_sparse!` nor `_solve_factorization`, instead of +# surfacing a raw MethodError. function _solve_multi_rhs!( K, B::SparseArrays.SparseMatrixCSC{Float64, Int}, @@ -38,6 +42,11 @@ function _solve_multi_rhs!( ) n = size(B, 1) col = zeros(n) + applicable(_solve_factorization, K, col) || error( + "Factorization backend $(typeof(K)) supports neither a batched " * + "`solve_sparse!` nor a single-RHS `_solve_factorization`; extend one of " * + "them to use `populate_cache` with this backend.", + ) @inbounds for j in axes(B, 2) fill!(col, 0.0) for p in SparseArrays.nzrange(B, j) @@ -121,19 +130,25 @@ function populate_cache(vptdf::VirtualPTDF, components) if !isempty(new_rows) sol = @lock core.solver_lock _solve_arc_columns(core, new_rows) valid_ix = core.valid_ix + # Build each row outside the cache lock (the scatter, dist-slack, and + # sparsify dominate); take cache_lock only for the double-check + insert, + # matching the `cached_row_lookup` pattern. + stored_rows = Vector{RowCacheValue}(undef, length(new_rows)) + full = zeros(buscount) + for (j, _) in enumerate(new_rows) + fill!(full, 0.0) + @inbounds for i in eachindex(valid_ix) + full[valid_ix[i]] = sol[i, j] + end + if use_dist_slack + full .-= dot(full, dist_slack_normalized) + end + stored_rows[j] = tol > eps() ? sparsify(full, tol) : copy(full) + end @lock cache_lock begin - full = zeros(buscount) for (j, r) in enumerate(new_rows) haskey(cache, r) && continue # lost a race; keep the winner - fill!(full, 0.0) - @inbounds for i in eachindex(valid_ix) - full[valid_ix[i]] = sol[i, j] - end - if use_dist_slack - full .-= dot(full, dist_slack_normalized) - end - stored = tol > eps() ? sparsify(full, tol) : copy(full) - set_persistent_row!(cache, r, stored) + set_persistent_row!(cache, r, stored_rows[j]) end end end @@ -189,12 +204,16 @@ function populate_cache(vlodf::VirtualLODF, components) @inbounds for (j, r) in enumerate(new_rows) lodf_cols[r, j] = -1.0 # self-element convention end + # Build each row outside the cache lock; take cache_lock only to insert. + stored_rows = Vector{RowCacheValue}(undef, length(new_rows)) + for j in eachindex(new_rows) + row = lodf_cols[:, j] + stored_rows[j] = tol > eps() ? sparsify(row, tol) : row + end @lock cache_lock begin for (j, r) in enumerate(new_rows) haskey(cache, r) && continue - row = lodf_cols[:, j] - stored = tol > eps() ? sparsify(row, tol) : row - set_persistent_row!(cache, r, stored) + set_persistent_row!(cache, r, stored_rows[j]) end end end @@ -228,6 +247,13 @@ function _resolve_modification(vmodf::VirtualMODF, uuid::Base.UUID) return contingency_cache[uuid].modification end +# Resolve a monitored arc identifier to its row index. Tuples go through +# `_monitored_arc_index` so an arc that was reduced away yields VirtualMODF's +# descriptive error instead of a raw KeyError. +_resolve_monitored_index(::VirtualMODF, m::Integer) = Int(m) +_resolve_monitored_index(vmodf::VirtualMODF, m::Tuple{Int, Int}) = + _monitored_arc_index(vmodf, m) + """ _woodbury_factors_from_base(base_full, BA, arc_sus, modifications, n_bus) -> WoodburyFactors @@ -353,7 +379,7 @@ $(TYPEDSIGNATURES) function populate_cache(vmodf::VirtualMODF, contingencies; monitored) mods = unique(NetworkModification[_resolve_modification(vmodf, c) for c in contingencies]) - mon_idx = unique(Int[_resolve_arc_index(vmodf, m) for m in monitored]) + mon_idx = unique(Int[_resolve_monitored_index(vmodf, m) for m in monitored]) (isempty(mods) || isempty(mon_idx)) && return nothing core = get_core(vmodf) diff --git a/src/row_cache.jl b/src/row_cache.jl index 31423721c..aee09dc3d 100644 --- a/src/row_cache.jl +++ b/src/row_cache.jl @@ -1,3 +1,7 @@ +# Value type stored per cached row: a dense row, or a sparsified row when a +# tolerance is applied. +const RowCacheValue = Union{Vector{Float64}, SparseArrays.SparseVector{Float64}} + """ Structure used for saving the rows of the Virtual PTDF and LODF matrix. From f2e442f32b45eeffeab598e145b8a064568683b0 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 2 Jun 2026 22:36:03 +0000 Subject: [PATCH 5/7] Use _get_arc_susceptances getter in tests instead of getproperty Replace the remaining `vmodf.arc_susceptances` field reads in the populate_cache tests with the `_get_arc_susceptances` accessor, so no code in this change reaches struct fields through the wrappers' getproperty forwarding. --- test/test_populate_cache.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_populate_cache.jl b/test/test_populate_cache.jl index 17e985b91..be9296591 100644 --- a/test/test_populate_cache.jl +++ b/test/test_populate_cache.jl @@ -92,7 +92,7 @@ end ctgs = ContingencySpec[] for e in 1:3 - b_e = v_pop.arc_susceptances[e] + b_e = PNM._get_arc_susceptances(v_pop)[e] uuid = Base.UUID(UInt128(7000 + e)) ctg = ContingencySpec( uuid, @@ -123,7 +123,7 @@ end vmodf = VirtualMODF(sys5) e = 1 - b_e = vmodf.arc_susceptances[e] + b_e = PNM._get_arc_susceptances(vmodf)[e] uuid = Base.UUID(UInt128(7100)) ctg = ContingencySpec( uuid, From f442dd7e73e9be22d2984ec5ac9948c1ed3bd5f3 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 09:26:16 +0000 Subject: [PATCH 6/7] Test populate_cache against the shared VirtualFactorCore MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add a test that backs a VirtualPTDF, VirtualLODF, and VirtualMODF with one shared VirtualFactorCore (the common ABA factorization) and runs populate_cache on each. Asserts the populated rows match independently-built matrices and that populate_cache reuses the single shared factorization object (get_core(...).K === core.K) rather than rebuilding it — confirming the feature is built on the common factorization backend introduced by the psy6 redesign. --- test/test_populate_cache.jl | 49 +++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/test/test_populate_cache.jl b/test/test_populate_cache.jl index be9296591..f948bcd74 100644 --- a/test/test_populate_cache.jl +++ b/test/test_populate_cache.jl @@ -144,3 +144,52 @@ end monitored = [1], ) end + +@testset "populate_cache: shared VirtualFactorCore (common ABA factorization)" begin + sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") + + # One factorization (VirtualFactorCore) shared across PTDF/LODF/MODF wrappers. + vptdf = VirtualPTDF(Ybus(sys5)) + vlodf = VirtualLODF(vptdf) + vmodf = VirtualMODF(vptdf, sys5; automatically_register_outages = false) + core = PNM.get_core(vptdf) + @test PNM.get_core(vlodf) === core + @test PNM.get_core(vmodf) === core + + # Independent references — each builds its own factorization. + vptdf_i = VirtualPTDF(sys5) + vlodf_i = VirtualLODF(sys5) + vmodf_i = VirtualMODF(sys5; automatically_register_outages = false) + + arc_lookup = PNM.get_arc_lookup(vptdf) + arcs = PNM.get_arc_axis(vptdf)[1:min(4, length(arc_lookup))] + + # Warm PTDF and LODF rows through the common factorization. + populate_cache(vptdf, arcs) + populate_cache(vlodf, arcs) + for a in arcs + @test isapprox(vptdf[a, :], vptdf_i[a, :]; atol = 1e-10) + @test isapprox(vlodf[a, :], vlodf_i[a, :]; atol = 1e-10) + end + + # MODF on the shared core: register one contingency in both and populate. + e = 1 + b_e = PNM._get_arc_susceptances(vmodf)[e] + uuid = Base.UUID(UInt128(7300)) + ctg = ContingencySpec( + uuid, + NetworkModification("shared_core_ctg", [ArcModification(e, -b_e)]), + ) + PNM.get_contingency_cache(vmodf)[uuid] = ctg + PNM.get_contingency_cache(vmodf_i)[uuid] = ctg + + monitored = collect(1:min(4, size(vmodf, 1))) + populate_cache(vmodf, [ctg]; monitored = monitored) + for m in monitored + @test isapprox(vmodf[m, ctg], vmodf_i[m, ctg]; atol = 1e-8) + end + + # populate_cache reused the single shared factorization — it never rebuilt it. + @test PNM.get_core(vlodf).K === core.K + @test PNM.get_core(vmodf).K === core.K +end From 04ba5351f36a69539cd85f0edd94b323aa91ae13 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 13:22:12 +0000 Subject: [PATCH 7/7] Address Copilot review: out-of-place solve, MODF docstring, lock scope - Generic _solve_multi_rhs! fallback now uses the value returned by _solve_factorization instead of assuming an in-place solve, so a backend that returns a fresh vector is handled correctly. - VirtualMODF populate_cache docstring no longer claims branch-name strings are valid monitored identifiers; _resolve_monitored_index accepts integer indices and bus-pair tuples (matching the MODF getindex API). - Resolve contingencies/monitored arcs inside core.solver_lock so the contingency_cache reads (UUID/Outage inputs) match getindex's locking and cannot race with clear_all_caches!. --- src/populate_cache.jl | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/populate_cache.jl b/src/populate_cache.jl index 1df8ae3d8..409f5cbee 100644 --- a/src/populate_cache.jl +++ b/src/populate_cache.jl @@ -52,8 +52,10 @@ function _solve_multi_rhs!( for p in SparseArrays.nzrange(B, j) col[SparseArrays.rowvals(B)[p]] = SparseArrays.nonzeros(B)[p] end - _solve_factorization(K, col) - copyto!(view(out, :, j), col) + # Capture the return: KLU/Accelerate solve in place and return `col`, + # but a future backend may return a fresh vector, so read the result. + result = _solve_factorization(K, col) + copyto!(view(out, :, j), result) end return out end @@ -360,7 +362,7 @@ batched multi-RHS solves. `contingencies` may mix `NetworkModification`, `ContingencySpec`, registered `PSY.Outage`, and registered outage UUIDs. `monitored` is an iterable of arc -identifiers (integer indices, bus-pair tuples, or branch names). +identifiers (integer indices or bus-pair tuples). The acceleration comes from a single batched solve over the union of all arcs referenced — every contingency's modified arcs plus every monitored arc. Because @@ -377,11 +379,6 @@ pinned, so subsequent `vmodf[monitored, contingency]` queries are cache hits. $(TYPEDSIGNATURES) """ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) - mods = - unique(NetworkModification[_resolve_modification(vmodf, c) for c in contingencies]) - mon_idx = unique(Int[_resolve_monitored_index(vmodf, m) for m in monitored]) - (isempty(mods) || isempty(mon_idx)) && return nothing - core = get_core(vmodf) row_caches = get_row_caches(vmodf) woodbury_cache = get_woodbury_cache(vmodf) @@ -392,6 +389,15 @@ function populate_cache(vmodf::VirtualMODF, contingencies; monitored) arc_sus = core.arc_susceptances @lock core.solver_lock begin + # Resolve under the lock: `_resolve_modification` reads `contingency_cache` + # for UUID/Outage inputs, so resolving here (not before the lock) matches + # `getindex` and avoids racing with `clear_all_caches!`. + mods = unique( + NetworkModification[_resolve_modification(vmodf, c) for c in contingencies], + ) + mon_idx = unique(Int[_resolve_monitored_index(vmodf, m) for m in monitored]) + (isempty(mods) || isempty(mon_idx)) && return nothing + # Union of all arcs needing a pre-contingency solve: monitored arcs and # every contingency's modified arcs. arc_set = Set{Int}(mon_idx)