Merge pull request #296 from Sienna-Platform/claude/klu-sparse-rhs-wrapper-jLAmX

Add KLUWrapper

Author: jd-lara

.claude/Sienna.md

@@ -72,6 +72,10 @@ Key rules:
   - Globals: `UPPER_CASE` for constants
   - Exports: all exports in main module file
   - Comments: complete sentences, describe why not how
+  - Nothing checks: use `isnothing(x)` / `!isnothing(x)`. Do **not** use the older `x === nothing` / `x !== nothing` patterns.
+  - Type checks: prefer multiple dispatch over `isa` (the only acceptable use of `isa` is filtering inside a `catch` block, where dispatch is unavailable).
+  - Conditionals: prefer `if/else` over the ternary `? :` operator for readability, especially in multi-line expressions.
+  - Cache lookups: `get!(dict, key) do ... end` (the closure form is **lazy** — only evaluates on a miss). Never use the 3-arg `get!(dict, key, default)` when `default` is expensive: Julia evaluates function arguments eagerly, so `default` runs on every call (including cache hits) and silently defeats the cache.
 
 ## Documentation Practices and Requirements
 

.github/workflows/main-tests.yml

@@ -31,6 +31,7 @@ jobs:
         continue-on-error: ${{ matrix.julia-version == 'nightly' }}
         env:
           PYTHON: ""
+          JULIA_NUM_THREADS: "2"
       - uses: julia-actions/julia-processcoverage@v1
         if: matrix.julia-version != 'nightly'
         with:

.github/workflows/pr_testing.yml

@@ -26,6 +26,7 @@ jobs:
       - uses: julia-actions/julia-runtest@latest
         env:
           PYTHON: ""
+          JULIA_NUM_THREADS: "2"
       - uses: julia-actions/julia-processcoverage@v1
         with:
           directories: src,ext

Project.toml

@@ -1,18 +1,18 @@
 name = "PowerNetworkMatrices"
 uuid = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
-version = "0.20.0"
+version = "0.21.0"
 authors = ["SIENNA TEAM"]
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
-KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SuiteSparse_jll = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
 
 [weakdeps]
 AppleAccelerate = "13e28ba4-7ad8-5781-acae-3021b1ed3924"
@@ -28,10 +28,10 @@ DataStructures = "^0.19"
 DocStringExtensions = "~0.8, ~0.9"
 HDF5 = "0.17"
 InfrastructureSystems = "3"
-KLU = "^0.6"
 LinearAlgebra = "1"
 Pardiso = "1"
-PowerSystems = "^5.7"
+PowerSystems = "^5.10"
 Preferences = "^1.5"
 SparseArrays = "1"
+SuiteSparse_jll = "5, 6, 7"
 julia = "^1.10"

docs/make.jl

@@ -36,6 +36,7 @@ pages = OrderedDict(
     "Reference" => Any[
         "Matrix Overview" => "reference/network_matrices_overview.md",
         "Public API" => "reference/public.md",
+        "Internals" => "reference/internals.md",
     ],
 )
 

docs/src/reference/internals.md

@@ -0,0 +1,16 @@
+# Internals
+
+The symbols documented on this page are **internal** to `PowerNetworkMatrices` and are
+not part of the public API. They are documented here so that the published manual covers
+every docstring shipped with the package, but they may change at any time without notice
+and should not be relied on by downstream packages.
+
+## `KLUWrapper`
+
+`PowerNetworkMatrices.KLUWrapper` is a thin, allocation-aware wrapper over `libklu`
+(provided by `SuiteSparse_jll`) used internally for sparse linear solves. None of these
+symbols are exported from `PowerNetworkMatrices`.
+
+```@autodocs
+Modules = [PowerNetworkMatrices.KLUWrapper]
+```

src/BA_ABA_matrices.jl

@@ -159,7 +159,7 @@ power flow analysis, sensitivity calculations, and linear power system studies.
         Mapping from reference bus numbers to their corresponding subnetwork axes
 - `ref_bus_position::Vector{Int}`:
         Vector containing the original indices of reference buses before matrix reduction
-- `K::F <: Union{Nothing, KLU.KLUFactorization{Float64, Int}}`:
+- `K::F <: Union{Nothing, KLULinSolveCache{Float64}}`:
         Optional KLU factorization object for efficient linear system solving. Nothing if unfactorized
 - `network_reduction_data::NetworkReductionData`:
         Container for network reduction information applied during matrix construction
@@ -179,7 +179,7 @@ power flow analysis, sensitivity calculations, and linear power system studies.
 struct ABA_Matrix{
     Ax <: NTuple{2, Vector},
     L <: NTuple{2, Dict},
-    F <: Union{Nothing, KLU.KLUFactorization{Float64, Int}},
+    F <: Union{Nothing, KLULinSolveCache{Float64}},
 } <: PowerNetworkMatrix{Float64}
     data::SparseArrays.SparseMatrixCSC{Float64, Int}
     axes::Ax
@@ -295,7 +295,7 @@ function ABA_Matrix(ybus::Ybus; factorize::Bool = false)
     bus_ax_ref = make_ax_ref(axes[1])
     lookup = (bus_ax_ref, bus_ax_ref)
     if factorize
-        K = klu(ABA)
+        K = klu_factorize(ABA)
     else
         K = nothing
     end
@@ -339,7 +339,7 @@ function factorize(ABA::ABA_Matrix{Ax, L, Nothing}) where {Ax, L <: NTuple{2, Di
         deepcopy(ABA.lookup),
         deepcopy(ABA.subnetwork_axes),
         deepcopy(ABA.ref_bus_position),
-        klu(ABA.data),
+        klu_factorize(ABA.data),
         deepcopy(ABA.network_reduction_data),
     )
     return ABA_lu
@@ -358,7 +358,7 @@ Check if an ABA_Matrix has been factorized (i.e., contains LU factorization matr
 """
 is_factorized(ABA::ABA_Matrix{Ax, L, Nothing}) where {Ax, L <: NTuple{2, Dict}} = false
 is_factorized(
-    ABA::ABA_Matrix{Ax, L, KLU.KLUFactorization{Float64, Int}},
+    ABA::ABA_Matrix{Ax, L, <:KLULinSolveCache{Float64}},
 ) where {Ax, L <: NTuple{2, Dict}} = true
 
 # get_index functions: BA_Matrix stores the transposed matrix, thus get index

src/BranchesParallel.jl

@@ -1,23 +1,45 @@
-mutable struct BranchesParallel{T <: PSY.ACTransmission} <: PSY.ACTransmission
+abstract type AbstractBranchesParallel <: PSY.ACTransmission end
+
+mutable struct BranchesParallel{T <: PSY.ACTransmission} <: AbstractBranchesParallel
     branches::Vector{T}
     equivalent_ybus::Union{Matrix{YBUS_ELTYPE}, Nothing}
+
+    function BranchesParallel{T}(
+        branches::Vector{T},
+        equivalent_ybus::Union{Matrix{YBUS_ELTYPE}, Nothing},
+    ) where {T <: PSY.ACTransmission}
+        if !isconcretetype(T)
+            error(
+                "BranchesParallel{T} requires a concrete branch type T. " *
+                "Use MixedBranchesParallel for groups with mixed branch types. Got T=$T.",
+            )
+        end
+        return new{T}(branches, equivalent_ybus)
+    end
 end
 
 function BranchesParallel(branches::Vector{T}) where {T <: PSY.ACTransmission}
-    BranchesParallel(branches, nothing)
+    return BranchesParallel{T}(branches, nothing)
 end
-# Constructor for the mixed types
-function BranchesParallel(branches::Vector{PSY.ACTransmission})
-    return BranchesParallel{PSY.ACTransmission}(branches, nothing)
+
+mutable struct MixedBranchesParallel <: AbstractBranchesParallel
+    branches::Vector{PSY.ACTransmission}
+    equivalent_ybus::Union{Matrix{YBUS_ELTYPE}, Nothing}
+end
+
+function MixedBranchesParallel(branches::Vector{<:PSY.ACTransmission})
+    return MixedBranchesParallel(Vector{PSY.ACTransmission}(branches), nothing)
 end
 
 function add_branch!(bp::BranchesParallel{T}, branch::T) where {T <: PSY.ACTransmission}
     push!(bp.branches, branch)
 end
 
-get_branch_type(::BranchesParallel{T}) where {T <: PSY.ACTransmission} = T
+function add_branch!(mbp::MixedBranchesParallel, branch::PSY.ACTransmission)
+    push!(mbp.branches, branch)
+end
 
-function get_name(bp::BranchesParallel{T}) where {T <: PSY.ACTransmission}
+function get_name(bp::AbstractBranchesParallel)
     base_string = _longest_starting_substring(PSY.get_name.(bp.branches)...)
     if isempty(base_string)
         base_string = join(PSY.get_name.(bp.branches), "_") * "_"
@@ -44,7 +66,7 @@ function _longest_starting_substring(branch_names...)
 end
 
 function compute_parallel_multiplier(
-    parallel_branch_set::BranchesParallel,
+    parallel_branch_set::AbstractBranchesParallel,
     branch_name::String,
 )
     b_total = 0.0
@@ -58,45 +80,82 @@ function compute_parallel_multiplier(
     return b_branch / b_total
 end
 
-function get_series_susceptance(segment::BranchesParallel)
+function get_series_susceptance(segment::AbstractBranchesParallel)
     return sum(get_series_susceptance(branch) for branch in segment.branches)
 end
 
-function get_equivalent_physical_branch_parameters(bp::BranchesParallel)
+function get_equivalent_physical_branch_parameters(bp::AbstractBranchesParallel)
     if isnothing(bp.equivalent_ybus)
         populate_equivalent_ybus!(bp)
     end
     equivalent_ybus = bp.equivalent_ybus
     return _get_equivalent_physical_branch_parameters(equivalent_ybus)
 end
 
-function populate_equivalent_ybus!(bp::BranchesParallel)
+function populate_equivalent_ybus!(bp::AbstractBranchesParallel)
     Y11, Y12, Y21, Y22 = ybus_branch_entries(bp)
     bp.equivalent_ybus = YBUS_ELTYPE[Y11 Y12; Y21 Y22]
     return
 end
 
 """
-    get_equivalent_rating(bp::BranchesParallel)
+    get_sum_of_max_rating(bp::AbstractBranchesParallel)
 
-Calculate the total rating for branches in parallel.
-For parallel circuits, the rating is the sum of individual ratings divided by the number of circuits.
-This provides a conservative estimate that accounts for potential overestimation of total capacity.
+Sum of the individual branch ratings, treating each circuit as independently loadable
+up to its own thermal limit. This is the least conservative aggregate and assumes
+unconstrained flow steering across the parallel group.
+"""
+function get_sum_of_max_rating(bp::AbstractBranchesParallel)
+    return sum(get_equivalent_rating(branch) for branch in bp.branches)
+end
+
+"""
+    get_single_element_contingency_rating(bp::AbstractBranchesParallel)
+
+N-1 rating for the parallel group: the surviving capacity after the largest-rated
+circuit trips, ``\\sum_i S_i - \\max_i S_i``. For a group of one branch this is zero.
 """
-function get_equivalent_rating(bp::BranchesParallel)
-    # Sum of ratings divided by number of circuits
-    return sum(get_equivalent_rating(branch) for branch in bp.branches) /
-           length(bp.branches)
+function get_single_element_contingency_rating(bp::AbstractBranchesParallel)
+    ratings = get_equivalent_rating.(bp.branches)
+    return sum(ratings) - maximum(ratings)
 end
 
 """
-    get_equivalent_emergency_rating(bp::BranchesParallel)
+    get_impedance_averaged_rating(bp::AbstractBranchesParallel)
+
+Susceptance-weighted average of individual branch ratings,
+``\\sum_i f_i \\cdot S_i`` with ``f_i = b_i / \\sum_k b_k``. Reflects how DC flow
+physically splits across a parallel group. Throws `ArgumentError` if the total
+series susceptance is zero or non-finite.
+"""
+function get_impedance_averaged_rating(bp::AbstractBranchesParallel)
+    b_total = sum(PSY.get_series_susceptance, bp.branches)
+    if !isfinite(b_total) || iszero(b_total)
+        throw(
+            ArgumentError(
+                "Cannot compute impedance-averaged rating: total series susceptance across the parallel group must be finite and non-zero.",
+            ),
+        )
+    end
+    return sum(
+        PSY.get_series_susceptance(br) / b_total * get_equivalent_rating(br)
+        for br in bp.branches
+    )
+end
+
+# Series-chain rating contribution for a parallel block: dispatch arm for
+# `get_equivalent_rating(::BranchesSeries)` defined in BranchesSeries.jl.
+_series_member_rating(bp::AbstractBranchesParallel) =
+    get_single_element_contingency_rating(bp)
+
+"""
+    get_equivalent_emergency_rating(bp::AbstractBranchesParallel)
 
 Calculate the total emergency rating for branches in parallel.
 For parallel circuits, the emergency rating is the sum of individual emergency ratings divided by the number of circuits.
 This provides a conservative estimate that accounts for potential overestimation of total capacity.
 """
-function get_equivalent_emergency_rating(bp::BranchesParallel)
+function get_equivalent_emergency_rating(bp::AbstractBranchesParallel)
     equivalent_rating = 0.0
     for branch in bp.branches
         rating_b = get_equivalent_emergency_rating(branch)
@@ -106,36 +165,38 @@ function get_equivalent_emergency_rating(bp::BranchesParallel)
 end
 
 """
-    get_equivalent_available(bp::BranchesParallel)
+    get_equivalent_available(bp::AbstractBranchesParallel)
 
 Get the availability status for parallel branches.
 All branches in parallel must be available for the parallel circuit to be available.
 """
-function get_equivalent_available(bp::BranchesParallel)
+function get_equivalent_available(bp::AbstractBranchesParallel)
     # All branches must be available
     return all(PSY.get_available(branch) for branch in bp.branches)
 end
 
+PSY.get_available(bp::AbstractBranchesParallel) = get_equivalent_available(bp)
+
 """
-    get_equivalent_α(bp::BranchesParallel)
+    get_equivalent_α(bp::AbstractBranchesParallel)
 
 Get the phase angle shift for parallel branches.
 Returns the average phase angle shift across all parallel branches.
 Returns 0.0 if branches don't support phase angle shift (e.g., lines).
 """
-function get_equivalent_α(bp::BranchesParallel)
+function get_equivalent_α(bp::AbstractBranchesParallel)
     # Need to check the PS books
 end
 
-function Base.iterate(bp::BranchesParallel)
+function Base.iterate(bp::AbstractBranchesParallel)
     return iterate(bp.branches)
 end
 
-function Base.iterate(bp::BranchesParallel, state)
+function Base.iterate(bp::AbstractBranchesParallel, state)
     return iterate(bp.branches, state)
 end
 
-function Base.length(bp::BranchesParallel)
+function Base.length(bp::AbstractBranchesParallel)
     return length(bp.branches)
 end
 
@@ -144,33 +205,32 @@ function add_to_map(
     filters::Dict,
 ) where {T <: PSY.ACTransmission}
     isempty(filters) && return true
-    if isabstracttype(T)
-        @warn "Parallel circuit contains mixed branch types, filters might be applied to more components than intended. Use Logging.Debug for additional information."
-        @debug "Parallel circuit branch types: $(typeof.(double_circuit.branches))"
-        @debug "Parallel circuit branch names: $(PSY.get_name.(double_circuit.branches))"
-        for branch in double_circuit.branches
-            filter = get(filters, typeof(branch), x -> true)
-            if !filter(branch)
-                return false
-            end
-        end
+    if !haskey(filters, T)
         return true
     end
+    return any(filters[T](device) for device in double_circuit)
+end
 
-    if !haskey(filters, T)
-        return true
-    else
-        return any([filters[T](device) for device in double_circuit])
+function add_to_map(double_circuit::MixedBranchesParallel, filters::Dict)
+    isempty(filters) && return true
+    @warn "Parallel circuit contains mixed branch types, filters might be applied to more components than intended. Use Logging.Debug for additional information."
+    @debug "Parallel circuit branch types: $(typeof.(double_circuit.branches))"
+    @debug "Parallel circuit branch names: $(PSY.get_name.(double_circuit.branches))"
+    for branch in double_circuit.branches
+        filter = get(filters, typeof(branch), x -> true)
+        if !filter(branch)
+            return false
+        end
     end
-    error("Invalid condition reached in add_to_map for BranchesParallel")
+    return true
 end
 
-function Base.:(==)(a::BranchesParallel, b::BranchesParallel)
+function Base.:(==)(a::AbstractBranchesParallel, b::AbstractBranchesParallel)
     return a.branches == b.branches
 end
 
-function Base.show(io::IO, x::MIME{Symbol("text/plain")}, y::BranchesParallel)
+function Base.show(io::IO, x::MIME{Symbol("text/plain")}, y::AbstractBranchesParallel)
     show(io, x, y.branches)
 end
 
-is_a_reduction(::BranchesParallel) = true
+is_a_reduction(::AbstractBranchesParallel) = true