Merge pull request #6 from janbruedigam/adapt_to_dojo

janbruedigam · web-flow · commit 8bc68d3c9cc0 · 2023-01-21T16:31:53.000+01:00
Comply with Dojo
diff --git a/Project.toml b/Project.toml
@@ -5,6 +5,7 @@ version = "0.2.0"
 
 [deps]
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
diff --git a/README.md b/README.md
@@ -20,7 +20,7 @@ dimensions = [2; 3; 0; 1] # The dimension of row/column
 
 system = System{Float64}(graph_matrix, dimensions; symmetric=false) # The resulting linear system. Set symmetric=true for symmetric systems
 
-randomize!(system) # Randomize all system entries
+initialize!(system) # initialize! all system entries
 system.matrix_entries[1,2].value = rand(2,3) # Directly set the value of a matrix entry
 system.vector_entries[4].value = rand(1) # Directly set the value of a vector entry
 
diff --git a/benchmark/example_benchmark.jl b/benchmark/example_benchmark.jl
@@ -77,8 +77,8 @@ A[13,21] = A[21,13] = 1
 A[16,30] = A[30,16] = 1
 A[20,36] = A[36,20] = 1
 
-function randomize_posdef!(system)
-    randomize!(system,rand)
+function initialize!_posdef!(system)
+    initialize!(system,rand)
     for i=1:size(A)[1]
         system.matrix_entries[i,i].value += 1000*I 
     end
@@ -89,10 +89,10 @@ system = System{Float64}(A, ones(Int,size(A)[1])*3)
 systemldlt = System{Float64}(A, ones(Int,size(A)[1])*3, symmetric=true)
 systemllt = System{Float64}(A, ones(Int,size(A)[1])*3, symmetric=true)
 
-SUITE["ldu"] = @benchmarkable ldu_solve!($system) samples=2 setup=(randomize!($system))
-SUITE["lu"] = @benchmarkable lu_solve!($system) samples=2 setup=(randomize!($system))
-SUITE["ldlt"] = @benchmarkable ldlt_solve!($systemldlt) samples=2 setup=(randomize!($systemldlt))
-SUITE["llt"] = @benchmarkable llt_solve!($systemllt) samples=2 setup=(randomize_posdef!($systemllt))
+SUITE["ldu"] = @benchmarkable ldu_solve!($system) samples=2 setup=(initialize!($system))
+SUITE["lu"] = @benchmarkable lu_solve!($system) samples=2 setup=(initialize!($system))
+SUITE["ldlt"] = @benchmarkable ldlt_solve!($systemldlt) samples=2 setup=(initialize!($systemldlt))
+SUITE["llt"] = @benchmarkable llt_solve!($systemllt) samples=2 setup=(initialize!_posdef!($systemllt))
 
 # A = [
 #     0 1 1 1 1
diff --git a/src/GraphBasedSystems.jl b/src/GraphBasedSystems.jl
@@ -1,23 +1,35 @@
 module GraphBasedSystems
 
+using LinearAlgebra
 using SparseArrays
 using StaticArrays
 using Graphs
 
 
 export System,
+    Entry,
+
     full_matrix,
     full_vector,
-    randomize!,
+    initialize!,
+    reset_inverse_diagonals!,
     
     children,
     connections,
     parents,
 
     ldu_solve!,
+    ldu_factorization!,
+    ldu_backsubstitution!,
     lu_solve!,
+    lu_factorization!,
+    lu_backsubstitution!,
     ldlt_solve!,
-    llt_solve!
+    ldlt_factorization!,
+    ldlt_backsubstitution!,
+    llt_solve!,
+    llt_factorization!,
+    llt_backsubstitution!
 
 
 include(joinpath("util", "custom_static.jl"))
@@ -26,6 +38,9 @@ include(joinpath("system", "entry.jl"))
 include(joinpath("system", "system.jl"))
 include(joinpath("system", "setup_functions.jl"))
 
+include(joinpath("system", "interface.jl"))
+include(joinpath("system", "dense.jl"))
+
 include(joinpath("solvers", "lu.jl"))
 include(joinpath("solvers", "llt.jl"))
 include(joinpath("solvers", "ldlt.jl"))
diff --git a/src/solvers/ldlt.jl b/src/solvers/ldlt.jl
@@ -24,6 +24,8 @@ function ldlt_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             ldlt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
@@ -83,7 +85,6 @@ function ldlt_backsubstitution!(system::System)
 end
 
 function ldlt_solve!(system::System)
-    reset_inverse_diagonals!(system)
     ldlt_factorization!(system)
     ldlt_backsubstitution!(system)
     return
diff --git a/src/solvers/ldu.jl b/src/solvers/ldu.jl
@@ -24,6 +24,8 @@ function ldu_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
@@ -84,7 +86,6 @@ function ldu_backsubstitution!(system::System)
 end
 
 function ldu_solve!(system::System)
-    reset_inverse_diagonals!(system)
     ldu_factorization!(system)
     ldu_backsubstitution!(system)
     return
diff --git a/src/solvers/llt.jl b/src/solvers/llt.jl
@@ -28,6 +28,8 @@ function llt_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             llt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
@@ -92,7 +94,6 @@ function llt_backsubstitution!(system::System)
 end
 
 function llt_solve!(system::System)
-    reset_inverse_diagonals!(system)
     llt_factorization!(system)
     llt_backsubstitution!(system)
     return
diff --git a/src/solvers/lu.jl b/src/solvers/lu.jl
@@ -23,6 +23,8 @@ function lu_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             lu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
@@ -83,7 +85,6 @@ function lu_backsubstitution!(system::System)
 end
 
 function lu_solve!(system::System)
-    reset_inverse_diagonals!(system)
     lu_factorization!(system)
     lu_backsubstitution!(system)
     return
diff --git a/src/system/dense.jl b/src/system/dense.jl
@@ -0,0 +1,19 @@
+full_matrix(system::System) = full_matrix(system.matrix_entries, issymmetric(system), system.dims, system.dims)
+# There probably exists a smarter way of getting the dense matrix from the spares one 
+function full_matrix(matrix_entries::SparseMatrixCSC, symmetric::Bool, dimensions_rows, dimensions_cols)
+    range_dict_rows = ranges(dimensions_rows)
+    range_dict_cols = ranges(dimensions_cols)
+    A = zeros(sum(dimensions_rows), sum(dimensions_cols))
+
+    for (i,row) in enumerate(matrix_entries.rowval)
+        col = findfirst(x -> i<x, matrix_entries.colptr)-1
+        A[range_dict_rows[row],range_dict_cols[col]] = matrix_entries[row,col].value
+        if symmetric && col != row
+            A[range_dict_cols[col],range_dict_rows[row]] = matrix_entries[row,col].value'
+        end
+    end
+    return A
+end
+
+full_vector(system::System) = full_vector(system.vector_entries, system.dims)
+full_vector(vector_entries::AbstractVector, dimensions_rows) = vcat([getfield(vector_entries[i],:value) for i=1:size(dimensions_rows)[1]]...)
diff --git a/src/system/entry.jl b/src/system/entry.jl
@@ -25,10 +25,10 @@ function Base.zero(::Type{Entry{ET}}) where ET
     return Entry{ET.parameters[2]}(dims...)
 end
 function Base.zero(::Type{Entry})
-    return 0
+    return nothing
 end
 
-function randomize!(entry::Entry, rand_function = randn)
+function initialize!(entry::Entry, init_function = randn)
     value = entry.value
-    entry.value = rand_function(eltype(value), size(value))
+    entry.value = init_function(eltype(value), size(value))
 end
diff --git a/src/system/interface.jl b/src/system/interface.jl
@@ -0,0 +1,54 @@
+@inline children(system::System, v) = outneighbors(system.dfs_graph, v) # all direct children of v
+@inline connections(system::System, v) = neighbors(system.graph, v)     # all connected nodes of v
+@inline parents(system::System, v) = inneighbors(system.dfs_graph, v)   # same elements as system.parents[v], but potentially different order
+
+LinearAlgebra.issymmetric(::System{N, Symmetric}) where N = true
+LinearAlgebra.issymmetric(::System{N, Unsymmetric}) where N = false
+
+function ranges(system::System{N}) where N
+    dims = system.dims
+    range_dict = Dict(1=>1:dims[1])
+    for i=2:N
+        range_dict[i] = last(range_dict[i-1])+1:sum(dims[1:i])
+    end
+
+    return range_dict
+end
+function ranges(dims::AbstractVector)
+    range_dict = Dict(1=>1:dims[1])
+    for i=2:size(dims)[1]
+        range_dict[i] = last(range_dict[i-1])+1:sum(dims[1:i])
+    end
+
+    return range_dict
+end
+
+function initialize!(system::System, init_function = randn)
+    for entry in system.matrix_entries.nzval
+        initialize!(entry, init_function)
+    end
+    for entry in system.vector_entries
+        initialize!(entry, init_function)
+    end
+end
+
+function initialize!(system::System{N,<:Symmetric}, init_function = randn) where N
+    matrix_entries = system.matrix_entries
+
+    for entry in matrix_entries.nzval
+        initialize!(entry, init_function)
+    end
+    for i=1:N
+        matrix_entries[i,i].value += matrix_entries[i,i].value'
+    end
+
+    for entry in system.vector_entries
+        initialize!(entry, init_function)
+    end
+end
+
+function reset_inverse_diagonals!(system::System)
+    for entry in system.diagonal_inverses
+        entry.isinverted = false
+    end
+end
diff --git a/src/system/system.jl b/src/system/system.jl
@@ -11,6 +11,7 @@ struct System{N,T}
     dfs_list::SVector{N,Int64}      # depth-first search list of nodes [last-found node, ..., first-found node]
     graph::SimpleGraph{Int64}       # the graph built from the adjacency matrix
     dfs_graph::SimpleDiGraph{Int64} # the directed graph built from the depth-first search
+    dims::Vector{Int64} # Dimensions of the matrix entries
 
     function System{T}(A, dims; force_static = false, symmetric = false) where T
         N = length(dims)
@@ -79,7 +80,7 @@ struct System{N,T}
         cyclic_children = [unique(vcat(cycles[i]...)) for i=1:N]
         cyclic_children = [intersect(dfs_list, cyclic_children[i]) for i=1:N]
 
-        new{N,S}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph)
+        new{N,S}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph, dims)
     end
 
     System(A, dims; force_static = false, symmetric = false) = System{Float64}(A, dims; force_static, symmetric)
@@ -94,79 +95,3 @@ function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, system::System{N,S}
     println(io, "System with "*string(N)*" nodes.")
     SparseArrays._show_with_braille_patterns(io, system.matrix_entries)
 end
-
-
-@inline children(system::System, v) = outneighbors(system.dfs_graph, v) # all direct children of v
-@inline connections(system::System, v) = neighbors(system.graph, v)     # all connected nodes of v
-@inline parents(system::System, v) = inneighbors(system.dfs_graph, v)   # same elements as system.parents[v], but potentially different order
-
-dimensions(system::System{N}) where N = [size(system.vector_entries[i].value)[1] for i=1:N]
-function ranges(system::System{N}) where N
-    dims = dimensions(system)
-    range_dict = Dict(1=>1:dims[1])
-    for i=2:N
-        range_dict[i] = last(range_dict[i-1])+1:sum(dims[1:i])
-    end
-
-    return range_dict
-end
-
-# There probably exists a smarter way of getting the dense matrix from the spares one 
-function full_matrix(system::System{N}) where N
-    dims = dimensions(system)
-    range_dict = ranges(system)
-    A = zeros(sum(dims), sum(dims))
-
-    for (i,row) in enumerate(system.matrix_entries.rowval)
-        col = findfirst(x -> i<x, system.matrix_entries.colptr)-1
-        A[range_dict[row],range_dict[col]] = system.matrix_entries[row,col].value
-    end
-    return A
-end
-
-function full_matrix(system::System{N,<:Symmetric}) where N
-    dims = dimensions(system)
-    range_dict = ranges(system)
-    A = zeros(sum(dims),sum(dims))
-
-    for (i,row) in enumerate(system.matrix_entries.rowval)
-        col = findfirst(x -> i<x, system.matrix_entries.colptr)-1
-        A[range_dict[row],range_dict[col]] = system.matrix_entries[row,col].value
-        if col != row
-            A[range_dict[col],range_dict[row]] = system.matrix_entries[row,col].value'
-        end
-    end
-    return A
-end
-
-full_vector(system::System) = vcat(getfield.(system.vector_entries,:value)...)
-
-function randomize!(system::System, rand_function = randn)
-    for entry in system.matrix_entries.nzval
-        randomize!(entry, rand_function)
-    end
-    for entry in system.vector_entries
-        randomize!(entry, rand_function)
-    end
-end
-
-function randomize!(system::System{N,<:Symmetric}, rand_function = randn) where N
-    matrix_entries = system.matrix_entries
-
-    for entry in matrix_entries.nzval
-        randomize!(entry, rand_function)
-    end
-    for i=1:N
-        matrix_entries[i,i].value += matrix_entries[i,i].value'
-    end
-
-    for entry in system.vector_entries
-        randomize!(entry, rand_function)
-    end
-end
-
-function reset_inverse_diagonals!(system::System)
-    for entry in system.diagonal_inverses
-        entry.isinverted = false
-    end
-end
diff --git a/test/interface.jl b/test/interface.jl
@@ -0,0 +1,27 @@
+using GraphBasedSystems
+using LinearAlgebra
+
+include("adjacency_matrix.jl")
+
+system = System{Float64}(A, rand(0:3, size(A)[1]))
+
+children(system, 1)
+connections(system, 1)
+parents(system, 1)
+children(system, 10)
+connections(system, 10)
+parents(system, 10)
+@test true
+
+@test GraphBasedSystems.ranges(system) == GraphBasedSystems.ranges(system.dims)
+
+initialize!(system, rand)
+initialize!(system, randn)
+initialize!(system, ones)
+initialize!(system, zeros)
+
+@test true
+
+reset_inverse_diagonals!(system)
+
+@test true
diff --git a/test/ldlt_test.jl b/test/ldlt_test.jl
@@ -7,7 +7,7 @@ include("adjacency_matrix.jl")
 system = System{Float64}(A, rand(0:3, size(A)[1]), symmetric=true)
 
 for i=1:10
-    randomize!(system)
+    initialize!(system)
 
     F = full_matrix(system)
     f = full_vector(system)
diff --git a/test/ldu_test.jl b/test/ldu_test.jl
@@ -7,7 +7,7 @@ include("adjacency_matrix.jl")
 system = System{Float64}(A, rand(0:3, size(A)[1]))
 
 for i=1:10
-    randomize!(system)
+    initialize!(system)
 
     F = full_matrix(system)
     f = full_vector(system)
diff --git a/test/llt_test.jl b/test/llt_test.jl
diff --git a/test/lu_test.jl b/test/lu_test.jl
diff --git a/test/runtests.jl b/test/runtests.jl
