Merge pull request #48 from NREL-Sienna/add_LODF_functions

Add lodf functions

src/BA_ABA_matrices.jl

@@ -3,16 +3,16 @@ Structure containing the BA matrix and other relevant data.
 
 # Arguments
 - `data::SparseArrays.SparseMatrixCSC{Float64, Int}`:
-        the BA matrix coming from the product between the Incidence Matrix A and
-        the Matrix of Susceptance B
+        the transposed BA matrix coming from the product between the Incidence 
+        Matrix A and the Matrix of Susceptance B
 - `axes<:NTuple{2, Dict}`:
         Tuple containing two vectors, the first one contains the names of each
-        line of the network (each one related to a row of the Matrix in "data"),
-        the second one contains the names of each bus of the network (each one
+        buse of the network (each one related to a row of the Matrix in "data"),
+        the second one contains the names of each line of the network (each one
         related to a column of the Matrix in "data")
 - `lookup<:NTuple{2, Dict}`:
         Tuple containing 2 Dictionaries mapping the number of rows and columns
-        with the names of branches and buses
+        with the names of buses and branches
 - `ref_bus_positions::Set{Int}`:
         Vector containing the indexes of the columns of the BA matrix corresponding
         to the refence buses
@@ -34,8 +34,8 @@ function BA_Matrix(sys::PSY.System)
     bus_lookup = make_ax_ref(buses)
     line_ax = [PSY.get_name(branch) for branch in branches]
     bus_ax = [PSY.get_number(bus) for bus in setdiff(buses, ref_bus_positions)]
-    axes = (line_ax, bus_ax)
-    lookup = (make_ax_ref(line_ax), make_ax_ref(bus_ax))
+    axes = (bus_ax, line_ax)
+    lookup = (make_ax_ref(bus_ax), make_ax_ref(line_ax))
     data = calculate_BA_matrix(branches, bus_lookup)
     return BA_Matrix(data, axes, lookup, ref_bus_positions)
 end
@@ -46,15 +46,14 @@ Structure containing the ABA matrix and other relevant data.
 # Arguments
 - `data::SparseArrays.SparseMatrixCSC{Float64, Int}`:
         the ABA matrix coming from the product between the Incidence Matrix A and
-        the Matrix BA
+        the Matrix BA.
 - `axes<:NTuple{2, Dict}`:
-        Tuple containing two vectors, the first one contains the names of each
-        line of the network (each one related to a row of the Matrix in "data"),
-        the second one contains the names of each bus of the network (each one
-        related to a column of the Matrix in "data")
+        Tuple containing two identical vectors, both containing the number of 
+        each bus of the network (each one related to a row/column of the Matrix
+        in "data"), excluding the slack buses.
 - `lookup<:NTuple{2, Dict}`:
         Tuple containing 2 Dictionaries mapping the number of rows and columns
-        with the names of branches and buses
+        with the number of the buses.
 - `ref_bus_positions::Set{Int}`:
         Vector containing the indexes of the columns of the BA matrix corresponding
         to the refence buses
@@ -95,7 +94,7 @@ function ABA_Matrix(sys::PSY.System; factorize = false)
     bus_ax = [PSY.get_number(bus) for bus in buses]
     bus_ax_ = setdiff(bus_ax, ref_bus_positions)
     axes = (bus_ax_, bus_ax_)
-    bus_ax_ref = make_ax_ref(bus_ax)
+    bus_ax_ref = make_ax_ref(bus_ax_)
     lookup = (bus_ax_ref, bus_ax_ref)
     if factorize
         K = klu(ABA)
@@ -134,3 +133,36 @@ is_factorized(ABA::ABA_Matrix{Ax, L, Nothing}) where {Ax, L <: NTuple{2, Dict}}
 is_factorized(
     ABA::ABA_Matrix{Ax, L, KLU.KLUFactorization{Float64, Int}},
 ) where {Ax, L <: NTuple{2, Dict}} = true
+
+# get_index functions: BA_Matrix stores the transposed matrix, thus get index
+# must export values according to [branch, bus] indexing.
+function Base.getindex(A::BA_Matrix, line, bus)
+    i, j = to_index(A, bus, line)
+    return A.data[i, j]
+end
+
+function Base.getindex(
+    A::BA_Matrix,
+    line_number::Union{Int, Colon},
+    bus_number::Union{Int, Colon},
+)
+    return A.data[bus_number, line_number]
+end
+
+# get_index functions: ABA_Matrix stores a square matrix whose number of rows 
+# and column is equal to the number of the system's buses minus the slack ones,
+# NOTE: bus_1, bus_2 are bus numbers not row and column indices!
+function Base.getindex(A::ABA_Matrix, bus_1, bus_2)
+    if bus_1 in A.ref_bus_positions || bus_2 in A.ref_bus_positions
+        err_msg = string(
+            " Rows and columns related to slack buses are not defined for the ABA matrix. \n",
+            "Indices must be referred to any bus number that is not a slack one. \n",
+            "For the current Systems the reference slack buses are: ",
+            collect(A.ref_bus_positions), ". \n")
+        error(err_msg)
+    else
+        i, j = to_index(A, bus_1, bus_2)
+        return A.data[i, j]
+    end
+    return
+end

src/PowerNetworkMatrices.jl

@@ -3,11 +3,11 @@ module PowerNetworkMatrices
 export ABA_Matrix
 export AdjacencyMatrix
 export BA_Matrix
-export drop_small_entries!
 export factorize
 export find_subnetworks
 export from_hdf5
 export get_ptdf_data
+export get_lodf_data
 export IncidenceMatrix
 export is_factorized
 export LODF

src/PowerNetworkMatrix.jl

@@ -224,11 +224,11 @@ function Base.summary(io::IO, A::PowerNetworkMatrix)
         show(IOContext(io, :limit => true), ax)
         println(io)
     end
-    print(io, "And data, a ", size(A))
+    print(io, "And data with size ", size(A))
     return
 end
 
-_summary(io::IO, A::PowerNetworkMatrix) = println(io, "PowerNetworkMatrix")
+_summary(io::IO, ::T) where {T <: PowerNetworkMatrix} = println(io, "$T")
 
 function Base.summary(
     io::IOContext{Base.GenericIOBuffer{Array{UInt8, 1}}},
@@ -303,7 +303,7 @@ function Base.show_nd(
     end
 end
 
-function Base.show(io::IO, array::PowerNetworkMatrix)
+function Base.show(io::IO, ::MIME{Symbol("text/plain")}, array::PowerNetworkMatrix)
     summary(io, array)
     isempty(array) && return
     println(io, ":")

src/common.jl

@@ -226,7 +226,7 @@ value is above a certain tolerance.
 function sparsify(dense_array::Matrix{Float64}, tol::Float64)
     m, n = size(dense_array)
     sparse_array = SparseArrays.spzeros(m, n)
-    for i in 1:m, j in 1:n
+    for j in 1:n, i in 1:m
         if abs(dense_array[i, j]) > tol
             sparse_array[i, j] = dense_array[i, j]
         end
@@ -235,57 +235,25 @@ function sparsify(dense_array::Matrix{Float64}, tol::Float64)
 end
 
 """
-Sets to zero every element of a Sparse matrix if absolute values is below a
-certain tolerance.
-
-# Arguments
-- `sparse_array::SparseArrays.SparseMatrixCSC{Float64, Int}`: input sparse array.
-- `tol::Float64`: tolerance for removing entries in the PTDF matrix.
-"""
-function make_entries_zero!(
-    sparse_array::SparseArrays.SparseMatrixCSC{Float64, Int},
-    tol::Float64,
-)
-    for i in 1:size(sparse_array, 1)
-        sparse_array[i, abs.(sparse_array[i, :]) .<= tol] .= 0.0
-    end
-    SparseArrays.dropzeros!(sparse_array)
-    return
-end
-
-"""
-Sets to zero every element of a Dense matrix if absolute values is below a
-certain tolerance.
-
-# Arguments
-- `dense_array::Matrix{Float64}`: input dense matrix.
-- `tol::Float64`: tolerance.
-"""
-function make_entries_zero!(
-    dense_array::Matrix{Float64},
-    tol::Float64,
-)
-    for i in 1:size(dense_array, 1)
-        dense_array[i, abs.(dense_array[i, :]) .<= tol] .= 0.0
-    end
-    return
-end
+Return a sparse vector given a dense one by dropping element whose absolute
+value is above a certain tolerance.
 
-"""
-Sets to zero every element of a Dense vector if absolute values is below a
-certain tolerance.
 
 # Arguments
-- `vector::Vector{Float64}`:input dense vector.
-- `tol::Float64`: tolerance.
+- dense_array::Vector{Float64}`:
+        input vector (e.g., PTDF row from VirtualPTDF).
+- `tol::Float64`:
+        tolerance.
 """
-function make_entries_zero!(vector::Vector{Float64}, tol::Float64)
-    for i in eachindex(vector)
-        if abs(vector[i]) <= tol
-            vector[i] = 0.0
+function sparsify(dense_array::Vector{Float64}, tol::Float64)
+    m = length(dense_array)
+    sparse_array = SparseArrays.spzeros(m)
+    for i in 1:m
+        if abs(dense_array[i]) > tol
+            sparse_array[i] = dense_array[i]
         end
     end
-    return
+    return sparse_array
 end
 
 """