make ptdf diag lazy

docs/src/how_to_guides/compute_network_matrices.md

@@ -79,6 +79,21 @@ vlodf_matrix = PNM.VirtualLODF(sys)
 vlodf_matrix[(1, 2), (3, 4)]
 ```
 
+## Computing Virtual MODF Matrix
+
+For post-contingency / post-modification PTDF rows using [`VirtualMODF`](@ref):
+
+```julia
+vmodf_matrix = PNM.VirtualMODF(sys)
+
+# Inspect contingencies auto-registered from PSY.Outage supplemental attributes
+PNM.get_registered_contingencies(vmodf_matrix)
+
+# Access the post-modification PTDF row for monitored arc (1, 2) under a contingency
+contingency = first(values(PNM.get_registered_contingencies(vmodf_matrix)))
+vmodf_matrix[(1, 2), contingency]
+```
+
 ## Computing Incidence and BA Matrices
 
 For the fundamental network topology matrices:
@@ -134,6 +149,7 @@ For matrices involving branches (IncidenceMatrix, PTDF, LODF), branches are repr
 | `Ybus`            | Bus numbers        | Bus numbers           |
 | `VirtualPTDF`     | Arc tuples         | Bus numbers           |
 | `VirtualLODF`     | Arc tuples         | Arc tuples            |
+| `VirtualMODF`     | Arc tuples         | Bus numbers           |
 
 !!! note
 

docs/src/reference/network_matrices_overview.md

@@ -93,6 +93,19 @@ LODF is derived from PTDF using matrix operations that model the effect of remov
 
 See [`VirtualLODF`](@ref) for cases where it is not possible to compute or store the full LODF matrix.
 
+### Modification Distribution Factors ([`VirtualMODF`](@ref))
+
+The MODF answers: "Under a registered contingency or network modification (one or more branch outages, partial-susceptance changes, or shunt changes), what is the resulting PTDF row for a given monitored arc?"
+
+**Key Characteristics:**
+
+  - On-demand computation of post-modification PTDF rows via the Woodbury matrix identity over a base PTDF.
+  - Caches Woodbury factors per modification and post-modification PTDF rows per `(monitored_arc, modification)` pair.
+  - Auto-registers `PSY.Outage` supplemental attributes attached to the source system.
+  - Rows are indexed by monitored arc tuples or arc indices; the second index is a [`ContingencySpec`](@ref), a [`NetworkModification`](@ref), or a `PSY.Outage`.
+
+**Use it when:** you need post-contingency or post-modification flow sensitivities for one or a few contingencies on a large system, and constructing the full LODF or recomputing PTDF per contingency is too expensive.
+
 ## Arc-Based Indexing
 
 All matrices that involve branches use **arc tuples** as identifiers instead of branch name strings. An arc tuple is a `Tuple{Int, Int}` of the form `(from_bus_number, to_bus_number)`, representing the directed connection between two buses.
@@ -119,6 +132,7 @@ Each matrix stores `axes` and `lookup` fields:
 | `Ybus`            | Bus numbers        | Bus numbers           |
 | `VirtualPTDF`     | Arc tuples         | Bus numbers           |
 | `VirtualLODF`     | Arc tuples         | Arc tuples            |
+| `VirtualMODF`     | Arc tuples         | Bus numbers           |
 
 !!! note
 

src/virtual_lodf_calculations.jl

@@ -17,8 +17,9 @@ every libklu solve runs under `_LIBKLU_LOCK` (process-wide) and the per-cache
 callers can issue requests concurrently; the libklu work runs one at a time.
 
 # Arguments
-- `K::KLULinSolveCache{Float64}`:
-        ABA factorization (single cache; solves serialize via the locks above).
+- `K`:
+        ABA factorization cache (`KLULinSolveCache{Float64}` or `AAFactorization`).
+        Solves serialize via the locks above.
 - `BA::SparseArrays.SparseMatrixCSC{Float64, Int}`:
         BA matrix.
 - `A::SparseArrays.SparseMatrixCSC{Int8, Int}`:
@@ -109,8 +110,24 @@ function Base.show(io::IO, ::MIME{Symbol("text/plain")}, array::VirtualLODF)
     return
 end
 
+# Map bus index (1:n_buses) to its position in `valid_ix`; ref buses map to 0.
+function _build_bus_to_valid_idx(n_buses::Int, valid_ix::Vector{Int})::Vector{Int}
+    bus_to_valid = zeros(Int, n_buses)
+    for (vi, bus_ix) in enumerate(valid_ix)
+        bus_to_valid[bus_ix] = vi
+    end
+    return bus_to_valid
+end
+
+"""
+    _get_PTDF_A_diag(K, BA, A, ref_bus_positions) -> Vector{Float64}
+
+Compute `diag(PTDF · A)`. Each row of `A` has exactly two nonzeros (+1 at the
+from-bus, -1 at the to-bus), so the per-arc dot product reduces to two indexed
+reads into the solved PTDF row after a one-time transpose of `A`.
+"""
 function _get_PTDF_A_diag(
-    K::KLULinSolveCache{Float64},
+    K,
     BA::SparseArrays.SparseMatrixCSC{Float64, Int},
     A::SparseArrays.SparseMatrixCSC{Int8, Int},
     ref_bus_positions::Set{Int},
@@ -119,35 +136,61 @@ function _get_PTDF_A_diag(
     n_buses = size(BA, 1)
     diag_ = zeros(n_branches)
 
-    # Pre-compute valid indices (non-reference buses)
     valid_ix = setdiff(1:n_buses, ref_bus_positions)
     n_valid = length(valid_ix)
+    bus_to_valid_idx = _build_bus_to_valid_idx(n_buses, valid_ix)
+
+    # Per-arc (from_valid, to_valid) via one transpose of A; 0 = ref bus.
+    A_T = SparseArrays.sparse(transpose(A))
+    arc_from_valid = Vector{Int}(undef, n_branches)
+    arc_to_valid = Vector{Int}(undef, n_branches)
+    at_rv = SparseArrays.rowvals(A_T)
+    at_nz = SparseArrays.nonzeros(A_T)
+    for i in 1:n_branches
+        f_valid = 0
+        t_valid = 0
+        @inbounds for k in SparseArrays.nzrange(A_T, i)
+            bus_ix = at_rv[k]
+            v = at_nz[k]
+            valid_i = bus_to_valid_idx[bus_ix]
+            if v > 0
+                f_valid = valid_i
+            elseif v < 0
+                t_valid = valid_i
+            end
+        end
+        arc_from_valid[i] = f_valid
+        arc_to_valid[i] = t_valid
+    end
 
-    # Pre-allocate work arrays for efficiency
     ba_col = zeros(n_valid)
-    ptdf_row = zeros(n_buses)
+    ba_rv = SparseArrays.rowvals(BA)
+    ba_nz = SparseArrays.nonzeros(BA)
 
-    # For each branch, compute PTDF row and dot with incidence column
     for i in 1:n_branches
-        # Extract BA column for valid indices
         fill!(ba_col, 0.0)
-        for idx in 1:n_valid
-            bus_idx = valid_ix[idx]
-            ba_col[idx] = BA[bus_idx, i]
+        @inbounds for k in SparseArrays.nzrange(BA, i)
+            valid_i = bus_to_valid_idx[ba_rv[k]]
+            valid_i > 0 || continue
+            ba_col[valid_i] = ba_nz[k]
         end
 
-        solve!(K, ba_col)
+        _solve_factorization(K, ba_col)
 
-        # Map back to full bus indices
-        fill!(ptdf_row, 0.0)
-        for idx in 1:n_valid
-            ptdf_row[valid_ix[idx]] = ba_col[idx]
+        # ba_col is now PTDF row i in valid-index space; H[e,e] = ptdf[from] - ptdf[to].
+        f = arc_from_valid[i]
+        t = arc_to_valid[i]
+        v_f = if f > 0
+            ba_col[f]
+        else
+            0.0
         end
-
-        # Compute diagonal element: sum of PTDF[i,j] * A[i,j] for all buses j
-        for j in 1:n_buses
-            diag_[i] += ptdf_row[j] * A[i, j]
+        v_t = if t > 0
+            ba_col[t]
+        else
+            0.0
         end
+        @inbounds diag_[i] = v_f - v_t
     end
     return diag_
 end

src/virtual_modf_calculations.jl

@@ -30,7 +30,8 @@ cache and skips the recomputation.
 - `A::SparseArrays.SparseMatrixCSC{Int8, Int}`:
         Incidence matrix.
 - `PTDF_A_diag::Vector{Float64}`:
-        Raw diagonal elements of PTDF·A (H[e,e] values).
+        Diagonal of `PTDF·A` (H[e,e] values). Lazily populated on the first
+        read of `vmodf.PTDF_A_diag`; empty until then.
 - `arc_susceptances::Vector{Float64}`:
         Effective susceptance for each arc.
 - `branch_susceptances_by_arc::Vector{Vector{Float64}}`:
@@ -111,6 +112,49 @@ get_bus_axis(mat::VirtualMODF) = mat.axes[2]
 get_tol(mat::VirtualMODF) = mat.tol[]
 get_system_uuid(M::VirtualMODF) = M.system_uuid
 
+function Base.getproperty(vmodf::VirtualMODF, name::Symbol)
+    name === :PTDF_A_diag && return _get_ptdf_a_diag_lazy!(vmodf)
+    return getfield(vmodf, name)
+end
+
+# Use `getfield` for every field read here — `vmodf.PTDF_A_diag` would
+# recurse via `getproperty`.
+function _get_ptdf_a_diag_lazy!(vmodf::VirtualMODF)
+    diag = getfield(vmodf, :PTDF_A_diag)
+    !isempty(diag) && return diag
+    @lock getfield(vmodf, :solver_lock) begin
+        diag = getfield(vmodf, :PTDF_A_diag)
+        !isempty(diag) && return diag
+        n_arcs = length(getfield(vmodf, :axes)[1])
+        @info "Computing VirtualMODF.PTDF_A_diag on first access ($n_arcs arcs)."
+        t0 = time_ns()
+        K = getfield(vmodf, :K)
+        BA = getfield(vmodf, :BA)
+        A = getfield(vmodf, :A)
+        ref_bus_set = _ref_bus_positions(vmodf)
+        new_diag = _get_PTDF_A_diag(K, BA, A, ref_bus_set)
+        resize!(diag, length(new_diag))
+        copyto!(diag, new_diag)
+        elapsed = (time_ns() - t0) / 1e9
+        @info "Computed VirtualMODF.PTDF_A_diag in $(round(elapsed; digits = 2)) s (cached)."
+        return diag
+    end
+end
+
+function _ref_bus_positions(vmodf::VirtualMODF)
+    n_buses = length(getfield(vmodf, :axes)[2])
+    valid_ix = getfield(vmodf, :valid_ix)
+    return Set{Int}(setdiff(1:n_buses, valid_ix))
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Return `H[e, e]` for each arc `e`. Same as `vmodf.PTDF_A_diag`; both trigger
+the lazy compute on first call and return the cached vector thereafter.
+"""
+get_PTDF_A_diag(vmodf::VirtualMODF) = vmodf.PTDF_A_diag
+
 # Woodbury kernel accessors. The Woodbury kernel always goes through
 # `with_solver` so it picks up the `solver_lock` + `_LIBKLU_LOCK` chain.
 _get_BA(m::VirtualMODF) = m.BA
@@ -242,8 +286,8 @@ function VirtualMODF(
 
     valid_ix = setdiff(1:length(bus_ax), ref_bus_positions)
 
-    # PTDF diagonal precompute runs serially on the dispatcher thread.
-    PTDF_A_diag = _get_PTDF_A_diag(K, BA.data, A.data, Set(ref_bus_positions))
+    # Empty: populated lazily on first read of `vmodf.PTDF_A_diag`.
+    PTDF_A_diag = Float64[]
     arc_susceptances = _extract_arc_susceptances(BA.data)
     branch_susceptances_by_arc =
         _extract_branch_susceptances_by_arc(BA.data, arc_ax, Ymatrix.network_reduction_data)

test/test_modf_lodf_reductions.jl

@@ -62,6 +62,90 @@ function verify_modf_lodf_identity(
     return
 end
 
+"""
+Verify the N-2 simultaneous-outage identity for VirtualMODF.
+
+For full outage of arcs `e1` and `e2` simultaneously, checks:
+
+    VirtualMODF[m, ctg] ≈ PTDF[m,:] + eff1*PTDF[e1,:] + eff2*PTDF[e2,:]
+
+for every monitored arc m, where eff1/eff2 come from the 2×2 LODF
+correction (Woodbury identity on the pre-contingency LODF matrix):
+
+    denom = 1 - LODF[e1,e2] * LODF[e2,e1]
+    eff1  = (LODF[m,e1] + LODF[e2,e1] * LODF[m,e2]) / denom
+    eff2  = (LODF[e1,e2] * LODF[m,e1] + LODF[m,e2]) / denom
+
+Raises an error if the pair {e1, e2} is an N-2 islanding pair
+(denom ≈ 0), since the closed-form formula is undefined in that case.
+"""
+function verify_modf_n2_lodf_identity(
+    vmodf::VirtualMODF,
+    vlodf::VirtualLODF,
+    ptdf::PTDF,
+    arc_idx1::Int,
+    arc_idx2::Int;
+    atol = 1e-6,
+)
+    n_arcs = size(vlodf, 1)
+    b_arc1 = vmodf.arc_susceptances[arc_idx1]
+    b_arc2 = vmodf.arc_susceptances[arc_idx2]
+
+    L12 = vlodf[arc_idx1, arc_idx2]
+    L21 = vlodf[arc_idx2, arc_idx1]
+    denom = 1 - L12 * L21
+    abs(denom) < 1e-6 &&
+        error(
+            "Arc pair ($arc_idx1, $arc_idx2) is an N-2 islanding pair (denom=$denom); use a non-islanding pair",
+        )
+
+    ctg_uuid = Base.UUID(UInt128(hash((arc_idx1, arc_idx2, :n2))))
+    ctg = ContingencySpec(
+        ctg_uuid,
+        NetworkModification(
+            "n2_test_$(arc_idx1)_$(arc_idx2)",
+            [ArcModification(arc_idx1, -b_arc1), ArcModification(arc_idx2, -b_arc2)],
+        ),
+    )
+    vmodf.contingency_cache[ctg_uuid] = ctg
+
+    for m in 1:n_arcs
+        Lm1 = vlodf[m, arc_idx1]
+        Lm2 = vlodf[m, arc_idx2]
+        eff1 = (Lm1 + L21 * Lm2) / denom
+        eff2 = (L12 * Lm1 + Lm2) / denom
+
+        modf_row = vmodf[m, ctg]
+        expected = ptdf[m, :] .+ eff1 .* ptdf[arc_idx1, :] .+ eff2 .* ptdf[arc_idx2, :]
+        @test isapprox(modf_row, expected; atol = atol)
+    end
+
+    empty!(vmodf.woodbury_cache)
+    return
+end
+
+"""
+Return two distinct non-islanding arcs from `branch_map` whose simultaneous
+outage does not island the network (L12*L21 ≠ 1 in the pre-contingency LODF).
+Raises an error if no such pair exists in the map.
+"""
+function _find_two_non_islanding_arcs(vmodf, vlodf, branch_map)
+    candidates = Tuple{Any, Int}[]
+    for (arc_tuple, _) in branch_map
+        arc_idx = vmodf.lookup[1][arc_tuple]
+        if abs(vmodf.PTDF_A_diag[arc_idx]) < 1.0 - 1e-6
+            push!(candidates, (arc_tuple, arc_idx))
+        end
+    end
+    for i in eachindex(candidates), j in (i + 1):length(candidates)
+        _, idx1 = candidates[i]
+        _, idx2 = candidates[j]
+        denom = 1 - vlodf[idx1, idx2] * vlodf[idx2, idx1]
+        abs(denom) >= 1e-6 && return candidates[i], candidates[j]
+    end
+    error("No non-islanding arc pair found in map ($(length(candidates)) candidates)")
+end
+
 """
 Find a representative non-islanding arc from a branch map.
 
@@ -185,3 +269,16 @@ end
         verify_modf_lodf_identity(vmodf, vlodf, ptdf, arc_idx, delta_b)
     end
 end
+
+@testset "MODF vs LODF N-2: direct branch pair" begin
+    sys = PSB.build_system(PSSEParsingTestSystems, "psse_14_network_reduction_test_system")
+    reductions = NetworkReduction[]
+    ptdf = PTDF(sys; network_reductions = reductions)
+    vlodf = VirtualLODF(sys; network_reductions = reductions)
+    vmodf = VirtualMODF(sys; network_reductions = reductions)
+    nrd = get_network_reduction_data(vmodf)
+
+    (_, arc_idx1), (_, arc_idx2) =
+        _find_two_non_islanding_arcs(vmodf, vlodf, nrd.direct_branch_map)
+    verify_modf_n2_lodf_identity(vmodf, vlodf, ptdf, arc_idx1, arc_idx2)
+end

test/test_virtual_lodf.jl

@@ -147,3 +147,45 @@ end
     end
     @test test_value
 end
+
+@testset "_get_PTDF_A_diag: matches reference implementation" begin
+    # Reference: full-bus dot product. Slow, unambiguously correct.
+    function _reference_ptdf_a_diag(K, BA, A, ref_bus_positions::Set{Int})
+        n_branches = size(BA, 2)
+        n_buses = size(BA, 1)
+        valid_ix = setdiff(1:n_buses, ref_bus_positions)
+        n_valid = length(valid_ix)
+        diag_ = zeros(n_branches)
+        ba_col = zeros(n_valid)
+        ptdf_row = zeros(n_buses)
+        for i in 1:n_branches
+            fill!(ba_col, 0.0)
+            for idx in 1:n_valid
+                bus_idx = valid_ix[idx]
+                ba_col[idx] = BA[bus_idx, i]
+            end
+            PNM._solve_factorization(K, ba_col)
+            fill!(ptdf_row, 0.0)
+            for idx in 1:n_valid
+                ptdf_row[valid_ix[idx]] = ba_col[idx]
+            end
+            for j in 1:n_buses
+                diag_[i] += ptdf_row[j] * A[i, j]
+            end
+        end
+        return diag_
+    end
+
+    for case in ("c_sys5", "c_sys14")
+        sys = PSB.build_system(PSB.PSITestSystems, case)
+        A = PNM.IncidenceMatrix(sys)
+        BA = PNM.BA_Matrix(sys)
+        ref_pos = Set(PNM.get_ref_bus_position(A))
+        ABA = PNM.calculate_ABA_matrix(A.data, BA.data, ref_pos)
+        K = PNM.klu_factorize(ABA)
+
+        fast = PNM._get_PTDF_A_diag(K, BA.data, A.data, ref_pos)
+        ref = _reference_ptdf_a_diag(K, BA.data, A.data, ref_pos)
+        @test fast ≈ ref atol = 1e-12 rtol = 0
+    end
+end