diff --git a/Project.toml b/Project.toml index 4c97cada..2d8af182 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.7.1" [deps] AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" @@ -22,7 +23,6 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" AMD = "0.4, 0.5" GenericLinearAlgebra = "0.3" MathOptInterface = "1.20" -QDLDL = "0.4" Requires = "1" SnoopPrecompile = "1" StaticArrays = "1" diff --git a/src/Clarabel.jl b/src/Clarabel.jl index 2e099a31..8a782c92 100644 --- a/src/Clarabel.jl +++ b/src/Clarabel.jl @@ -4,7 +4,9 @@ module Clarabel using SparseArrays, LinearAlgebra, Printf, Requires const DefaultFloat = Float64 const DefaultInt = LinearAlgebra.BlasInt - const IdentityMatrix = UniformScaling{Bool} + + # Rust-like Option type + const Option{T} = Union{Nothing,T} #internal constraint RHS limits. This let block #hides the INFINITY field in the module and makes @@ -29,11 +31,12 @@ module Clarabel include("./cones/compositecone_type.jl") #core solver components + include("./abstract_types.jl") include("./settings.jl") - include("./conicvector.jl") include("./statuscodes.jl") + include("./chordal/include.jl") + include("./types.jl") include("./presolver.jl") - include("./types.jl") include("./variables.jl") include("./residuals.jl") include("./problemdata.jl") diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index a60007ba..9d93aa2c 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -495,7 +495,7 @@ function process_constraints( rowranges = dest.rowranges m = mapreduce(length, +, values(rowranges), init=0) - b = Vector{T}(undef, m) + b = zeros(T, m) #these will be used for a triplet representation of A I = Int[] diff --git a/src/abstract_types.jl b/src/abstract_types.jl new file mode 100644 index 00000000..6bc4d075 --- /dev/null +++ b/src/abstract_types.jl @@ -0,0 +1,9 @@ +abstract type AbstractVariables{T <: AbstractFloat} end +abstract type AbstractEquilibration{T <: AbstractFloat} end +abstract type AbstractResiduals{T <: AbstractFloat} end +abstract type AbstractProblemData{T <: AbstractFloat} end +abstract type AbstractKKTSystem{T <: AbstractFloat} end +abstract type AbstractKKTSolver{T <: AbstractFloat} end +abstract type AbstractInfo{T <: AbstractFloat} end +abstract type AbstractSolution{T <: AbstractFloat} end +abstract type AbstractSolver{T <: AbstractFloat} end \ No newline at end of file diff --git a/src/chordal/chordal_info.jl b/src/chordal/chordal_info.jl new file mode 100644 index 00000000..033894d5 --- /dev/null +++ b/src/chordal/chordal_info.jl @@ -0,0 +1,304 @@ +# ------------------------------------- +# Chordal Decomposition Information +# ------------------------------------- +struct ConeMapEntry + orig_index::Int + tree_and_clique::Option{Tuple{Int,Int}} +end + +mutable struct ChordalInfo{T} + + # sketch of the original problem + init_dims::Tuple{Int64, Int64} # (n,m) dimensions of the original problem + init_cones::Vector{SupportedCone} # original cones of the problem + + # decomposed problem data + spatterns::Vector{SparsityPattern} # sparsity patterns for decomposed cones + + # `H' matrix for the standard chordal problem transformation + # remains as nothing if the `compact' transform is used + H::Option{SparseMatrixCSC{T}} + + # mapping from each generated cone to its original cone + # index, plus its tree and clique information if it has + # been generated as part of a chordal decomposition + # remains as nothing if the `standard' transform is used + cone_maps::Option{Vector{ConeMapEntry}} + + function ChordalInfo( + A::SparseMatrixCSC{T}, + b::Vector{T}, + cones::Vector{SupportedCone}, + settings::Clarabel.Settings{T} + ) where {T} + + # initial problem data + init_dims = (size(A,2),size(A,1)) + + chordal_info = new{T}( + init_dims, + SupportedCone[], + SparsityPattern[], + nothing, # no H to start + nothing # no cone_maps to start + ) + + find_sparsity_patterns!( + chordal_info, + A, b, + cones, + settings.chordal_decomposition_merge_method) + + # Only copy the generating cones if we have decomposition, + # since otherwise this object is going to be dropped anyway + if is_decomposed(chordal_info) + chordal_info.init_cones = deepcopy(cones) + end + + return chordal_info + + end +end + + +function find_sparsity_patterns!( + chordal_info::ChordalInfo{T}, + A::SparseMatrixCSC{T}, + b::Vector{T}, + cones::Vector{SupportedCone}, + merge_method::Symbol +) where {T} + + rng_cones = rng_cones_iterator(cones) + + # aggregate sparsity pattern across the rows of [A;b] + nz_mask = find_aggregate_sparsity_mask(A, b) + + # find the sparsity patterns of the PSD cones + for (coneidx, (cone,rowrange)) in enumerate(zip(cones,rng_cones)) + + if isa(cone, PSDTriangleConeT) + analyse_psdtriangle_sparsity_pattern!( + chordal_info, + view(nz_mask,rowrange), + cone.dim, + coneidx, + merge_method) + end + end +end + + +function analyse_psdtriangle_sparsity_pattern!( + chordal_info::ChordalInfo{T}, + nz_mask::AbstractVector{Bool}, + conedim::Int, + coneidx::Int, + merge_method::Symbol +) where {T} + + # Force the diagonal entries to be marked, otherwise + # the symbolic LDL step will fail. + for i = 1:conedim + nz_mask[triangular_index(i)] = true + end + + if all(nz_mask) + return #dense / decomposable + end + + L, ordering = find_graph!(nz_mask) + spattern = SparsityPattern(L, ordering, coneidx, merge_method) + + if spattern.sntree.n_cliques == 1 + return #not decomposed, or everything re-merged + end + + push!(chordal_info.spatterns, spattern) + +end + + +# did any PSD cones get decomposed? +function is_decomposed(chordal_info::ChordalInfo{T}) where {T} + return !isempty(chordal_info.spatterns) +end + + +# total number of cones we started with +function init_cone_count(chordal_info::ChordalInfo{T}) where {T} + length(chordal_info.init_cones) +end + +# total number of PSD cones we started with +function init_psd_cone_count(chordal_info::ChordalInfo{T}) where {T} + count(c -> isa(c,PSDTriangleConeT), chordal_info.init_cones) +end + +# total number of cones of all types after decomp and merge +function final_cone_count(chordal_info::ChordalInfo{T}) where {T} + return init_cone_count(chordal_info) + final_psd_cones_added(chordal_info) +end + +# total number of PSD cones after decomp and merge +function final_psd_cone_count(chordal_info::ChordalInfo{T}) where {T} + return init_psd_cone_count(chordal_info) + final_psd_cones_added(chordal_info) +end + +# total number of PSD cones after decomp and merge +function premerge_psd_cone_count(chordal_info::ChordalInfo{T}) where {T} + return init_psd_cone_count(chordal_info) + premerge_psd_cones_added(chordal_info) +end + +# total number of PSD cones that were decomposable +function decomposable_cone_count(chordal_info::ChordalInfo{T}) where {T} + length(chordal_info.spatterns) +end + +# Net PSD cones added after decomposition and merging +function final_psd_cones_added(chordal_info::ChordalInfo{T}) where {T} + # sum the number of cliques in each spattern + ncliques = sum([spattern.sntree.n_cliques for spattern in chordal_info.spatterns]) + ndecomposable = decomposable_cone_count(chordal_info) + #subtract to avoid double counting + ncliques - ndecomposable +end + +# PSD cones added by decomposition, and before merging +function premerge_psd_cones_added(chordal_info::ChordalInfo{T}) where {T} + # sum the number of cliques in each spattern + ncones = sum([length(spattern.sntree.snode) for spattern in chordal_info.spatterns]) + ndecomposable = decomposable_cone_count(chordal_info) + #subtract to avoid double counting + ncones - ndecomposable +end + + + +function get_decomposed_dim_and_overlaps(chordal_info::ChordalInfo{T}) where{T} + + cones = chordal_info.init_cones + sum_cols = 0 + sum_overlaps = 0 + patterns_iter = Iterators.Stateful(chordal_info.spatterns) + + for (coneidx,cone) in enumerate(cones) + cols, overlap = let + if !isempty(patterns_iter) && peek(patterns_iter).orig_index == coneidx + get_decomposed_dim_and_overlaps(first(patterns_iter).sntree) + else + nvars(cone), 0 + end + end + sum_cols += cols + sum_overlaps += overlap + end + + sum_cols, sum_overlaps + +end + + + +# ------------------------------------- +# FUNCTION DEFINITIONS +# ------------------------------------- + +function find_aggregate_sparsity_mask( + A::SparseMatrixCSC{T}, + b::Vector{T}, +) where {T <: AbstractFloat} + + # returns true in every row in which [A;b] has a nonzero + + active = falses(length(b)) + active[A.rowval] .= true + @. active |= (b != zero(T)) + active +end + + +function find_graph!(nz_mask::AbstractVector{Bool}) + + nz = count(nz_mask) + rows = sizehint!(Int[], nz) + cols = sizehint!(Int[], nz) + + # check final row/col to get matrix dimension + (m,n) = upper_triangular_index_to_coord(length(nz_mask)) + @assert m == n + + for (linearidx, isnonzero) in enumerate(nz_mask) + if isnonzero + (row,col) = upper_triangular_index_to_coord(linearidx) + push!(rows, row) + push!(cols, col) + end + end + + # QDLDL doesn't currently allow for logical-only decomposition + # on a matrix of Bools, so pattern must be a Float64 matrix here + pattern = sparse(rows, cols, ones(Float64, length(rows)), n, n) + + F = QDLDL.qdldl(pattern, logical = true) + + L = F.L + ordering = F.perm + + # this takes care of the case that QDLDL returns an unconnected adjacency matrix L + connect_graph!(L) + + return L, ordering +end + + +function connect_graph!(L::SparseMatrixCSC{Float64}) + + # unconnected blocks don't have any entries below the diagonal in their right-most columns + n = size(L, 2) + for j = 1:n-1 + + row_val = L.rowval + col_ptr = L.colptr + + connected = false + for k in col_ptr[j]:col_ptr[j+1]-1 + if row_val[k] > j + connected = true + break + end + end + + #this insertion can happen in a midrange column, as long as + #that column is the last one for a given block + if !connected + L[j+1, j] = 1 + end + end + +end + + +# ---------------------------------------------- +# Iterator for the range of indices of the cones + +struct RangeSupportedConesIterator + cones::Vector{SupportedCone} +end + +function rng_cones_iterator(cones::Vector{SupportedCone}) + RangeSupportedConesIterator(cones) +end + +Base.length(C::RangeSupportedConesIterator) = length(C.cones) + +function Base.iterate(C::RangeSupportedConesIterator, state=(1, 1)) + (coneidx, start) = state + if coneidx > length(C.cones) + return nothing + else + stop = start + nvars(C.cones[coneidx]) - 1 + state = (coneidx + 1, stop + 1) + return (start:stop, state) + end +end diff --git a/src/chordal/decomposition/augment_compact.jl b/src/chordal/decomposition/augment_compact.jl new file mode 100644 index 00000000..1fe38864 --- /dev/null +++ b/src/chordal/decomposition/augment_compact.jl @@ -0,0 +1,521 @@ +# ----------------------------------------- +# Functions related to clique tree based transformation (compact decomposition) +# see: Kim: Exploiting sparsity in linear and nonlinear matrix inequalities via +# positive semidefinite matrix completion (2011), p.53 +# ----------------------------------------- + +const BlockOverlapTriplet = Tuple{Int, Int, Bool} + +function decomp_augment_compact!( + chordal_info::ChordalInfo, + P::SparseMatrixCSC{T}, + q::Vector{T}, + A::SparseMatrixCSC{T}, + b::Vector{T}, + ) where {T} + + A_new, b_new, cones_new = find_compact_A_b_and_cones(chordal_info, A, b) + + # how many variables did we add? + nadd = A_new.n - A.n + + P_new = blockdiag(P, spzeros(T, nadd, nadd)) + + q_new = zeros(T, length(q) + nadd) + q_new[1:length(q)] .= q + + return P_new, q_new, A_new, b_new, cones_new + +end + +function find_compact_A_b_and_cones( + chordal_info::ChordalInfo{T}, + A::SparseMatrixCSC{T}, + b::Vector{T}, +) where{T} + + # the cones that we used to form the decomposition + cones = chordal_info.init_cones + + # determine number of final augmented matrix and number of overlapping entries + Aa_m, Aa_n, n_overlaps = find_A_dimension(chordal_info, A) + + # allocate sparse components for the augmented A + Aa_nnz = nnz(A) + 2 * n_overlaps + Aa_I = zeros(Int, Aa_nnz) + Aa_J = extra_columns(Aa_nnz, nnz(A) + 1, A.n + 1) + Aa_V = alternating_sequence(T, Aa_nnz, nnz(A) + 1) + + findnz!(Aa_J, Aa_V, A) + + # allocate sparse components for the augmented b + bs = sparse(b) + ba_I = zeros(Int, length(bs.nzval)) + ba_V = bs.nzval + + # preallocate the decomposed cones and the mapping + # from decomposed cones back to the originals + n_decomposed = final_cone_count(chordal_info) + cones_new = sizehint!(SupportedCone[], n_decomposed) + cone_maps = sizehint!(ConeMapEntry[], n_decomposed) + + # an enumerate-like mutable iterator for the patterns. We will expand cones + # assuming that they are non-decomposed until we reach an index that agrees + # the internally stored orig_index of the next pattern. + + patterns_iter = Iterators.Stateful(chordal_info.spatterns) + patterns_count = Iterators.Stateful(eachindex(chordal_info.spatterns)) + row_ranges = rng_cones_iterator(cones) + + row_ptr = 1 # index to start of next cone in A_I + overlap_ptr = nnz(A) + 1 # index to next row for +1, -1 overlap entries + + for (coneidx, (cone, row_range)) in enumerate(zip(cones,row_ranges)) + + if !isempty(patterns_iter) && peek(patterns_iter).orig_index == coneidx + + @assert(isa(cone, PSDTriangleConeT)) + row_ptr, overlap_ptr = add_entries_with_sparsity_pattern!( + Aa_I, ba_I, cones_new, cone_maps, A, bs, row_range, + first(patterns_iter), first(patterns_count), + row_ptr, overlap_ptr) + else + row_ptr, overlap_ptr = add_entries_with_cone!( + Aa_I, ba_I, cones_new, cone_maps, A, bs, row_range, cone, row_ptr, overlap_ptr) + end + + end + + # save the cone_maps for use when reconstructing + # solution to the original problem + chordal_info.cone_maps = cone_maps + + A_new = allocate_sparse_matrix(Aa_I, Aa_J, Aa_V, Aa_m, Aa_n) + b_new = Vector(SparseArrays._sparsevector!(ba_I, ba_V, Aa_m)) + + A_new, b_new, cones_new +end + +# find the dimension of the `compact' form `A' matrix and its number of overlaps +function find_A_dimension(chordal_info::ChordalInfo{T}, A::SparseMatrixCSC{T}) where {T} + + dim, num_overlaps = get_decomposed_dim_and_overlaps(chordal_info) + + rows = dim + cols = A.n + num_overlaps + + return rows, cols, num_overlaps + +end + + +# Given the row, column, and nzval vectors and dimensions, assembles the sparse matrix `Aa` +# of the decomposed problem in a slightly more memory efficient way. +function allocate_sparse_matrix( + Aa_I::Vector{Int}, + Aa_J::Vector{Int}, + Aa_V::Vector{T}, + mA::Int, nA::Int +) where {T} + csrrowptr = zeros(Int, mA + 1) + csrcolval = zeros(Int, length(Aa_I)) + csrnzval = zeros(T, length(Aa_I)) + klasttouch = zeros(Int, nA + 1) + csccolptr = zeros(Int, nA + 1) + + SparseArrays.sparse!(Aa_I, Aa_J, Aa_V, mA, nA, +, klasttouch, csrrowptr, csrcolval, csrnzval, csccolptr, Aa_I, Aa_V ) + +end + + +# Handles all cones that are not decomposed by a sparsity pattern +function add_entries_with_cone!( + Aa_I::Vector{Int}, + ba_I::Vector{Int}, + cones_new::Vector{SupportedCone}, + cone_maps::Vector{ConeMapEntry}, + A::SparseMatrixCSC{T}, + b::SparseVector{T}, + row_range::UnitRange{Int}, + cone::SupportedCone, + row_ptr::Int, overlap_ptr::Int +) where {T} + + n = A.n + offset = row_ptr - row_range.start + + + # populate b + row_range_col = get_rows_vec(b, row_range) + if !isnothing(row_range_col) + for k in row_range_col + ba_I[k] = b.nzind[k] + offset + end + end + + # populate A + for col = 1:n + # indices that store the rows in column col in A + row_range_col = get_rows_mat(A, col, row_range) + if !isnothing(row_range_col) + for k in row_range_col + Aa_I[k] = A.rowval[k] + offset + end + end + end + + # here we make a copy of the cone + push!(cones_new, deepcopy(cone)) + + # since this cone is standalone and not decomposed, the index + # of its origin cone must be either one more than the previous one, + # or 1 (zero in rust) if it's the first + + orig_index = isempty(cone_maps) ? 1 : cone_maps[end].orig_index + 1 + push!(cone_maps, ConeMapEntry(orig_index, nothing)) + + return row_ptr + nvars(cone), overlap_ptr + +end + + + +# Handle decomposable cones with a SparsityPattern. The row vectors A_I and b_I +# have to be edited in such a way that entries of one clique appear contiguously. +function add_entries_with_sparsity_pattern!( + A_I::Vector{Int}, + b_I::Vector{Int}, + cones_new::Vector{SupportedCone}, + cone_maps::Vector{ConeMapEntry}, + A::SparseMatrixCSC{T}, + b::SparseVector{T}, + row_range::UnitRange{Int}, + spattern::SparsityPattern, + spattern_index::Int, + row_ptr::Int, overlap_ptr::Int +) where {T} + + sntree = spattern.sntree + ordering = spattern.ordering + + (_, n) = size(A) + + # determine the row ranges for each of the subblocks + clique_to_rows = clique_rows_map(row_ptr, sntree) + + # loop over cliques in descending topological order + for i in (sntree.n_cliques):-1:1 + + # get supernodes and separators and undo the reordering + # NB: these are now Vector, not VertexSet + separator = sort!([spattern.ordering[v] for v in get_separators(sntree, i)]) + snode = sort!([spattern.ordering[v] for v in get_snode(sntree, i)]) + + # compute sorted block indices (i, j, flag) for this clique with an + # information flag whether an entry (i, j) is an overlap + block_indices = get_block_indices(snode, separator, length(ordering)) + + # If we encounter an overlap with a parent clique we have to be able to find the + # location of the overlapping entry. Therefore load and reorder the parent clique + if i == sntree.n_cliques + parent_rows = 0:0 + parent_clique = Int[] + else + parent_index = get_clique_parent(sntree, i) + parent_rows = clique_to_rows[parent_index] + parent_clique = [spattern.ordering[v] for v in get_clique_by_index(sntree, parent_index)] + sort!(parent_clique) + end + + # Loop over all the columns and shift the rows in A_I and b_I according to the clique structure + for col = 1:n + + row_range_col = get_rows_mat(A, col, row_range) + row_range_b = col == 1 ? get_rows_vec(b, row_range) : 0:0 + + #PJG this method for producing empty ranges is gross + #should just be Nothing + row_range_col = isnothing(row_range_col) ? (1:0) : row_range_col + row_range_b = isnothing(row_range_b) ? (1:0) : row_range_b + + overlap_ptr = add_clique_entries!( + A_I, b_I, A.rowval, b.nzind, block_indices, + parent_clique, parent_rows, col, + row_ptr, overlap_ptr, row_range, row_range_col, row_range_b) + + end + + # create new PSD cones for the subblocks, and tag them + # with their tree and clique number + + cone_dim = get_nblk(sntree, i) + push!(cones_new, PSDTriangleConeT(cone_dim)) + push!(cone_maps, ConeMapEntry(spattern.orig_index, (spattern_index, i))) + + row_ptr += triangular_number(cone_dim) + + end + + return row_ptr, overlap_ptr + +end + + + +# Loop over all entries (i, j) in the clique and either set the correct row in `A_I` and `b_I` +# if (i, j) is not an overlap,or add an overlap column with (-1 and +1) in the correct positions. + +function add_clique_entries!( + A_I::Vector{Int}, + b_I::Vector{Int}, + A_rowval::Vector{Int}, + b_nzind::Vector{Int}, + block_indices::Vector{BlockOverlapTriplet}, + parent_clique::Vector{Int}, + parent_rows::UnitRange{Int}, + col::Int, + row_ptr::Int, + overlap_ptr::Int, + row_range::UnitRange{Int}, + row_range_col::UnitRange{Int}, + row_range_b::UnitRange{Int} +) + + counter = 0 + for block_idx in block_indices + new_row_val = row_ptr + counter + (i,j,is_overlap) = block_idx + # a block index that corresponds to an overlap + if is_overlap + if col == 1 + # this creates the +1 entry + A_I[overlap_ptr] = new_row_val + # this creates the -1 entry + A_I[overlap_ptr + 1] = parent_rows.start + parent_block_indices(parent_clique, i, j) - 1 + overlap_ptr += 2 + end + else + k = coord_to_upper_triangular_index((i, j)) + modify_clique_rows!(A_I, k, A_rowval, new_row_val, row_range, row_range_col) + col == 1 && modify_clique_rows!(b_I, k, b_nzind, new_row_val, row_range, row_range_b) + + end + counter += 1 + + end + return overlap_ptr +end + + +# Given the nominal entry position `k = linearindex(i, j)` find and modify with `new_row_val` +# the actual location of that entry in the global row vector `rowval`. + +function modify_clique_rows!( + v::Vector{Int}, + k::Int, + rowval::Vector{Int}, + new_row_val::Int, + row_range::UnitRange{Int}, + row_range_col::UnitRange{Int} +) + + row_0 = get_row_index(k, rowval, row_range, row_range_col) + # row_0 happens when (i, j) references an edge that was added by merging cliques, + # the corresponding value will be zero and can be disregarded + if row_0 != 0 + v[row_0] = new_row_val + end + return nothing +end + + +# Given the svec index `k` and an offset `row_range_col.start`, return the location of the +# (i, j)th entry in the row vector `rowval`. +function get_row_index( + k::Int, + rowval::Vector{Int}, + row_range::UnitRange{Int}, + row_range_col::UnitRange{Int} + ) + + #PJG: possible logic error here, since we don't pass down a Union type + #this far. Fix this and use Option everywhere, or short circuit + #higher up the call stack. Should also return Nothing instead of 0 + isnothing(row_range_col) && return 0 + + k_shift = row_range.start + k - 1 + + # determine upper set boundary of where the row could be + u = min(row_range_col.stop, row_range_col.start + k_shift - 1) + + # find index of first entry >= k, starting in the interval [l, u] + # if no, entry is >= k, returns u + 1 + r = searchsortedfirst(rowval, k_shift, row_range_col.start, u, Base.Order.Forward) + + # if no r s.t. rowval[r] = k_shift was found, that means that the + # (i, j)th entry represents an edded edge (zero) from clique merging + if r > u || rowval[r] != k_shift + return 0 + else + return r + end + + +end + + + +# Find the index of k=svec(i, j) in the parent clique `par_clique`.# + +function parent_block_indices(parent_clique::Vector{Int}, i::Int, j::Int) + ir = searchsortedfirst(parent_clique, i) + jr = searchsortedfirst(parent_clique, j) + return coord_to_upper_triangular_index((ir, jr)) +end + + + +# Given a cliques supernodes and separators, compute all the indices (i, j) of the corresponding matrix block +# in the format (i, j, flag), where flag is equal to false if entry (i, j) corresponds to an overlap of the +# clique and true otherwise. + +# `nv` is the number of vertices in the graph that we are trying to decompose. + +function get_block_indices(snode::Array{Int}, separator::Array{Int}, nv::Int) + + N = length(separator) + length(snode) + + block_indices = sizehint!(BlockOverlapTriplet[],triangular_number(N)) + + for j in separator, i in separator + if i <= j + push!(block_indices,(i, j, true)) + end + end + + for j in snode, i in snode + if i <= j + push!(block_indices,(i, j, false)) + end + end + + for i in snode + for j in separator + push!(block_indices,(min(i, j), max(i, j), false)) + end + end + + sort!(block_indices, by = x -> x[2] * nv + x[1] ) + return block_indices +end + + + +# Return the row ranges of each clique after the decomposition, shifted by `row_start`. + +function clique_rows_map(row_start::Int, sntree::SuperNodeTree) + + n_cliques = sntree.n_cliques + + # PJG: rows / inds not necessary here. Should just size hint + # on the output Dict and push the entries directly to it + rows = sizehint!(UnitRange{Int}[], n_cliques) + inds = sizehint!(Int[], n_cliques) + + for i in n_cliques:-1:1 + num_rows = triangular_number(get_nblk(sntree, i)) + push!(rows, row_start:row_start+num_rows-1) + push!(inds, sntree.snode_post[i]) + row_start += num_rows + end + + return Dict(inds .=> rows) +end + +function get_rows_subset(rows, row_range) + + if length(rows) == 0 + return nothing + end + + s = searchsortedfirst(rows, row_range.start) + if s == length(rows) + 1 + return nothing + end + + if rows[s] > row_range.stop || s == 0 + return nothing + else + e = searchsortedlast(rows, row_range.stop) + return s:e + end + +end + +function get_rows_vec(b::SparseVector, row_range::UnitRange{Int}) + + get_rows_subset(b.nzind, row_range) + + end + + function get_rows_mat(A::SparseMatrixCSC, col::Int, row_range::UnitRange{Int}) + + colrange = A.colptr[col]:(A.colptr[col + 1]-1) + rows = view(A.rowval, colrange) + se = get_rows_subset(rows, row_range) + + if isnothing(se) + return nothing + end + + colrange[se.start]:colrange[se.stop] + + end + + +# Returns the appropriate amount of memory for `A.nzval`, including, starting from `n_start`, +# the (+1 -1) entries for the overlaps. + +function alternating_sequence(T, total_length::Int, n_start::Int) + v = ones(T, total_length) + for i= (n_start + 1):2:length(v) + v[i] = -one(T) + end + return v +end + + +# Returns the appropriate amount of memory for the columns of the augmented problem matrix `A`, +# including, starting from `n_start`, the columns for the (+1 -1) entries for the overlaps. + +function extra_columns(total_length::Int, n_start::Int, start_val::Int) + v = zeros(Int, total_length) + for i = n_start:2:length(v)-1 + v[i] = start_val + v[i+1] = start_val + start_val += 1 + end + return v +end + + +# Given sparse matrix components, write the columns and non-zero values into the first `numnz` entries of `J` and `V`. + +function findnz!(J::Vector{Ti}, V::Vector{Tv}, S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} + count = 1 + @inbounds for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1) + J[count] = col + V[count] = S.nzval[k] + count += 1 + end +end + + +# Intentionally defined here separately from the other SuperNodeTree functions. +# This returns a clique directly from index i, rather than accessing the +# snode and separators via the postordering. Only used in the compact +# problem decomposition functions + +function get_clique_by_index(sntree::SuperNodeTree, i::Int) + return union(sntree.snode[i], sntree.separators[i]) +end diff --git a/src/chordal/decomposition/augment_standard.jl b/src/chordal/decomposition/augment_standard.jl new file mode 100644 index 00000000..5668b90a --- /dev/null +++ b/src/chordal/decomposition/augment_standard.jl @@ -0,0 +1,138 @@ +# ----------------------------------------- +# Functions related to `traditional' decomposition +# ----------------------------------------- + +function decomp_augment_standard!( + chordal_info::ChordalInfo{T}, + P::SparseMatrixCSC{T}, + q::Vector{T}, + A::SparseMatrixCSC{T}, + b::Vector{T}, +) where {T} + + # allocate H and new decomposed cones. H will be stored + # in chordal_info and is required for the reversal step + H, cones_new = find_standard_H_and_cones!(chordal_info) + + P_new = blockdiag(P, spzeros(T, H.n, H.n)) + + q_new = zeros(T, length(q) + H.n) + q_new[1:length(q)] .= q + + A_new = [A H; spzeros(T, H.n, A.n) -sparse(1.0I, H.n, H.n)] + + b_new = zeros(T, length(b) + H.n) + b_new[1:length(b)] .= b + + # save the H matrix for use when reconstructing + # solution to the original problem + chordal_info.H = H + + return P_new, q_new, A_new, b_new, cones_new + +end + + +# Find the transformation matrix `H` and its associated cones for the standard decomposition. + +function find_standard_H_and_cones!( + chordal_info::ChordalInfo{T}, +) where {T} + + # the cones that we used to form the decomposition + cones = chordal_info.init_cones + + # preallocate H and new decomposed cones + lenH = find_H_col_dimension(chordal_info) + H_I = sizehint!(Int[], lenH) + + # ncones from decomposition, plus one for an additional equality constraint + cones_new = sizehint!(SupportedCone[], final_cone_count(chordal_info) + 1) + + # +1 cone count above is for this equality constraint + (_,m) = chordal_info.init_dims + push!(cones_new, ZeroConeT(m)) + + # an iterator for the patterns. We will expand cones assuming that + # they are non-decomposed until we reach an index that agrees with + # the internally stored coneidx of the next pattern. + patterns_iter = Iterators.Stateful(chordal_info.spatterns) + row = 1 + + for (coneidx,cone) in enumerate(cones) + + if !isempty(patterns_iter) && peek(patterns_iter).orig_index == coneidx + @assert(isa(cone, PSDTriangleConeT)) + decompose_with_sparsity_pattern!(H_I, cones_new, first(patterns_iter), row) + else + decompose_with_cone!(H_I, cones_new, cone, row) + end + + row += nvars(cone) + end + + H = sparse(H_I, collect(1:lenH), ones(T,lenH)) + + return H, cones_new +end + +function find_H_col_dimension(chordal_info::ChordalInfo{T}) where {T} + + cols, _ = get_decomposed_dim_and_overlaps(chordal_info) + return cols +end + + +function decompose_with_cone!( + H_I::Vector{Int}, + cones_new::Vector{SupportedCone}, + cone::SupportedCone, + row::Int +) + for i = 1:nvars(cone) + push!(H_I,row + i - 1) + end + + push!(cones_new, deepcopy(cone)) +end + + +function decompose_with_sparsity_pattern!( + H_I::Vector{Int}, + cones_new::Vector{SupportedCone}, + spattern::SparsityPattern, + row::Int +) + + sntree = spattern.sntree + + for i in 1:sntree.n_cliques + + clique = get_clique(sntree, i) + + # the graph and tree algorithms determined the clique vertices of an + # AMD-permuted matrix. Since the location of the data hasn't changed + # in reality, we have to map the clique vertices back + c = sort!([spattern.ordering[v] for v in clique]) + + add_subblock_map!(H_I, c, row) + + # add a new cone for this subblock + cdim = get_nblk(sntree, i) + push!(cones_new, PSDTriangleConeT(cdim)) + end + +end + + +function add_subblock_map!(H_I::Vector{Int}, clique_vertices::Vector{Int}, row_start::Int) + + v = clique_vertices + + for j in eachindex(v), i in 1:j + row = coord_to_upper_triangular_index((v[i], v[j])) + push!(H_I, row_start + row - 1) + end +end + + diff --git a/src/chordal/decomposition/decomp.jl b/src/chordal/decomposition/decomp.jl new file mode 100644 index 00000000..cb008e29 --- /dev/null +++ b/src/chordal/decomposition/decomp.jl @@ -0,0 +1,59 @@ +function decomp_augment!( + chordal_info::ChordalInfo{T}, + P::SparseMatrixCSC{T}, + q::Vector{T}, + A::SparseMatrixCSC{T}, + b::Vector{T}, + settings::Settings{T} +) where{T} + + if settings.chordal_decomposition_compact + decomp_augment_compact!(chordal_info, P, q, A, b) + else + decomp_augment_standard!(chordal_info, P, q, A, b) + end + +end + + +#AbstractVariables here should be DefaultVariables only, but this +# is not enforced since it causes a circular inclusion issue. + +function decomp_reverse!( + chordal_info::ChordalInfo{T}, + old_vars::AbstractVariables{T}, + old_cones::Vector{SupportedCone}, + settings::Settings{T} +) where{T} + + # We should have either H (for standard decomp) or cone_maps (for compact decomp) + # but never both, and they should be consistent with the settings + @assert settings.chordal_decomposition_compact == isnothing(chordal_info.H) + @assert settings.chordal_decomposition_compact != isnothing(chordal_info.cone_maps) + + # here `old_cones' are the ones that were used internally + # in the solver, producing internal solution in `old_vars' + # the cones for the problem as provided by the user or the + # upstream presolver are held internally in chordal_info.cones + + (n,m) = chordal_info.init_dims + new_vars = typeof(old_vars)(n,m) + + new_vars.x .= old_vars.x[1:n] + + # reassemble the original variables s and z + if settings.chordal_decomposition_compact + decomp_reverse_compact!( + chordal_info, new_vars, old_vars, old_cones) + else + decomp_reverse_standard!( + chordal_info, new_vars, old_vars, old_cones) + end + + if settings.chordal_decomposition_complete_dual + psd_completion!(chordal_info, new_vars) + end + + return new_vars +end + diff --git a/src/chordal/decomposition/psd_completion.jl b/src/chordal/decomposition/psd_completion.jl new file mode 100644 index 00000000..d72e763c --- /dev/null +++ b/src/chordal/decomposition/psd_completion.jl @@ -0,0 +1,106 @@ +# ----------------------------------- +# psd completion +# ----------------------------------- + +# The psd entries of z that correspond to the zeros in s are not constrained by the problem. +# however, in order to make the dual psd cone positive semidefinite we have to do a +# positive semidefinite completion routine to choose the values + +function psd_completion!( + chordal_info::ChordalInfo{T}, + variables::AbstractVariables{T} +) where {T} + + # working now with the cones from the original + # problem, not the decomposed ones + cones = chordal_info.init_cones + + # loop over psd cones + row_ranges = collect(rng_cones_iterator(cones)) + + # loop over just the patterns + for pattern in chordal_info.spatterns + row_range = row_ranges[pattern.orig_index] + z = @view variables.z[row_range] + complete!(z, pattern) + end + + return nothing +end + +function complete!(z::AbstractVector{T},pattern::SparsityPattern) where{T} + + n = length(pattern.ordering) + Z = zeros(n,n) + svec_to_mat!(Z,z) + psd_complete!(Z,pattern) + mat_to_svec!(z,Z) + +end + +# positive semidefinite completion (from Vandenberghe - Chordal Graphs..., p. 362) +# input: A - positive definite completable matrix +function psd_complete!(A::AbstractMatrix{T}, pattern::SparsityPattern) where {T <: AbstractFloat} + + sntree = pattern.sntree + p = pattern.ordering + ip = invperm(p) + N = size(A,2) + + # PJG: not clear if this copy of A is required, or + # whether I can operate directly on A by permuting the + # the indices in the loops below. Only worth doing that + # if copying in or out of A is expensive. + + # permutate matrix based on ordering p + # W is in the order that the cliques are based on + W = A[p, p] + + # go through supernode tree in descending order (given a post-ordering). + # This is ensured in the get_snode, get_separators functions + + for j = (sntree.n_cliques - 1):-1:1 + + # in order to obtain ν, α the vertex numbers of the supernode are + # mapped to the new position of the permuted matrix index + # set of snd(i) sorted using the numerical ordering i,i+1,...i+ni + ν = get_snode(sntree, j) + + # index set containing the elements of col(i) \ snd(i) + #sorted using numerical ordering σ(i) + α = get_separators(sntree, j) + + # index set containing the row indices of the lower-triangular zeros in + # column i (i: representative index) sorted by σ(i) + i = ν[1] + + ν = collect(ν) #unborks indexing below + α = collect(α) #unborks indexing below + + η = collect(i+1:1:N) + filter!(x -> !insorted(x, α) && !insorted(x, ν), η) + + Waa = W[α, α] + Wαν = view(W, α, ν) + Wηα = view(W, η, α) + + Y = zeros(length(α), length(ν)) + try + Y[:, :] = Waa \ Wαν + catch + Waa_pinv = pinv(Waa) + Y[:, :] = Waa_pinv * Wαν + end + + W[η, ν] = Wηα * Y + + # symmetry condition + W[ν, η] = view(W, η, ν)' + + end + + # invert the permutation + A[:, :] = W[ip, ip] + +end + diff --git a/src/chordal/decomposition/reverse_compact.jl b/src/chordal/decomposition/reverse_compact.jl new file mode 100644 index 00000000..bf107955 --- /dev/null +++ b/src/chordal/decomposition/reverse_compact.jl @@ -0,0 +1,120 @@ +#AbstractVariables here should be DefaultVariables only, but this +# is not enforced since it causes a circular inclusion issue. + +# ----------------------------------- +# reverse the compact decomposition +# ----------------------------------- + +function decomp_reverse_compact!( + chordal_info::ChordalInfo{T}, + new_vars::AbstractVariables{T}, + old_vars::AbstractVariables{T}, + old_cones::Vector{SupportedCone} +) where {T} + + old_s = old_vars.s + old_z = old_vars.z + new_s = new_vars.s + new_z = new_vars.z + + # the cones for the originating problem, i.e. the cones + # that are compatible with the new_vars we want to populate, + # are held in chordal_info.init_cones + + cone_maps = chordal_info.cone_maps + row_ranges = collect(rng_cones_iterator(chordal_info.init_cones)) + row_ptr = 1 + + # add the blocks for each cone requires a buffer in which + # to hold sorted cliques. Allocate it here to avoid + # repeated allocations. The size of each clique is + # never bigger than the associated nblk + clique_buffer = zeros(Int,largest_nblk(chordal_info)) + + for (cone, cone_map) in zip(old_cones,cone_maps) + + row_range = row_ranges[cone_map.orig_index] + + if isnothing(cone_map.tree_and_clique) + row_ptr = add_blocks_with_cone!( + new_s, old_s, new_z, old_z, row_range, cone, row_ptr) + + else + @assert isa(cone, PSDTriangleConeT) + @assert !isnothing(cone_map.tree_and_clique) + + (tree_index, clique_index) = cone_map.tree_and_clique + pattern = chordal_info.spatterns[tree_index] + + row_ptr = add_blocks_with_sparsity_pattern!( + new_s, old_s, new_z, old_z, row_range, pattern, + clique_index, clique_buffer, row_ptr) + end + end + +end + +function add_blocks_with_sparsity_pattern!( + new_s::Vector{T}, + old_s::Vector{T}, + new_z::Vector{T}, + old_z::Vector{T}, + row_range::UnitRange{Int}, + spattern::SparsityPattern, + clique_index::Int, + clique_buffer::Vector{Int}, + row_ptr::Int +) where {T} + + sntree = spattern.sntree + ordering = spattern.ordering + + # load the clique into the buffer provided + clique = get_clique(sntree, clique_index) + resize!(clique_buffer, length(clique)) + for (i, v) = enumerate(clique) + clique_buffer[i] = ordering[v] + end + sort!(clique_buffer) + + counter = 0 + for j in clique_buffer, i in clique_buffer + if i <= j + offset = coord_to_upper_triangular_index((i, j)) - 1 + @views new_s[row_range.start + offset] += old_s[row_ptr + counter] + # notice: z overwrites (instead of adding) the overlapping entries + @views new_z[row_range.start + offset] = old_z[row_ptr + counter] + counter += 1 + end + end + + row_ptr + triangular_number(length(clique)) +end + +function add_blocks_with_cone!( + new_s::Vector{T}, + old_s::Vector{T}, + new_z::Vector{T}, + old_z::Vector{T}, + row_range::UnitRange{Int}, + cone::SupportedCone, + row_ptr::Int +) where {T} + + src_range = row_ptr:(row_ptr + nvars(cone) - 1) + @views new_s[row_range] .= old_s[src_range] + @views new_z[row_range] .= old_z[src_range] + row_ptr += nvars(cone) + +end + + +# the largest nblk across all spatterns + +function largest_nblk(chordal_info::ChordalInfo{T}) where {T} + max_block = 0 + for sp in chordal_info.spatterns + max_block = max(max_block, maximum(sp.sntree.nblk)) + end + return max_block +end \ No newline at end of file diff --git a/src/chordal/decomposition/reverse_standard.jl b/src/chordal/decomposition/reverse_standard.jl new file mode 100644 index 00000000..e77d2c0f --- /dev/null +++ b/src/chordal/decomposition/reverse_standard.jl @@ -0,0 +1,39 @@ +#----------------------------------- +# reverse the standard decomposition +#----------------------------------- + +function decomp_reverse_standard!( + chordal_info::ChordalInfo{T}, + new_vars::AbstractVariables{T}, + old_vars::AbstractVariables{T}, + _old_cones::Vector{SupportedCone} +) where {T} + + H = chordal_info.H + (_,m) = variables_dims(new_vars) + + mul!(new_vars.s, H, @view old_vars.s[(1+m):end]) + mul!(new_vars.z, H, @view old_vars.z[(1+m):end]) + + # to remove the overlaps we take the average of the values for + # each overlap by dividing by the number of blocks that overlap + #in a particular entry, i.e. number of 1s in each row of H + + rows, nnzs = number_of_overlaps_in_rows(H) + + for (ri, nnz) in zip(rows,nnzs) + new_vars.z[ri] /= nnz + end + +end + +function number_of_overlaps_in_rows( + A::SparseMatrixCSC{T} +) where {T} + + # sum the entries row-wise + n_overlaps = sum(A, dims = 2) + ri = findall(x -> x > 1, n_overlaps) + return ri, n_overlaps[ri] + + end \ No newline at end of file diff --git a/src/chordal/include.jl b/src/chordal/include.jl new file mode 100644 index 00000000..ce9f64a3 --- /dev/null +++ b/src/chordal/include.jl @@ -0,0 +1,17 @@ +include("./types.jl") +include("./supernode_tree.jl") + +include("./merge_strategy/defaults.jl") +include("./merge_strategy/nomerge.jl") +include("./merge_strategy/parent_child.jl") +include("./merge_strategy/clique_graph.jl") + +include("./sparsity_pattern.jl") +include("./chordal_info.jl") + +include("./decomposition/decomp.jl") +include("./decomposition/augment_standard.jl") +include("./decomposition/augment_compact.jl") +include("./decomposition/reverse_standard.jl") +include("./decomposition/reverse_compact.jl") +include("./decomposition/psd_completion.jl") \ No newline at end of file diff --git a/src/chordal/merge_strategy/clique_graph.jl b/src/chordal/merge_strategy/clique_graph.jl new file mode 100644 index 00000000..83fbdb4a --- /dev/null +++ b/src/chordal/merge_strategy/clique_graph.jl @@ -0,0 +1,735 @@ +# The (default) merge strategy based on the *reduced* clique graph ``\\mathcal{G}(\\mathcal{B}, \\xi)``, +# for a set of cliques ``\\mathcal{B} = \\{ \\mathcal{C}_1, \\dots, \\mathcal{C}_p\\}``, where the edge +# set ``\\xi`` is obtained by taking the edges of the union of clique trees. + +# Moreover, given an edge weighting function ``e(\\mathcal{C}_i,\\mathcal{C}_j) = w_{ij}``, we compute a +# weight for each edge that quantifies the computational savings of merging the two cliques. +# +# After the initial weights are computed, we merge cliques in a loop: +# +# **while** clique graph contains positive weights: +# - select two permissible cliques with the highest weight ``w_{ij}`` +# - merge cliques ``\\rightarrow`` update clique graph +# - recompute weights for updated clique graph +# +# See also: *Garstka, Cannon, Goulart - A clique graph based merging strategy for decomposable SDPs (2019)* +# +# NB: edges is currently an integer valued matrix since the weights are taken as +# powers of the cardinality of the intersection of the cliques. This needs to change +# to floats if empirical edge weight functions are to be supported. + +mutable struct CliqueGraphMergeStrategy <: AbstractMergeStrategy + stop::Bool # a flag to indicate that merging should be stopped + edges::SparseMatrixCSC{Int,Int} # the edges and weights of the reduced clique graph + p::Vector{Int} # as a workspace variable to store the sorting of weights + adjacency_table::Dict{Int, VertexSet} # a double structure of edges, to allow fast lookup of neighbors + edge_weight::EdgeWeightMethod # used to dispatch onto the correct scoring function + + function CliqueGraphMergeStrategy(; edge_weight = CUBIC::EdgeWeightMethod) + new(false, spzeros(Int, 0, 0), Int[], Dict{Int, VertexSet}(), edge_weight) + end +end + + +function initialise!(strategy::CliqueGraphMergeStrategy, t::SuperNodeTree) + + # this merge strategy is clique-graph based, we give up the tree structure and add + # the seperators to the supernodes. The supernodes then represent the full clique. + # after clique merging a new clique tree will be computed in post_process_merge! + # for this type + + for (snode,separator) in zip(t.snode, t.separators) + union!(snode, separator) + end + + for i in eachindex(t.snode_parent) + t.snode_parent[i] = INACTIVE_NODE + t.snode_children[i] = VertexSet() + end + + # compute the edges and intersections of cliques in the reduced clique graph + rows, cols = compute_reduced_clique_graph!(t.separators, t.snode) + + weights = compute_weights!(rows, cols, t.snode, strategy.edge_weight) + + strategy.edges = sparse(rows, cols, weights, t.n_cliques, t.n_cliques) + strategy.p = zeros(Int, length(strategy.edges.nzval)) + strategy.adjacency_table = compute_adjacency_table(strategy.edges, t.n_cliques) + return nothing + +end + +function is_done(strategy::CliqueGraphMergeStrategy) + strategy.stop +end + + +# Find the next two cliques in the clique graph `t` to merge. +function traverse(strategy::CliqueGraphMergeStrategy, t::SuperNodeTree) + + p = strategy.p + # find edge with highest weight, if permissible return cliques + edge = max_elem(strategy.edges) + + if ispermissible(edge, strategy.adjacency_table, t.snode) + return edge + end + + # sort the weights in edges.nzval to find the permutation p + sortperm!(view(p, 1:length(strategy.edges.nzval)), strategy.edges.nzval, alg = QuickSort, rev = true) + + # try edges with decreasing weight and check if the edge is permissible + for k = 2:length(strategy.edges.nzval) + edge = edge_from_index(strategy.edges, p[k]) + + if ispermissible(edge, strategy.adjacency_table, t.snode) + return (edge[1], edge[2]) + end + end + + # PJG: it seems possible to exit here and return "nothing". Handled + # as an option in Rust. Is it actually possible to reach this point? + # Maybe try debugging by using a crazy scoring function that is always + # terrible value + +end + + +function evaluate(strategy::CliqueGraphMergeStrategy, _t::SuperNodeTree, cand::Tuple{Int, Int}) + + (c1,c2) = cand + + do_merge = (strategy.edges[c1, c2] >= 0) + + if !do_merge + strategy.stop = true + end + return do_merge +end + + +function merge_two_cliques!(strategy::CliqueGraphMergeStrategy, t::SuperNodeTree, cand::Tuple{Int, Int}) + + (c1,c2) = cand + + # merge clique c2 into c1 + union!(t.snode[c1], t.snode[c2]) + empty!(t.snode[c2]) + + # decrement number of mergeable / nonempty cliques in graph + t.n_cliques -= 1 + + return nothing +end + + +# After a merge happened, update the reduced clique graph. + +function update_strategy!( + strategy::CliqueGraphMergeStrategy, + t::SuperNodeTree, + cand::Tuple{Int, Int}, + do_merge::Bool +) + do_merge || return; + + # After a merge operation update the information of the strategy + c_1_ind, c_removed = cand + + edges = strategy.edges + n = size(edges, 2) + adjacency_table = strategy.adjacency_table + + c_1 = t.snode[c_1_ind] + neighbors = adjacency_table[c_1_ind] + # neighbors exclusive to the removed clique (and not c1) + new_neighbors = setdiff(adjacency_table[c_removed], neighbors, c_1_ind) + + # recalculate edge values of all of c_1's neighbors + for n_ind in neighbors + if n_ind != c_removed + neighbor = t.snode[n_ind] + row = max(c_1_ind, n_ind) + col = min(c_1_ind, n_ind) + val = edge_metric(c_1, neighbor, strategy.edge_weight) + edges[row, col] = val + end + end + + # point edges exclusive to removed clique to surviving clique 1 + for n_ind in new_neighbors + neighbor = t.snode[n_ind] + row = max(c_1_ind, n_ind) + col = min(c_1_ind, n_ind) + val = edge_metric(c_1, neighbor, strategy.edge_weight) + edges[row, col] = val + end + + # overwrite the weight to any removed edges that still contain a link to c_removed + edges[c_removed+1:n, c_removed] .= 0 + edges[c_removed, 1:c_removed] .= 0 + + dropzeros!(edges) + + # update adjacency table in a similar manner + union!(adjacency_table[c_1_ind], new_neighbors) + for new_neighbor in new_neighbors + push!(adjacency_table[new_neighbor], c_1_ind) + end + delete!(adjacency_table, c_removed) + for (_, set) in adjacency_table + delete!(set, c_removed) + end + + return nothing +end + + +function post_process_merge!(strategy::CliqueGraphMergeStrategy, t::SuperNodeTree) + + # since for now we have a graph, not a tree, a post ordering or a parent structure + #does not make sense. Therefore just number the non-empty supernodes in t.snd + + t.snode_post = findall(x -> !isempty(x), t.snode) + t.snode_parent = fill(INACTIVE_NODE, length(t.snode)) + + # recompute a clique tree from the clique graph + t.n_cliques > 1 && clique_tree_from_graph!(strategy, t) + + # PJG: This seems unnecessary because the next operation on this + # object is the call to reorder_snode_consecutively, which overwrites + # the snode anyway. Treatment of separators possibly ends up different. + # Seems to work without, but keep for now for consistency with COSMO. + + foreach(s->sort!(s), t.snode) + foreach(s->sort!(s), t.separators) + + return nothing +end + + +# Given the cliques and edges of a clique graph, this function computes a valid clique tree. +# This is necessary to perform the psd completion step of the dual variable after solving the problem. + +function clique_tree_from_graph!(strategy::CliqueGraphMergeStrategy, t::SuperNodeTree) + + # a clique tree is a maximum weight spanning tree of the clique graph, where the edge weight is the + # cardinality of the intersection between two cliques compute intersection value for each edge + # in the clique graph + clique_intersections!(strategy.edges, t.snode) + + # find a maximum weight spanning tree of the clique graph using Kruskal's algorithm + kruskal!(strategy.edges, t.n_cliques) + + # determine the root clique of the clique tree (it can be any clique, but we use the + # clique that contains the vertex with the highest order) + determine_parent_cliques!(t.snode_parent, t.snode_children, t.snode, t.post, strategy.edges) + + # recompute a postorder for the supernodes (NB: snode_post will shrink + # to the possibly reduced length n_cliques after the merge) + post_order!(t.snode_post, t.snode_parent, t.snode_children, t.n_cliques) + + # Clear the (graph) separators. They will be rebuilt in the split_cliques + map(set -> empty!(set), t.separators) + + # split clique sets back into separators and supernodes + split_cliques!(t.snode, t.separators, t.snode_parent, t.snode_post, t.n_cliques) + + return nothing + +end + + + + +#------------------- internal utilities -------------------# + +# Compute the reduced clique graph (union of all clique trees) given an initial clique tree defined by its +# supernodes and separator sets. + +# We are using the algorithm described in **Michel Habib and Juraj Stacho - Polynomial-time algorithm for the +# leafage of chordal graphs (2009)**, which +# computes the reduced clique graph in the following way: +# 1. Sort all minimal separators by size +# 2. Initialise graph CG(R) with cliques as nodes and no edges +# 3. for largest unprocessed separator S and +# | add an edge between any two cliques C1 and C2 if they both contain S and are in different connected +# components of CG(R) and store in `edges`. +# | Compute an edge weight used for merge decision and store in `val`. +# | Store the index of the separator which is the intersection C1 ∩ C2 in `iter` +# end + +function compute_reduced_clique_graph!(separators::Vector{VertexSet}, snode::Vector{VertexSet}) + + # loop over separators by decreasing cardinality + sort!(separators, by = x -> length(x), rev = true) + + rows = Int[] + cols = Int[] + + for separator in separators + # find cliques that contain the separator + clique_indices = findall(x -> separator ⊆ x, snode) + + # Compute the separator graph (see Habib, Stacho - Reduced clique graphs of chordal graphs) + # to analyse connectivity. We represent the separator graph H by a hashtable + H = separator_graph(clique_indices, separator, snode) + # find the connected components of H + components = find_components(H, clique_indices) + + # for each pair of cliques that contain the separator, add an edge to the reduced + # clique tree if they are in unconnected components + + ncliques = length(clique_indices) + for i in 1:ncliques, j in (i+1):ncliques + pair = (clique_indices[i], clique_indices[j]) + if is_unconnected(pair, components) + push!(rows, max(pair...)) + push!(cols, min(pair...)) + end + end + end + + return rows, cols + +end + +# Find the separator graph H given a separator and the relevant index-subset of cliques. + +function separator_graph(clique_ind::Vector{Int}, separator::VertexSet, snd::Vector{VertexSet}) + + # make the separator graph using a hash table + # key: clique_ind --> edges to other clique indices + H = Dict{Int, Array{Int, 1}}() + + nindex = length(clique_ind) + for i in 1:nindex, j in (i+1):nindex + ca = clique_ind[i] + cb = clique_ind[j] + # if intersect_dim(snd[ca], snd[cb]) > length(separator) + if !inter_equal(snd[ca], snd[cb], separator) + if haskey(H, ca) + push!(H[ca], cb) + else + H[ca] = [cb] + end + if haskey(H, cb) + push!(H[cb], ca) + else + H[cb] = [ca] + end + end + end + # add unconnected cliques + for v in clique_ind + !haskey(H, v) && (H[v] = Int[]) + end + return H +end + + +# Find connected components in undirected separator graph represented by `H`. +function find_components(H::Dict{Int, Vector{Int}}, clique_ind::Vector{Int}) + visited = Dict{Int, Bool}(v => false for v in clique_ind) + components = Vector{VertexSet}() + for v in clique_ind + if visited[v] == false + component = VertexSet() + DFS_hashtable!(component, v, visited, H) + push!(components, component) + end + end + return components +end + + +# Check whether the `pair` of cliques are in different `components`. +function is_unconnected(pair::Tuple{Int, Int}, components::Vector{VertexSet}) + component_ind = findfirst(x -> pair[1] ∈ x, components) + return pair[2] ∉ components[component_ind] +end + + +# Depth first search on a hashtable `H`. +function DFS_hashtable!( + component::VertexSet, + v::Int, + visited::Dict{Int, Bool}, + H::Dict{Int, Vector{Int}} +) + + visited[v] = true + push!(component, v) + for n in H[v] + if visited[n] == false + DFS_hashtable!(component, n, visited, H) + end + end +end + + +# Check if s ∩ s2 == s3. +function inter_equal(s1::T, s2::T, s3::T) where {T <: AbstractSet} + dim = 0 + + len_s1 = length(s1) + len_s2 = length(s2) + len_s3 = length(s3) + + # maximum possible intersection size + max_intersect = len_s1 + len_s2 + + # abort if there's no way the intersection can be the same + if max_intersect < len_s3 + return false + end + + if len_s1 < len_s2 + sa = s1 + sb = s2 + else + sa = s2 + sb = s1 + end + for e in sa + if e in sb + dim += 1 + dim > len_s3 && return false + e in s3 || return false + end + max_intersect -= 1 + max_intersect < len_s3 && return false + end + return dim == len_s3 +end + + +# Given a list of edges, return an adjacency hash-table `table` with nodes from 1 to `num_vertices`. + +function compute_adjacency_table(edges::SparseMatrixCSC{T}, num_vertices::Int) where {T} + + table = Dict(i => VertexSet() for i = 1:num_vertices) + r = edges.rowval + c = edges.colptr + for col = 1:num_vertices + for i in c[col]:c[col+1]-1 + row = r[i] + push!(table[row], col) + push!(table[col], row) + end + end + return table +end + +# Check whether `edge` is permissible for a merge. An edge is permissible if for every common neighbor N, +# C_1 ∩ N == C_2 ∩ N or if no common neighbors exist. + +function ispermissible( + edge::Tuple{Integer, Integer}, + adjacency_table::Dict{Int, VertexSet}, + snode::Vector{VertexSet} +) + + c_1 = edge[1] + c_2 = edge[2] + common_neighbors = intersect(adjacency_table[c_1], adjacency_table[c_2]) + + # N.B. This is allocating and could be made more efficient + # Here OrderedSets return equality if they have the same + # elements, regardless of insertion order. + for neighbor in common_neighbors + intersect(snode[c_1], snode[neighbor]) != intersect(snode[c_2], snode[neighbor]) && return false + end + return true +end + + +# Find the matrix indices (i, j) of the first maximum element among the elements stored in A.nzval + +function max_elem(A::SparseMatrixCSC{T}) where {T} + length(A.nzval) == 0 && throw(DomainError("Sparse matrix A doesn't contain any entries")) + n = size(A, 2) + + ~, ind = findmax(A.nzval) + row = A.rowval[ind] + + col = 0 + for c = 1:n + col_indices = A.colptr[c]:A.colptr[c+1]-1 + if in(ind, col_indices) + col = c + break; + end + end + return (row, col) +end + +function edge_from_index(A::SparseMatrixCSC{T, Int}, ind::Int) where {T} + index_to_coord(A,ind) + end + + + +function clique_intersections!(E::SparseMatrixCSC{T}, snd::Vector{VertexSet}) where {T} + # iterate over the nonzeros of the connectivity matrix E which represents the + # clique graph and replace the value by |C_i ∩ C_j| + rows = rowvals(E) + for col in 1:size(E, 2) + for j in nzrange(E, col) + row = rows[j] + E.nzval[j] = intersect_dim(snd[row], snd[col]) + end + end + return nothing + end + +# Return the number of elements in s ∩ s2. +function intersect_dim(s1::AbstractSet, s2::AbstractSet) + if length(s1) < length(s2) + sa = s1 + sb = s2 + else + sa = s2 + sb = s1 + end + dim = 0 + for e in sa + e in sb && (dim += 1) + end + return dim +end + +# Find the size of the set `A ∪ B` under the assumption that `A` and `B` only have unique elements. +function union_dim(s1::T, s2::T) where {T <: AbstractSet} + + length(s1) + length(s2) - intersect_dim(s1, s2) + +end + + +# Kruskal's algorithm to find a maximum weight spanning tree from the clique intersection graph. +# +# `E[i,j]` holds the cardinalities of the intersection between two cliques (i, j). Changes the entries in the +# connectivity matrix `E` to a negative balue if an edge between two cliques is included in the max spanning tree. +# +# This is a modified version of https://github.com/JuliaGraphs/LightGraphs.jl/blob/master/src/spanningtrees/kruskal.jl + + +function kruskal!(E::SparseMatrixCSC{T}, num_cliques::Int) where{T} + + num_initial_cliques = size(E, 2) + + connected_c = DataStructures.IntDisjointSets(num_initial_cliques) + + I, J, V = findnz(E) + + # sort the weights and edges from maximum to minimum value + p = sortperm(V, rev = true) + I = I[p] + J = J[p] + num_edges_found = 0 + + # iterate through edges (I -- J) with decreasing weight + for k in eachindex(I) + row = I[k] + col = J[k] + + if !in_same_set(connected_c, row, col) + union!(connected_c, row, col) + # indicate an edge in the MST with a negative value in E (all other values are >= 0) + # PJG: I think here E[row,col] == E.nzval[p[k]] and could be assigned that way + E[row, col] = -1.0 + num_edges_found += 1 + # break when all cliques are connected in one tree + num_edges_found >= num_cliques - 1 && break + end + end + return nothing +end + +# Given the maximum weight spanning tree represented by `E`, determine a parent +# structure `snd_par` for the clique tree. + +function determine_parent_cliques!( + snode_parent::Vector{Int}, + snode_children::Vector{VertexSet}, + cliques::Vector{VertexSet}, + post::Vector{Int}, + E::SparseMatrixCSC{T} + ) where {T} + + # vertex with highest order + v = post[end] + c = 0 + # Find clique that contains that vertex + for (k, clique) in enumerate(cliques) + if v ∈ clique + # set that clique to the root + snode_parent[k] = NO_PARENT + c = k + break + end + end + + # assign children to cliques along the MST defined by E + assign_children!(snode_parent, snode_children, c, E) + + return nothing +end + + +# function assign_children!( +# snode_parent::Vector{Int}, +# snode_children::Vector{VertexSet}, +# c::Int, +# edges::SparseMatrixCSC{T} +# ) where {T} + +# # determine neighbors +# neighbors = find_neighbors(edges, c) + +# for n in neighbors +# # conditions that there is a edge in the MST and that n is not the parent of c +# if edges[max(c, n), min(c, n)] == -1.0 && snode_parent[c] != n +# snode_parent[n] = c +# push!(snode_children[c], n) +# assign_children!(snode_parent, snode_children, n, edges) +# end +# end +# return nothing +# end + +#NB: non-recursive version of the COSMO.jl implementation +function assign_children!( + snode_parent::Vector{Int}, + snode_children::Vector{VertexSet}, + c::Int, + edges::SparseMatrixCSC{T} +) where {T} + + stack = [c] + while !isempty(stack) + c = pop!(stack) + neighbors = find_neighbors(edges, c) + + for n in neighbors + # conditions that there is an edge in the MST and that n is not the parent of c + if edges[max(c, n), min(c, n)] == -1.0 && snode_parent[c] != n + snode_parent[n] = c + push!(snode_children[c], n) + push!(stack, n) # Add child to stack for processing + end + end + end +end + + +# Find all the cliques connected to `c` which are given by the nonzeros in `(c, 1:c-1)` and `(c+1:n, c)`. + +function find_neighbors(edges::SparseMatrixCSC, c::Int) + neighbors = zeros(Int, 0) + _, n = size(edges) + # find all nonzero columns in row c up to column c + if c > 1 + neighbors = vcat(neighbors, findall(x -> x != 0, edges[c, 1:c-1])) + end + # find all nonzero entries in column c below c + if c < n + rows = @view edges.rowval[edges.colptr[c]:edges.colptr[c+1]-1] + if edges.colptr[c] <= edges.colptr[c+1] - 1 + neighbors = vcat(neighbors, rows) + end + end + return neighbors +end + + +# Traverse the clique tree in descending topological order and split the clique sets into supernodes and separators. + +function split_cliques!( + snode::Vector{VertexSet}, + separators::Vector{VertexSet}, + snode_parent::Vector{Int}, + snode_post::Vector{Int}, + num_cliques::Int +) + + # travese in topological decending order through the clique tree and split the clique + # into supernodes and separators + for j = 1:1:(num_cliques - 1) + c_ind = snode_post[j] + p_ind = snode_parent[c_ind] + + # find intersection of clique with parent + separators[c_ind] = intersect(snode[c_ind], snode[p_ind]) + snode[c_ind] = filter!(x -> x ∉ separators[c_ind], snode[c_ind]) + end + return nothing +end + + +# ------------------- +# functions relating to edge weights +# ------------------- + + +# Compute the edge weight between all cliques specified by the edges (rows, cols). +# weights on the edges currently defined as integer values, but could be changed +# to floats to allow emperical edge weight functions. + +function compute_weights!( + rows::Vector{Int}, + cols::Vector{Int}, + snode::Vector{VertexSet}, + edge_weight::EdgeWeightMethod +) + + weights = zeros(Int, length(rows)) + + for k = 1:length(rows) + c_1 = snode[rows[k]] + c_2 = snode[cols[k]] + weights[k] = edge_metric(c_1, c_2, edge_weight) + end + return weights + +end + + + +#Given two cliques `c_a` and `c_b` return a value for their edge weight. + +function _edge_metric(c_a::VertexSet, c_b::VertexSet, edge_weight::EdgeWeightMethod) + n_1 = length(c_a) + n_2 = length(c_b) + + # merged block size + n_m = union_dim(c_a, c_b) + + if edge_weight == CUBIC::EdgeWeightMethod + return (n_1^3 + n_2^3 - n_m^3)::Int + else + throw(ArgumentError("Unknown weight metric not implemented")) + end +end + + +#PJG: this function appears to give better performance, but can be deactived +#if favor of the above to agree with COSMO results for testing + +function edge_metric(c_a::VertexSet, c_b::VertexSet, edge_weight::EdgeWeightMethod) + + n_1 = triangular_number(length(c_a)) + n_2 = triangular_number(length(c_b)) + + # merged block size + n_m = triangular_number(union_dim(c_a, c_b)) + + if edge_weight == CUBIC::EdgeWeightMethod + return (n_1^3 + n_2^3 - n_m^3)::Int + else + throw(ArgumentError("Unknown weight metric not implemented")) + end +end + + diff --git a/src/chordal/merge_strategy/defaults.jl b/src/chordal/merge_strategy/defaults.jl new file mode 100644 index 00000000..b820b4a8 --- /dev/null +++ b/src/chordal/merge_strategy/defaults.jl @@ -0,0 +1,72 @@ +# Main clique merging routine, common to all strategies + +function merge_cliques!(strategy::AbstractMergeStrategy, t::SuperNodeTree) + + initialise!(strategy, t) + + while !is_done(strategy) + + # find merge candidates + cand = traverse(strategy, t) + + # evaluate wether to merge the candidates + do_merge = evaluate(strategy, t, cand) + if do_merge + merge_two_cliques!(strategy, t, cand) + end + + # update strategy information after the merge + update_strategy!(strategy, t, cand, do_merge) + t.n_cliques == 1 && break + + end + + post_process_merge!(strategy, t) + return nothing +end + +#All strategies must then implement the following interface: + +# initialise! - tree and strategy +# is_done - merging complete, so can stop the merging process +# traverse - find the next merge candidates +# evalute - evaluate wether to merge them or not +# merge_two_cliques! - execute a merge +# update_strategy! - update the tree/graph and strategy +# post_process_merge! - do any post-processing of the tree/graph + +function is_done(strategy::AbstractMergeStrategy) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return nothing +end + +function initialise!(strategy::AbstractMergeStrategy, t::SuperNodeTree) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return nothing +end + +function traverse(strategy::AbstractMergeStrategy, t::SuperNodeTree) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return nothing +end + +function evaluate(strategy::AbstractMergeStrategy, t::SuperNodeTree, cand::Tuple{Int, Int}) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return bool +end + +function merge_two_cliques!(strategy::AbstractMergeStrategy, t::SuperNodeTree, cand::Tuple{Int, Int}) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return nothing +end + +function update_strategy!(strategy::AbstractMergeStrategy, t::SuperNodeTree, cand::Tuple{Int, Int}, do_merge::Bool) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return nothing +end + +function post_process_merge!(strategy::AbstractMergeStrategy, t::SuperNodeTree) + error("Incomplete merge strategy specification: ",typeof(strategy)) + #return nothing +end + diff --git a/src/chordal/merge_strategy/nomerge.jl b/src/chordal/merge_strategy/nomerge.jl new file mode 100644 index 00000000..f5fcbaea --- /dev/null +++ b/src/chordal/merge_strategy/nomerge.jl @@ -0,0 +1,17 @@ +struct NoMergeStrategy <: AbstractMergeStrategy end + +function initialise!(strategy::NoMergeStrategy, t::SuperNodeTree) + return nothing +end + +function is_done(strategy::NoMergeStrategy) + true +end + +function post_process_merge!(strategy::NoMergeStrategy, t::SuperNodeTree) + return nothing +end + + +# All other functions should be unreachable this is not that nice, +# and should revert the the default / error implementation if called \ No newline at end of file diff --git a/src/chordal/merge_strategy/parent_child.jl b/src/chordal/merge_strategy/parent_child.jl new file mode 100644 index 00000000..74685651 --- /dev/null +++ b/src/chordal/merge_strategy/parent_child.jl @@ -0,0 +1,150 @@ +mutable struct ParentChildMergeStrategy <: AbstractMergeStrategy + stop::Bool + clique_index::Int + t_fill::Int + t_size::Int + + # PJG: fill and size need to be settable + function ParentChildMergeStrategy(; t_fill = 8, t_size = 8) + new(false, 0, t_fill, t_size) + end +end + +function initialise!(strategy::ParentChildMergeStrategy, t::SuperNodeTree) + # start with node that has second highest order + strategy.clique_index = length(t.snode) - 1 +end + +function is_done(strategy::ParentChildMergeStrategy) + strategy.stop +end + +# Traverse tree `t` in descending topological order and return parent and clique (root has highest order). + +function traverse(strategy::ParentChildMergeStrategy, t::SuperNodeTree) + + c = t.snode_post[strategy.clique_index] + + return (t.snode_parent[c], c) +end + + + +# Decide whether to merge the two clique candidates. + +function evaluate(strategy::ParentChildMergeStrategy, t::SuperNodeTree, cand::Tuple{Int, Int}) + + strategy.stop && return false + + (parent, child) = cand + + dim_parent_snode, dim_parent_sep = clique_dim(t, parent) + dim_clique_snode, dim_clique_sep = clique_dim(t, child) + + fill = fill_in(dim_clique_snode, dim_clique_sep, dim_parent_snode, dim_parent_sep) + max_snode = max(dim_clique_snode, dim_parent_snode) + + return fill <= strategy.t_fill || max_snode <= strategy.t_size +end + +# Given the clique tree `t` merge the two cliques with indices in `cand` as parent-child. + +function merge_two_cliques!( + strategy::ParentChildMergeStrategy, + t::SuperNodeTree, + cand::Tuple{Int, Int} +) + + # determine which clique is the parent + (p, ch) = determine_parent(t, cand...); + + # merge child's vertex sets into parent's vertex set + union!(t.snode[p], t.snode[ch]) + empty!(t.snode[ch]) + empty!(t.separators[ch]) + + # update parent structure + for grandch in t.snode_children[ch] + t.snode_parent[grandch] = p + end + t.snode_parent[ch] = INACTIVE_NODE + + # update children structure + delete!(t.snode_children[p], ch) + union!(t.snode_children[p], t.snode_children[ch]) + empty!(t.snode_children[ch]) + + # decrement number of cliques in tree + t.n_cliques -= 1 + return nothing + +end + + +# After a merge attempt, update the strategy information. +function update_strategy!( + strategy::ParentChildMergeStrategy, + t::SuperNodeTree, + cand::Tuple{Int, Int}, + do_merge::Bool +) + # try to merge last node of order 1, then stop + if strategy.clique_index == 1 + strategy.stop = true + # otherwise decrement node index + else + strategy.clique_index -= 1 + end +end + + + +function post_process_merge!(strategy::ParentChildMergeStrategy, t::SuperNodeTree) + + # the merging creates empty supernodes and seperators, recalculate a + # post order for the supernodes (shrinks to length t.n_cliques ) + post_order!(t.snode_post, t.snode_parent, t.snode_children, t.n_cliques) + + return nothing +end + + + +#-------------------- utilities -------------------- + + +# Given two cliques `c1` and `c2` in the tree `t`, return the parent clique first. + +# Not implemented as part of the general SuperNodeTree interface +# since this should only be called when we can guarantee that we +# are acting on a parent-child pair. + +function determine_parent(t::SuperNodeTree, c1::Int, c2::Int) + if in(c2, t.snode_children[c1]) + return c1, c2 + else + return c2, c1 + end +end + +# not implemented as part of the main SuperNodeTree interface since the +# index is not passed through the post ordering +function clique_dim(t::SuperNodeTree, i::Int) + return length(t.snode[i]), length(t.separators[i]) +end + + +# Compute the amount of fill-in created by merging two cliques with the +# respective supernode and separator dimensions. + +function fill_in( + dim_clique_snode::Int, + dim_clique_sep::Int, + dim_parent_snode::Int, + dim_parent_sep::Int +) + dim_parent = dim_parent_snode + dim_parent_sep + dim_clique = dim_clique_snode + dim_clique_sep + + return (dim_parent - dim_clique_sep) * (dim_clique - dim_clique_sep) +end \ No newline at end of file diff --git a/src/chordal/sparsity_pattern.jl b/src/chordal/sparsity_pattern.jl new file mode 100644 index 00000000..031765cd --- /dev/null +++ b/src/chordal/sparsity_pattern.jl @@ -0,0 +1,49 @@ +# --------------------------- +# Struct to hold clique and sparsity data for a constraint +# --------------------------- +mutable struct SparsityPattern + sntree::SuperNodeTree + ordering::Array{Int} + orig_index::Int # original index of the cone being decomposed + + # constructor for sparsity pattern + function SparsityPattern( + L::SparseMatrixCSC{T}, + ordering::Array{Int, 1}, + orig_index::Int, + merge_method::Symbol, + + ) where {T} + + sntree = SuperNodeTree(L) + + # clique merging only if more than one clique present + + if sntree.n_cliques > 1 + + if merge_method == :none + merge_cliques!(NoMergeStrategy(), sntree) + + elseif merge_method == :parent_child + merge_cliques!(ParentChildMergeStrategy(), sntree) + + elseif merge_method == :clique_graph + merge_cliques!(CliqueGraphMergeStrategy(), sntree) + + else + error("Unknown merge strategy: ", merge_method) + end + + end + + # reorder vertices in supernodes to have consecutive order + # necessary for equal column structure for psd completion + reorder_snode_consecutively!(sntree, ordering) + + # for each clique determine the number of entries of the block + #represented by that clique + calculate_block_dimensions!(sntree) + + return new(sntree, ordering, orig_index) + end +end \ No newline at end of file diff --git a/src/chordal/supernode_tree.jl b/src/chordal/supernode_tree.jl new file mode 100644 index 00000000..adcf9f71 --- /dev/null +++ b/src/chordal/supernode_tree.jl @@ -0,0 +1,395 @@ + + +# this value used to mark root notes, i.e. ones with no parent +const NO_PARENT = typemax(Int); + +# when cliques are merged, their vertices are marked thusly +const INACTIVE_NODE = typemax(Int) - 1; + +# A structure to represent and analyse the sparsity pattern of an LDL factor matrix L. +mutable struct SuperNodeTree + + # vertices of supernodes stored in one array (also called residuals) + snode::Vector{VertexSet} + # post order of supernodal elimination tree + snode_post::Vector{Int} + # parent of each supernodes + snode_parent::Vector{Int} + # children of each supernode + snode_children::Vector{VertexSet} + # post ordering of the vertices in elim tree σ(j) = v + post::Array{Int} + # vertices of clique seperators + separators::Vector{VertexSet} + + # sizes of submatrices defined by each clique, sorted by post-ordering, + # e.g. size of clique with order 3 => nblk[3]. Only populated + # after a post-merging call to `calculate_block_dimensions!` + nblk::Option{Vector{Int}} + + # number of nonempty supernodes / cliques in tree + n_cliques::Int + + function SuperNodeTree(L::SparseMatrixCSC{T}) where {T} + + parent = parent_from_L(L) + children = children_from_parent(parent) + post = zeros(Int, length(parent)) + post_order!(post, parent, children, length(parent)) + + degree = higher_degree(L) + snode, snode_parent = find_supernodes(parent, post, degree) + + snode_children = children_from_parent(snode_parent) + snode_post = zeros(Int,length(snode_parent)) + post_order!(snode_post, snode_parent, snode_children, length(snode_parent)) + + # Here we find separators in all cases, unlike COSMO which defers until + # after merging for the clique graph merging case. These are later + # modified in the clique-graph merge case in the call to add_separators + separators = find_separators(L, snode) + + # nblk will be allocated to the length of the *post-merging* + # supernode count in calculate_block_dimensions! + nblk = nothing + + # number of cliques / nonempty supernodes. decrements as supernodes are merged + n_cliques = length(snode) + + new(snode, snode_post, snode_parent, snode_children, + post, separators, nblk, n_cliques) + + end + +end + +# ------------------------------------------ +# functions implemented for the SuperNodeTree +# ------------------------------------------ + +function get_post_order(sntree::SuperNodeTree, i::Int) + return sntree.snode_post[i] +end + +function get_snode(sntree::SuperNodeTree, i::Int) + return sntree.snode[sntree.snode_post[i]] +end + +function get_separators(sntree::SuperNodeTree, i::Int) + return sntree.separators[sntree.snode_post[i]] +end + +function get_clique_parent(sntree::SuperNodeTree, clique_index::Int) + return sntree.snode_parent[sntree.snode_post[clique_index]] +end + +# the block sizes are stored in post order, e.g. if clique 4 (stored in pos 4) +# has order 2, then nblk[2] represents the cardinality of clique 4 +function get_nblk(sntree::SuperNodeTree, i::Int) + return sntree.nblk[i] +end + +function get_overlap(sntree::SuperNodeTree, i::Int) + return length(sntree.separators[sntree.snode_post[i]]) +end + +function get_clique(sntree::SuperNodeTree, i::Int) + c = sntree.snode_post[i] + union(sntree.snode[c], sntree.separators[c]) +end + +function get_decomposed_dim_and_overlaps(sntree::SuperNodeTree) + dim = 0 + overlaps = 0 + for i = 1:sntree.n_cliques + dim += triangular_number(get_nblk(sntree, i)) + overlaps += triangular_number(get_overlap(sntree, i)) + end + (dim, overlaps) +end + +# Takes a SuperNodeTree and reorders the vertices in each supernode (and separator) to have consecutive order. + +# The reordering is needed to achieve equal column structure for the psd completion of the dual variable `Y`. +# This also modifies `ordering` which maps the vertices in the `sntree` back to the actual location in the +# not reordered data, i.e. the primal constraint variable `S` and dual variables `Y`. + + +function reorder_snode_consecutively!(t::SuperNodeTree, ordering::Vector{Int}) + + # determine permutation vector p and permute the vertices in each snd + p = zeros(Int,length(t.post)) + + k = 1 + for i in t.snode_post + + snode = t.snode[i] + n = length(snode) + viewp = view(p, k:k+n-1) + viewp .= t.snode[i] + sort!(viewp) + + #assign k:(k+n)-1 to the OrderedSet snode, + #dropping the previous values + empty!(snode) + foreach(v->push!(snode,v), k:(k+n)-1) + + k += n + end + + # permute the separators as well + p_inv = invperm(p) + for sp in t.separators + + # use here the permutation vector p as scratch before flushing + # the separator set and repopulating. Assumes that the permutation + # will be at least as long as the largest separator set + + @assert length(p) >= length(sp) + tmp = view(p,1:length(sp)) + + tmp .= Iterators.map(x -> p_inv[x], sp) + empty!(sp) + foreach(v->push!(sp,v), tmp) + end + + # because I used 'p' as scratch space, I will + # ipermute using pinv rather than permute using p + + invpermute!(ordering, p_inv) + return nothing +end + +# The the block dimension vector is the size of the remaining +# (unemptied) supernodes. Before this call, t.nblk should be +# nothing + +function calculate_block_dimensions!(t::SuperNodeTree) +n = t.n_cliques +t.nblk = zeros(Int,n) + + for i = 1:n + c = t.snode_post[i] + t.nblk[i] = length(t.separators[c]) + length(t.snode[c]) + end +end + + +# --------------------------- +# utility functions for SuperNodeTree + + +function parent_from_L(L::SparseMatrixCSC{T,Int}) where {T} + + parent = fill(NO_PARENT, L.n) + # loop over vertices of graph + for i = 1:L.n + parent[i] = find_parent_direct(L, i) + end + return parent +end + +function find_parent_direct(L::SparseMatrixCSC{T}, v::Int) where{T} + v == size(L, 1) && return NO_PARENT + return L.rowval[L.colptr[v]] +end + + +function find_separators( + L::SparseMatrixCSC{T,Int}, + snode::Vector{VertexSet} +) where {T} + separators = new_vertex_sets(length(snode)) + + for (sn, sep) in zip(snode, separators) + vrep = minimum(sn) + adjplus = find_higher_order_neighbors(L, vrep) + + for neighbor in adjplus + if neighbor ∉ sn + push!(sep, neighbor) + end + end + end + + return separators + +end + +function find_higher_order_neighbors(L::SparseMatrixCSC, v::Int) + col_ptr = L.colptr + row_val = L.rowval + return view(row_val, col_ptr[v]:col_ptr[v + 1] - 1) +end + +# findall the cardinality of adj+(v) for all v in V +function higher_degree(L::SparseMatrixCSC{T}) where{T} + + degree = zeros(Int, L.n) + for v = 1:(L.n-1) + degree[v] = L.colptr[v + 1] - L.colptr[v] + end + return degree +end + + +function children_from_parent(parent::Vector{Int}) + + children = new_vertex_sets(length(parent)) + for (i,pi) = enumerate(parent) + pi != NO_PARENT && push!(children[pi], i) + end + return children +end + + +# This could be faster for the case that merges happened, i.e. nc != length(parent) + +function post_order!(post::Vector{Int}, parent::Vector{Int}, children::Vector{VertexSet}, nc::Int) + + order = (nc + 1) * ones(Int, length(parent)) + + root = findfirst(x -> x == NO_PARENT, parent) + + stack = sizehint!(Int[], length(parent)) + push!(stack, root) + + resize!(post, length(parent)) + post .= 1:length(parent) + + i = nc + + while !isempty(stack) + v = pop!(stack) + order[v] = i + i -= 1 + + # maybe faster to append to the stack vector and then + # sort a view of what was added, but this way gets + # the children sorted and keeps everything consistent + # with the COSMO implementation for reference + append!(stack, sort!(children[v])) + end + + sort!(post, by = x-> order[x]) + + # if merges happened, remove the entries pointing to empty arrays / cliques + nc != length(parent) && resize!(post, nc) + +end + + + + + +function find_supernodes(parent::Vector{Int}, post::Vector{Int}, degree::Vector{Int}) + + snode = new_vertex_sets(length(parent)) + + snode_parent, snode_index = pothen_sun(parent, post, degree) + + for (i, f) in enumerate(snode_index) + if f < 0 + push!(snode[i], i) + else + push!(snode[f], i) + end + end + filter!(x -> !isempty(x), snode) + return snode, snode_parent + +end + + +# Algorithm from A. Poten and C. Sun: Compact Clique Tree Data Structures in Sparse Matrix Factorizations (1989) + +function pothen_sun(parent::Vector{Int}, post::Vector{Int}, degree::Vector{Int}) + + n = length(parent) + + # if snode_index[v] < 0 then v is a rep vertex, otherwise v ∈ supernode[snode_index[v]] + snode_index = fill(-one(Int), n) + snode_parent = fill(NO_PARENT, n) + + # This also works as array of Int[], which might be faster + # note this arrays is local to the function, not the one + # contained in the SuperNodeTree + children = new_vertex_sets(length(parent)) + + # find the root + root_index = findfirst(x -> x == NO_PARENT, parent) + + # go through parents of vertices in post_order + for v in post + + if parent[v] == NO_PARENT + push!(children[root_index], v) + else + push!(children[parent[v]], v) + end + + # parent is not the root. + if parent[v] != NO_PARENT + if degree[v] - 1 == degree[parent[v]] && snode_index[parent[v]] == -1 + # Case A: v is a representative vertex + if snode_index[v] < 0 + snode_index[parent[v]] = v + snode_index[v] -= 1 + # Case B: v is not representative vertex, add to sn_ind[v] instead + else + snode_index[parent[v]] = snode_index[v] + snode_index[snode_index[v]] -= 1 + end + else + if snode_index[v] < 0 + snode_parent[v] = v + else + snode_parent[snode_index[v]] = snode_index[v] + end + end + end + + # k: rep vertex of the snd that v belongs to + if snode_index[v] < 0 + k = v + else + k = snode_index[v] + end + # loop over v's children + v_children = children[v] + if !isempty(v_children) + for w in v_children + if snode_index[w] < 0 + l = w + else + l = snode_index[w] + end + if l != k + snode_parent[l] = k + end + end + end + end # loop over vertices + + repr_vertex = findall(x-> x < 0, snode_index) + # vertices that are the parent of representative vertices + repr_parent = snode_parent[repr_vertex] + + # resize and reset snode_parent to take into account that all + # non-representative arrays are removed from the parent structure + resize!(snode_parent, length(repr_vertex)) + snode_parent .= NO_PARENT + + for (i, rp) in enumerate(repr_parent) + rpidx = findfirst(x -> x == rp, repr_vertex) + isnothing(rpidx) && (rpidx = NO_PARENT) + snode_parent[i] = rpidx + end + + return snode_parent, snode_index +end + + +function new_vertex_sets(n::Int) + [VertexSet() for _ = 1:n] +end \ No newline at end of file diff --git a/src/chordal/types.jl b/src/chordal/types.jl new file mode 100644 index 00000000..c0dcb69a --- /dev/null +++ b/src/chordal/types.jl @@ -0,0 +1,15 @@ +# DataStructures provides OrderedSet, but only because it itself includes +# OrderedCollections. Could just use OrderedCollections directly, +# but the kruskal! function also uses a function from DataStructures. +# If that is removed then use the lighter weight option here + +using DataStructures + +abstract type AbstractMergeStrategy end +VertexSet = OrderedSet{Int} + +#PJG: make a settable option +@enum EdgeWeightMethod begin + CUBIC = 1 +end + diff --git a/src/cones/compositecone_type.jl b/src/cones/compositecone_type.jl index 87767712..6f04ecc4 100644 --- a/src/cones/compositecone_type.jl +++ b/src/cones/compositecone_type.jl @@ -15,18 +15,19 @@ struct CompositeCone{T} <: AbstractCone{T} numel::DefaultInt degree::DefaultInt - #a vector showing the overall index of the - #first element in each cone. For convenience - headidx::Vector{Int} + #ranges for the indices of the constituent cones + rng_cones::Vector{UnitRange{Int64}} + + #ranges for the indices of the constituent Hs blocks + #associated with each cone + rng_blocks::Vector{UnitRange{Int64}} # the flag for symmetric cone check _is_symmetric::Bool function CompositeCone{T}(cone_specs::Vector{SupportedCone}) where {T} - ncones = length(cone_specs) - cones = AbstractCone{T}[] - sizehint!(cones,ncones) + cones = sizehint!(AbstractCone{T}[],length(cone_specs)) type_counts = Dict{Type,Int}() @@ -52,12 +53,14 @@ struct CompositeCone{T} <: AbstractCone{T} numel = sum(cone -> Clarabel.numel(cone), cones; init = 0) degree = sum(cone -> Clarabel.degree(cone), cones; init = 0) - #headidx gives the index of the first element - #of each constituent cone - headidx = Vector{Int}(undef,length(cones)) - _make_headidx!(headidx,cones) + #ranges for the subvectors associated with each cone, + #and the range for the corresponding entries + #in the Hs sparse block + + rng_cones = collect(rng_cones_iterator(cones)); + rng_blocks = collect(rng_blocks_iterator(cones)); - return new(cones,type_counts,numel,degree,headidx,_is_symmetric) + obj = new(cones,type_counts,numel,degree,rng_cones,rng_blocks,_is_symmetric) end end @@ -85,14 +88,56 @@ function get_type_count(cones::CompositeCone{T}, type::Type) where {T} end end -function _make_headidx!(headidx,cones) - if(length(cones) > 0) - #index of first element in each cone - headidx[1] = 1 - for i = 2:length(cones) - headidx[i] = headidx[i-1] + numel(cones[i-1]) - end - end - return nothing + +# ------------------------------------- +# iterators to generate indices into vectors +# in a cone or cone-related blocks in the Hessian +struct RangeConesIterator{T} + cones::Vector{AbstractCone{T}} +end +struct RangeBlocksIterator{T} + cones::Vector{AbstractCone{T}} end + +function rng_cones_iterator(cones::Vector{AbstractCone{T}}) where{T} + RangeConesIterator(cones) +end + +function rng_blocks_iterator(cones::Vector{AbstractCone{T}}) where{T} + RangeBlocksIterator(cones) +end + +Base.length(iter::RangeConesIterator) = length(iter.cones) +Base.length(iter::RangeBlocksIterator) = length(iter.cones) + +function Base.iterate(iter::RangeConesIterator, state=(1, 1)) + (coneidx, start) = state + if coneidx > length(iter.cones) + return nothing + else + nvars = numel(iter.cones[coneidx]) + stop = start + nvars - 1 + state = (coneidx + 1, stop + 1) + return (start:stop, state) + end +end + +function Base.iterate(iter::RangeBlocksIterator, state=(1, 1)) + (coneidx, start) = state + if coneidx > length(iter.cones) + return nothing + else + cone = iter.cones[coneidx] + nvars = numel(cone) + if Hs_is_diagonal(cone) + stop = start + nvars - 1 + else + stop = start + triangular_number(nvars) - 1 + end + state = (coneidx + 1, stop + 1) + return (start:stop, state) + end +end + + diff --git a/src/cones/cone_types.jl b/src/cones/cone_types.jl index 9030217d..48ce140d 100644 --- a/src/cones/cone_types.jl +++ b/src/cones/cone_types.jl @@ -94,7 +94,7 @@ mutable struct SecondOrderCone{T} <: AbstractCone{T} η::T #sparse representation of W^2 - sparse_data::Union{Nothing,SecondOrderConeSparseData{T}} + sparse_data::Option{SecondOrderConeSparseData{T}} function SecondOrderCone{T}(dim::Integer) where {T} @@ -124,15 +124,13 @@ SecondOrderCone(args...) = SecondOrderCone{DefaultFloat}(args...) mutable struct PSDConeData{T} - chol1::Union{Nothing,Cholesky{T,Matrix{T}}} - chol2::Union{Nothing,Cholesky{T,Matrix{T}}} - SVD::Union{Nothing,SVD{T,T,Matrix{T}}} + chol1::Option{Cholesky{T,Matrix{T}}} + chol2::Option{Cholesky{T,Matrix{T}}} + SVD::Option{SVD{T,T,Matrix{T}}} λ::Vector{T} Λisqrt::Diagonal{T,Vector{T}} R::Matrix{T} Rinv::Matrix{T} - kronRR::Matrix{T} - B::Matrix{T} Hs::Matrix{T} #workspace for various internal uses @@ -151,9 +149,7 @@ mutable struct PSDConeData{T} Λisqrt = Diagonal(zeros(T,n)) R = zeros(T,n,n) Rinv = zeros(T,n,n) - kronRR = zeros(T,n^2,n^2) - B = zeros(T,triangular_number(n),n^2) - Hs = zeros(T,size(B,1),size(B,1)) + Hs = zeros(T,triangular_number(n),triangular_number(n)) workmat1 = zeros(T,n,n) workmat2 = zeros(T,n,n) @@ -161,7 +157,7 @@ mutable struct PSDConeData{T} workvec = zeros(T,triangular_number(n)) return new(chol1,chol2,SVD,λ,Λisqrt,R,Rinv, - kronRR,B,Hs,workmat1,workmat2,workmat3,workvec) + Hs,workmat1,workmat2,workmat3,workvec) end end diff --git a/src/cones/coneops_compositecone.jl b/src/cones/coneops_compositecone.jl index 2dbba9a0..13efb1a6 100644 --- a/src/cones/coneops_compositecone.jl +++ b/src/cones/coneops_compositecone.jl @@ -27,8 +27,8 @@ end function rectify_equilibration!( cones::CompositeCone{T}, - δ::ConicVector{T}, - e::ConicVector{T} + δ::Vector{T}, + e::Vector{T} ) where{T} any_changed = false @@ -37,7 +37,9 @@ function rectify_equilibration!( #from this function. default is to do nothing at all δ .= 1 - for (cone,δi,ei) in zip(cones,δ.views,e.views) + for (cone,rng) in zip(cones,cones.rng_cones) + δi = view(δ,rng) + ei = view(e,rng) @conedispatch any_changed |= rectify_equilibration!(cone,δi,ei) end @@ -46,13 +48,13 @@ end function margins( cones::CompositeCone{T}, - z::ConicVector{T}, + z::Vector{T}, pd::PrimalOrDualCone, ) where {T} α = typemax(T) β = zero(T) - for (cone,zi) in zip(cones,z.views) - @conedispatch (αi,βi) = margins(cone,zi,pd) + for (cone,rng) in zip(cones,cones.rng_cones) + @conedispatch (αi,βi) = margins(cone,view(z,rng),pd) α = min(α,αi) β += βi end @@ -62,13 +64,13 @@ end function scaled_unit_shift!( cones::CompositeCone{T}, - z::ConicVector{T}, + z::Vector{T}, α::T, pd::PrimalOrDualCone ) where {T} - for (cone,zi) in zip(cones,z.views) - @conedispatch scaled_unit_shift!(cone,zi,α,pd) + for (cone,rng) in zip(cones,cones.rng_cones) + @conedispatch scaled_unit_shift!(cone,view(z,rng),α,pd) end return nothing @@ -77,12 +79,12 @@ end # unit initialization for asymmetric solves function unit_initialization!( cones::CompositeCone{T}, - z::ConicVector{T}, - s::ConicVector{T} + z::Vector{T}, + s::Vector{T} ) where {T} - for (cone,zi,si) in zip(cones,z.views,s.views) - @conedispatch unit_initialization!(cone,zi,si) + for (cone,rng) in zip(cones,cones.rng_cones) + @conedispatch unit_initialization!(cone,view(z,rng),view(s,rng)) end return nothing end @@ -100,18 +102,16 @@ end function update_scaling!( cones::CompositeCone{T}, - s::ConicVector{T}, - z::ConicVector{T}, + s::Vector{T}, + z::Vector{T}, μ::T, scaling_strategy::ScalingStrategy ) where {T} - # update cone scalings by passing subview to each of - # the appropriate cone types. - for (cone,si,zi) in zip(cones,s.views,z.views) + for (cone,rng) in zip(cones,cones.rng_cones) + si = view(s,rng) + zi = view(z,rng) @conedispatch is_scaling_success = update_scaling!(cone,si,zi,μ,scaling_strategy) - # YC: currently, only check whether SOC variables are in the interior; - # we could extend the interior checkfor other cones if !is_scaling_success return is_scaling_success = false end @@ -122,11 +122,11 @@ end # The Hs block for each cone. function get_Hs!( cones::CompositeCone{T}, - Hsblocks::Vector{Vector{T}} + Hsblock::Vector{T} ) where {T} - for (cone, block) in zip(cones,Hsblocks) - @conedispatch get_Hs!(cone,block) + for (cone, rng) in zip(cones,cones.rng_blocks) + @conedispatch get_Hs!(cone,view(Hsblock,rng)) end return nothing end @@ -137,13 +137,13 @@ end function mul_Hs!( cones::CompositeCone{T}, - y::ConicVector{T}, - x::ConicVector{T}, - work::ConicVector{T} + y::Vector{T}, + x::Vector{T}, + work::Vector{T} ) where {T} - for (cone,yi,xi,worki) in zip(cones,y.views,x.views,work.views) - @conedispatch mul_Hs!(cone,yi,xi,worki) + for (cone,rng) in zip(cones,cones.rng_cones) + @conedispatch mul_Hs!(cone,view(y,rng),view(x,rng),view(work,rng)) end return nothing @@ -152,11 +152,13 @@ end # x = λ ∘ λ for symmetric cone and x = s for asymmetric cones function affine_ds!( cones::CompositeCone{T}, - ds::ConicVector{T}, - s::ConicVector{T} + ds::Vector{T}, + s::Vector{T} ) where {T} - for (cone,dsi,si) in zip(cones,ds.views,s.views) + for (cone,rng) in zip(cones,cones.rng_cones) + dsi = view(ds,rng) + si = view(s,rng) @conedispatch affine_ds!(cone,dsi,si) end return nothing @@ -164,15 +166,16 @@ end function combined_ds_shift!( cones::CompositeCone{T}, - shift::ConicVector{T}, - step_z::ConicVector{T}, - step_s::ConicVector{T}, + shift::Vector{T}, + step_z::Vector{T}, + step_s::Vector{T}, σμ::T ) where {T} - for (cone,shifti,step_zi,step_si) in zip(cones,shift.views,step_z.views,step_s.views) - - # compute the centering and the higher order correction parts in ds and save it in dz + for (cone,rng) in zip(cones,cones.rng_cones) + shifti = view(shift,rng) + step_zi = view(step_z,rng) + step_si = view(step_s,rng) @conedispatch combined_ds_shift!(cone,shifti,step_zi,step_si,σμ) end @@ -181,13 +184,17 @@ end function Δs_from_Δz_offset!( cones::CompositeCone{T}, - out::ConicVector{T}, - ds::ConicVector{T}, - work::ConicVector{T}, - z::ConicVector{T} + out::Vector{T}, + ds::Vector{T}, + work::Vector{T}, + z::Vector{T} ) where {T} - for (cone,outi,dsi,worki,zi) in zip(cones,out.views,ds.views,work.views,z.views) + for (cone,rng) in zip(cones,cones.rng_cones) + outi = view(out,rng) + dsi = view(ds,rng) + worki = view(work,rng) + zi = view(z,rng) @conedispatch Δs_from_Δz_offset!(cone,outi,dsi,worki,zi) end @@ -197,27 +204,31 @@ end # maximum allowed step length over all cones function step_length( cones::CompositeCone{T}, - dz::ConicVector{T}, - ds::ConicVector{T}, - z::ConicVector{T}, - s::ConicVector{T}, + dz::Vector{T}, + ds::Vector{T}, + z::Vector{T}, + s::Vector{T}, settings::Settings{T}, αmax::T, ) where {T} - α = αmax - dz = dz.views - ds = ds.views - z = z.views - s = s.views + α = αmax + + function innerfcn(α,symcond) + for (cone,rng) in zip(cones,cones.rng_cones) + if @conedispatch is_symmetric(cone) == symcond + continue + end + (dzi,dsi) = (view(dz,rng),view(ds,rng)) + (zi,si) = (view(z,rng),view(s,rng)) + @conedispatch (nextαz,nextαs) = step_length(cone,dzi,dsi,zi,si,settings,α) + α = min(α,nextαz,nextαs) + end + return α + end # Force symmetric cones first. - for (cone,dzi,dsi,zi,si) in zip(cones,dz,ds,z,s) - - if !is_symmetric(cone) continue end - @conedispatch (nextαz,nextαs) = step_length(cone,dzi,dsi,zi,si,settings,α) - α = min(α,nextαz,nextαs) - end + α = innerfcn(α,true) #if we have any nonsymmetric cones, then back off from full steps slightly #so that centrality checks and logarithms don't fail right at the boundaries @@ -225,13 +236,8 @@ function step_length( α = min(α,settings.max_step_fraction) end - # Force asymmetric cones last. - for (cone,dzi,dsi,zi,si) in zip(cones,dz,ds,z,s) - - if @conedispatch(is_symmetric(cone)) continue end - @conedispatch (nextαz,nextαs) = step_length(cone,dzi,dsi,zi,si,settings,α) - α = min(α,nextαz,nextαs) - end + # now the nonsymmetric cones + α = innerfcn(α,false) return (α,α) end @@ -239,21 +245,19 @@ end # compute the total barrier function at the point (z + α⋅dz, s + α⋅ds) function compute_barrier( cones::CompositeCone{T}, - z::ConicVector{T}, - s::ConicVector{T}, - dz::ConicVector{T}, - ds::ConicVector{T}, + z::Vector{T}, + s::Vector{T}, + dz::Vector{T}, + ds::Vector{T}, α::T ) where {T} - - dz = dz.views - ds = ds.views - z = z.views - s = s.views - barrier = zero(T) - for (cone,zi,si,dzi,dsi) in zip(cones,z,s,dz,ds) + for (cone,rng) in zip(cones,cones.rng_cones) + zi = view(z,rng) + si = view(s,rng) + dzi = view(dz,rng) + dsi = view(ds,rng) @conedispatch barrier += compute_barrier(cone,zi,si,dzi,dsi,α) end diff --git a/src/cones/coneops_genpowcone.jl b/src/cones/coneops_genpowcone.jl index fcd259ef..4790edd7 100644 --- a/src/cones/coneops_genpowcone.jl +++ b/src/cones/coneops_genpowcone.jl @@ -96,8 +96,9 @@ function get_Hs!( #NB: we are returning here the diagonal D = [d1; d2] block from the #sparse representation of W^TW, but not the #extra 3 entries at the bottom right of the block. - #The ConicVector for s and z (and its views) don't + #The vectors for s and z (and its views) don't #know anything about the 3 extra sparsifying entries + dim1 = Clarabel.dim1(K) data = K.data diff --git a/src/cones/coneops_psdtrianglecone.jl b/src/cones/coneops_psdtrianglecone.jl index 52c52536..a235817a 100644 --- a/src/cones/coneops_psdtrianglecone.jl +++ b/src/cones/coneops_psdtrianglecone.jl @@ -18,7 +18,7 @@ function margins( α = floatmax(T) else Z = K.data.workmat1 - _svec_to_mat!(Z,z) + svec_to_mat!(Z,z) e = eigvals!(Hermitian(Z)) #NB: GenericLinearAlgebra doesn't support eigvals!(::Symmetric(...)) α = minimum(e) #minimum eigenvalue. end @@ -92,7 +92,7 @@ function update_scaling!( f = K.data (S,Z) = (f.workmat1,f.workmat2) - map((M,v)->_svec_to_mat!(M,v),(S,Z),(s,z)) + map((M,v)->svec_to_mat!(M,v),(S,Z),(s,z)) #compute Cholesky factors (PG: this is allocating) f.chol1 = cholesky!(S, check = false) @@ -116,54 +116,32 @@ function update_scaling!( #f.R = L1*(f.SVD.V)*f.Λisqrt mul!(f.R,L1,f.SVD.V); - mul!(f.R,f.R,f.Λisqrt) #mul! can take Rinv twice because Λ is diagonal + mul!(f.R,f.R,f.Λisqrt) #mul! can take R twice because Λ is diagonal #f.Rinv .= f.Λisqrt*(f.SVD.U)'*L2' mul!(f.Rinv,f.SVD.U',L2') mul!(f.Rinv,f.Λisqrt,f.Rinv) #mul! can take Rinv twice because Λ is diagonal - - # PJG: The following steps force us to form Hs in memory - # in the scaling update, even if we aren't using a - # direct method and therefore never actually require - # the matrix Hs to be formed. The steps below should - # be simplified if possible and then only implemented - # within get_Hs, placing the matrix directly into the - # diagonal Hs block provided by the direct factorizer. - - # we should compute here the upper triangular part - # of the matrix Q* ((RR^T) ⨂ (RR^T)) * P. The operator - # P is a matrix that transforms a packed triangle to - # a vectorized full matrix. Q sends it back. - # - # See notes by Kathrin Schäcke, 2013: "On the Kronecker Product" - # for some useful identities, particularly section 3 on symmetric - # Kronecker product - - @inbounds kron!(f.kronRR,f.R,f.R) - - #B .= Q'*kRR, where Q' is the svec operator - #this could be substantially faster - for i = 1:size(f.B,2) - @views M = reshape(f.kronRR[:,i],size(f.R,1),size(f.R,1)) - b = view(f.B,:,i) - _mat_to_svec!(b,M) - end - - #compute Hs = triu(B*B') - # PJG: I pack this into triu form by calling - # pack_triu with get_Hs. Would be ideal - # if this could be done directly, but it's - # not clear how to do so via blas. + #compute R*R' (upper triangular part only) + RRt = f.workmat1; + RRt .= zero(T) if T <: LinearAlgebra.BlasFloat - LinearAlgebra.BLAS.syrk!('U', 'N', one(T), f.B, zero(T), f.Hs) + LinearAlgebra.BLAS.syrk!('U', 'N', one(T), f.R, zero(T), RRt) else - f.Hs .= f.B*f.B' + RRt .= f.R*f.R' end + # PJG: it is possibly faster to compute the whole of RRt, and not + # just the upper triangle using syrk!, because then skron! can be + # be called with a Matrix type instead of Symmetric. The internal + # indexing within skron! is then more straightforward and probably + # faster. + skron!(f.Hs,Symmetric(RRt)) + return is_scaling_success = true end + function Hs_is_diagonal( K::PSDTriangleCone{T} ) where{T} @@ -198,11 +176,10 @@ function mul_Hs!( # way and then *never* hold Hs in the cone work data. # Perhaps the calculation of scale factors that includes # the kron(R,R) onwards could be done in a more compact - # way since that internal Hs work variable is only really + # way since the internal Hs work variable is only really # needed to populate the KKT Hs block. For a direct # method that block is never needed, so better to only - # form it in memory if get_Hs is actually called and - # provides a place for it. + # form it in memory if get_Hs is actually called mul_W!(K,:N,work,x,one(T),zero(T)) #work = Wx mul_W!(K,:T,y,work,one(T),zero(T)) #y = Wᵀwork = W^TWx @@ -265,11 +242,11 @@ function step_length( #d = Δz̃ = WΔz mul_W!(K, :N, d, dz, one(T), zero(T)) - αz = _step_length_psd_component(workΔ,d,Λisqrt,αmax) + αz = step_length_psd_component(workΔ,d,Λisqrt,αmax) #d = Δs̃ = W^{-T}Δs mul_Winv!(K, :T, d, ds, one(T), zero(T)) - αs = _step_length_psd_component(workΔ,d,Λisqrt,αmax) + αs = step_length_psd_component(workΔ,d,Λisqrt,αmax) return (αz,αs) end @@ -297,7 +274,7 @@ function _logdet_barrier(K::PSDTriangleCone{T},x::AbstractVector{T},dx::Abstract q = f.workvec @. q = x + alpha*dx - _svec_to_mat!(Q,q) + svec_to_mat!(Q,q) # PG: this is allocating f.chol1 = cholesky!(Q, check = false) @@ -326,7 +303,7 @@ function mul_W!( β::T ) where {T} - _mul_Wx_inner( + mul_Wx_inner( is_transpose, y,x, α, @@ -347,7 +324,7 @@ function mul_Winv!( β::T ) where {T} - _mul_Wx_inner( + mul_Wx_inner( is_transpose, y,x, α, @@ -367,7 +344,7 @@ function λ_inv_circ_op!( ) where {T} (X,Z) = (K.data.workmat1,K.data.workmat2) - map((M,v)->_svec_to_mat!(M,v),(X,Z),(x,z)) + map((M,v)->svec_to_mat!(M,v),(X,Z),(x,z)) λ = K.data.λ for i = 1:K.n @@ -375,7 +352,7 @@ function λ_inv_circ_op!( X[i,j] = 2*Z[i,j]/(λ[i] + λ[j]) end end - _mat_to_svec!(x,X) + mat_to_svec!(x,X) return nothing end @@ -393,7 +370,7 @@ function circ_op!( ) where {T} (Y,Z) = (K.data.workmat1,K.data.workmat2) - map((M,v)->_svec_to_mat!(M,v),(Y,Z),(y,z)) + map((M,v)->svec_to_mat!(M,v),(Y,Z),(y,z)) X = K.data.workmat3; @@ -404,7 +381,7 @@ function circ_op!( else X .= (Y*Z + Z*Y)/2 end - _mat_to_svec!(x,Symmetric(X)) + mat_to_svec!(x,Symmetric(X)) return nothing end @@ -433,7 +410,7 @@ end # internal operations for SDP cones # ---------------------------------------- -function _mul_Wx_inner( +function mul_Wx_inner( is_transpose::Symbol, y::AbstractVector{T}, x::AbstractVector{T}, @@ -446,7 +423,7 @@ function _mul_Wx_inner( ) where {T} (X,Y,tmp) = (workmat1,workmat2,workmat3) - map((M,v)->_svec_to_mat!(M,v),(X,Y),(x,y)) + map((M,v)->svec_to_mat!(M,v),(X,Y),(x,y)) if is_transpose === :T #Y .= α*(R*X*R') #W^T*x or.... @@ -460,12 +437,12 @@ function _mul_Wx_inner( mul!(Y,tmp,Rx,α,β) end - _mat_to_svec!(y,Y) + mat_to_svec!(y,Y) return nothing end -function _step_length_psd_component( +function step_length_psd_component( workΔ::Matrix{T}, d::Vector{T}, Λisqrt::Diagonal{T}, @@ -476,7 +453,7 @@ function _step_length_psd_component( γ = floatmax(T) else # NB: this could be made faster since we only need to populate the upper triangle - _svec_to_mat!(workΔ,d) + svec_to_mat!(workΔ,d) lrscale!(Λisqrt.diag,workΔ,Λisqrt.diag) # GenericLinearAlgebra doesn't support eigvals!(::Symmetric(::Matrix)), # and doesn't support choosing a subset of values @@ -496,12 +473,12 @@ function _step_length_psd_component( end #make a matrix view from a vectorized input -function _svec_to_mat!( M::AbstractMatrix{T}, x::AbstractVector{T}) where {T} +function svec_to_mat!( M::AbstractMatrix{T}, x::AbstractVector{T}) where {T} ISQRT2 = inv(sqrt(T(2))) idx = 1 - for col = 1:size(M,2), row = 1:col + for col in axes(M,2), row in 1:col if row == col M[row,col] = x[idx] else @@ -513,12 +490,12 @@ function _svec_to_mat!( M::AbstractMatrix{T}, x::AbstractVector{T}) where {T} end -function _mat_to_svec!(x::AbstractVector{T},M::AbstractMatrix{T}) where {T} +function mat_to_svec!(x::AbstractVector{T},M::AbstractMatrix{T}) where {T} ISQRT2 = 1/sqrt(T(2)) idx = 1 - for col = 1:size(M,2), row = 1:col + for col in axes(M,2), row in 1:col @inbounds x[idx] = row == col ? M[row,col] : (M[row,col]+M[col,row])*ISQRT2 idx += 1 end @@ -526,3 +503,45 @@ function _mat_to_svec!(x::AbstractVector{T},M::AbstractMatrix{T}) where {T} return nothing end + +# produce the upper triangular part of the Symmetric Kronecker product of +# a symmtric matrix A with itself, i.e. triu(A ⊗_s A) +function skron!( + out::Matrix{T}, + A::Symmetric{T, Matrix{T}}, +) where {T} + + sqrt2 = sqrt(2) + n = size(A, 1) + + col = 1 + for l in 1:n + for k in 1:l + row = 1 + kl_eq = k == l + + @inbounds for j in 1:n + Ajl = A[j, l] + Ajk = A[j, k] + + @inbounds for i in 1:j + (row > col) && break + ij_eq = i == j + + if (ij_eq, kl_eq) == (false, false) + out[row, col] = A[i, k] * Ajl + A[i, l] * Ajk + elseif (ij_eq, kl_eq) == (true, false) + out[row, col] = sqrt2 * Ajl * Ajk + elseif (ij_eq, kl_eq) == (false, true) + out[row, col] = sqrt2 * A[i, l] * Ajk + else (ij_eq,kl_eq) == (true, true) + out[row, col] = Ajl * Ajl + end + + row += 1 + end # i + end # j + col += 1 + end # k + end # l +end \ No newline at end of file diff --git a/src/conicvector.jl b/src/conicvector.jl deleted file mode 100644 index f0f0e25d..00000000 --- a/src/conicvector.jl +++ /dev/null @@ -1,99 +0,0 @@ -# ------------------------------------- -# vectors defined w.r.t. to conic constraints -# get this type with views into the subcomponents -# --------------------------------------- - -const VectorView{T} = SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true} - -struct ConicVector{T<:AbstractFloat} <: AbstractVector{T} - - #contiguous array of source data - vec::Vector{T} - - #array of data views of type Vector{T} - views::Vector{VectorView{T}} - - function ConicVector{T}(S::CompositeCone{T}) where {T} - - #undef initialization would possibly result - #in Infs or NaNs, causing failure in gemv! - #style vector updates - vec = zeros(T,S.numel) - views = Vector{VectorView{T}}(undef, length(S)) - - # loop over the sets and create views - last = 0 - for i = eachindex(S) - first = last + 1 - last = last + numel(S[i]) - views[i] = view(vec, first:last) - end - - return new(vec, views) - - end - -end - -ConicVector(args...) = ConicVector{DefaultFloat}(args...) - - - - -@inline function Base.getindex(s::ConicVector{T},i) where{T} - @boundscheck checkbounds(s.vec,i) - @inbounds s.vec[i] -end -@inline function Base.setindex!(s::ConicVector{T},v,i) where{T} - @boundscheck checkbounds(s.vec,i) - @inbounds s.vec[i] = T(v) -end - -Base.adjoint(s::ConicVector{T}) where{T} = adjoint(s.vec) -Base.iterate(s::ConicVector{T}) where{T} = iterate(s.vec) -Base.eltype(s::ConicVector{T}) where{T} = eltype(s.vec) -Base.elsize(s::ConicVector{T}) where{T} = Base.elsize(s.vec) -Base.size(s::ConicVector{T}) where{T} = size(s.vec) -Base.length(s::ConicVector{T}) where{T} = length(s.vec) -Base.IndexStyle(s::ConicVector{T}) where{T} = IndexStyle(s.vec) - -#For maximum speed, it seems we need to explicitly define -#a bunch of functions that use the vec field directly, which -#will force calls to the BLAS specialized methods when possible -#Alternatively, we could subtype DenseArray{T}, but that -#seems less general and still fails to capture high -#performance sparse matrix vector multiply - - -#need this if we want to make calls directly to the BLAS functions -Base.unsafe_convert(::Type{Ptr{T}}, s::ConicVector{T}) where {T} = - Base.unsafe_convert(Ptr{T}, s.vec) - -#dot -LinearAlgebra.dot(x::AbstractVector{T},y::ConicVector{T}) where {T} = dot(x,y.vec) -LinearAlgebra.dot(x::ConicVector{T},y::AbstractVector{T}) where {T} = dot(x.vec,y) -LinearAlgebra.dot(x::ConicVector{T},y::ConicVector{T}) where {T} = dot(x.vec,y.vec) - -#mul! -LinearAlgebra.mul!( - C::ConicVector, - A::AbstractVecOrMat, - B::AbstractVector, α::Number, β::Number) = mul!(C.vec, A, B, α, β) - -LinearAlgebra.mul!( - C::AbstractVector, - A::AbstractVecOrMat, - B::ConicVector, α::Number, β::Number) = mul!(C, A, B.vec, α, β) - -LinearAlgebra.mul!( - C::ConicVector, - A::Adjoint{<:Any, <:AbstractVecOrMat}, - B::AbstractVector, α::Number, β::Number) = mul!(C.vec, A, B, α, β) - -LinearAlgebra.mul!( - C::AbstractVector, - A::Adjoint{<:Any, <:AbstractVecOrMat}, - B::ConicVector, α::Number, β::Number) = mul!(C, A, B.vec, α, β) - -LinearAlgebra.norm(x::ConicVector{T}) where {T} = norm(x.vec) -LinearAlgebra.norm(x::ConicVector{T},p::Real) where {T} = norm(x.vec,p) \ No newline at end of file diff --git a/src/data_updating.jl b/src/data_updating.jl index e18e5ee9..f540c6a5 100644 --- a/src/data_updating.jl +++ b/src/data_updating.jl @@ -18,14 +18,15 @@ const VectorProblemDataUpdate{T} = Union{ update_data!(solver,P,q,A,b) Overwrites internal problem data structures in a solver object with new data, avoiding new memory -allocations. See [`update_P!`](@ref), [`update_q!`](@ref), [`update_A!`](@ref), [`update_b!`](@ref) for allowable input types. +allocations. See [`update_P!`](@ref), [`update_q!`](@ref), [`update_A!`](@ref), [`update_b!`](@ref) +for allowable input types. """ function update_data!( s::Solver{T}, P::VectorProblemDataUpdate{T} , - q::Union{Vector{T},Nothing}, + q::Option{Vector{T}}, A::MatrixProblemDataUpdate{T}, b::VectorProblemDataUpdate{T} ) where{T} @@ -58,7 +59,7 @@ function update_P!( ) where{T} isnothing(data) && return - _check_presolve_disabled(s) + _check_update_allowed(s) d = s.data.equilibration.d _update_matrix(data,s.data.P,d,d) # overwrite KKT data @@ -87,9 +88,9 @@ function update_A!( ) where{T} isnothing(data) && return - _check_presolve_disabled(s) + _check_update_allowed(s) d = s.data.equilibration.d - e = s.data.equilibration.e.vec #e is a conic vector type + e = s.data.equilibration.e _update_matrix(data,s.data.A,e,d) # overwrite KKT data kkt_update_A!(s.kktsystem,s.data.A) @@ -110,7 +111,7 @@ function update_q!( ) where{T} isnothing(data) && return - _check_presolve_disabled(s) + _check_update_allowed(s) d = s.data.equilibration.d dinv = s.data.equilibration.dinv _update_vector(data,s.data.q,d) @@ -134,9 +135,9 @@ function update_b!( ) where{T} isnothing(data) && return - _check_presolve_disabled(s) - e = s.data.equilibration.e.vec #e is a conic vector type - einv = s.data.equilibration.einv.vec #einv is a conic vector type + _check_update_allowed(s) + e = s.data.equilibration.e + einv = s.data.equilibration.einv _update_vector(data,s.data.b,e) # flush unscaled norm. Will be recalculated during solve @@ -145,12 +146,23 @@ function update_b!( return nothing end -function _check_presolve_disabled(s) - # Fail if presolve is enabled even if the sparsity is the same. - # Not strictly necessary but may avoid confusion about expectations. - if s.settings.presolve_enable - error("Disable presolve to allow data updates.") +function _check_update_allowed(s) + + # Fail if presolve / chordal decomp is enabled. + # Not strictly necessary since the presolve and chordal decomp + # might not do anything, but may avoid confusion about expectations. + + # checks both settings and existence of presolve objects, since otherwise + # it would be possible to presolve and then disable the settings. + + if s.settings.presolve_enable || + s.settings.chordal_decomposition_enable || + !isnothing(s.data.presolver) || + !isnothing(s.data.chordal_info) + + error("Disable presolve and chordal decomposition to allow data updates.") end + end function _update_matrix( diff --git a/src/info.jl b/src/info.jl index 9b3e9e11..0774c41e 100644 --- a/src/info.jl +++ b/src/info.jl @@ -193,11 +193,10 @@ function info_get_solve_time!( end -function info_finalize!( +function info_post_process!( info::DefaultInfo{T}, residuals::DefaultResiduals{T}, settings::Settings{T}, - timers::TimerOutput ) where {T} # if there was an error or we ran out of time @@ -210,6 +209,13 @@ function info_finalize!( _check_convergence_almost(info,residuals,settings) end +end + +function info_finalize!( + info::DefaultInfo{T}, + timers::TimerOutput +) where {T} + # final check of timers info_get_solve_time!(info,timers) return nothing diff --git a/src/info_print.jl b/src/info_print.jl index 5a517a9e..0c2cba77 100644 --- a/src/info_print.jl +++ b/src/info_print.jl @@ -21,10 +21,14 @@ function info_print_configuration( if(settings.verbose == false) return end - if(is_reduced(data.presolver)) + if(!isnothing(data.presolver)) @printf(io, "\npresolve: removed %i constraints\n", count_reduced(data.presolver)) end + if(!isnothing(data.chordal_info)) + print_chordal_decomposition(io, data.chordal_info, settings) + end + @printf(io, "\nproblem:\n") @printf(io, " variables = %i\n", data.n) @printf(io, " constraints = %i\n", data.m) @@ -43,6 +47,7 @@ function info_print_configuration( return nothing end + info_print_configuration(info,settings,data,cones) = info_print_configuration(stdout,info,settings,data,cones) @@ -176,6 +181,26 @@ function print_settings(io::IO, settings::Settings{T}) where {T} return nothing end +function print_chordal_decomposition(io::IO, chordal_info::ChordalInfo{T}, settings::Settings{T}) where {T} + + @printf(io, "\nchordal decomposition:\n") + @printf(io, " compact format = %s, dual completion = %s\n", + bool_on_off(settings.chordal_decomposition_compact), + bool_on_off(settings.chordal_decomposition_complete_dual)) + + @printf(io, " merge method = %s\n", settings.chordal_decomposition_merge_method) + + @printf(io, " PSD cones initial = %d\n", + init_psd_cone_count(chordal_info)) + @printf(io, " PSD cones decomposable = %d\n", + decomposable_cone_count(chordal_info)) + @printf(io, " PSD cones after decomposition = %d\n", + premerge_psd_cone_count(chordal_info)) + @printf(io, " PSD cones after merges = %d\n", + final_psd_cone_count(chordal_info)) + +end + get_precision_string(T::Type{<:Real}) = string(T) get_precision_string(T::Type{<:BigFloat}) = string(T," (", precision(T), " bit)") diff --git a/src/kktsolvers/direct-ldl/directldl_cholmod.jl b/src/kktsolvers/direct-ldl/directldl_cholmod.jl index 2242f6c3..0a79aad8 100644 --- a/src/kktsolvers/direct-ldl/directldl_cholmod.jl +++ b/src/kktsolvers/direct-ldl/directldl_cholmod.jl @@ -1,7 +1,7 @@ using SuiteSparse mutable struct CholmodDirectLDLSolver{T} <: AbstractDirectLDLSolver{T} - F::Union{SuiteSparse.CHOLMOD.Factor,Nothing} + F::Option{SuiteSparse.CHOLMOD.Factor} function CholmodDirectLDLSolver{T}(KKT::SparseMatrixCSC{T},Dsigns,settings) where {T} diff --git a/src/kktsolvers/direct-ldl/directldl_datamaps.jl b/src/kktsolvers/direct-ldl/directldl_datamaps.jl index 71a96cf7..911c7979 100644 --- a/src/kktsolvers/direct-ldl/directldl_datamaps.jl +++ b/src/kktsolvers/direct-ldl/directldl_datamaps.jl @@ -10,8 +10,8 @@ struct SOCExpansionMap <: SparseExpansionMap v::Vector{Int} #off diag dense columns v D::MVector{2, Int} #diag D function SOCExpansionMap(cone::SecondOrderCone) - u = Vector{Int}(undef,numel(cone)) - v = Vector{Int}(undef,numel(cone)) + u = zeros(Int,numel(cone)) + v = zeros(Int,numel(cone)) D = MVector(0,0) new(u,v,D) end @@ -86,9 +86,9 @@ struct GenPowExpansionMap <: SparseExpansionMap D::MVector{3, Int} #diag D function GenPowExpansionMap(cone::GenPowerCone) - p = Vector{Int}(undef,numel(cone)) - q = Vector{Int}(undef,dim1(cone)) - r = Vector{Int}(undef,dim2(cone)) + p = zeros(Int,numel(cone)) + q = zeros(Int,dim1(cone)) + r = zeros(Int,dim2(cone)) D = MVector(0,0,0) new(p,q,r,D) end @@ -171,7 +171,7 @@ struct LDLDataMap P::Vector{Int} A::Vector{Int} - Hsblocks::Vector{Vector{Int}} #indices of the lower RHS blocks (by cone) + Hsblocks::Vector{Int} #indices of the lower RHS blocks (by cone) sparse_maps::Vector{SparseExpansionMap} #sparse cone expansion terms #all of above terms should be disjoint and their union @@ -198,8 +198,7 @@ struct LDLDataMap #now do the sparse cone expansion pieces nsparse = count(cone->(@conedispatch is_sparse_expandable(cone)),cones) - sparse_maps = Vector{SparseExpansionMap}(); - sizehint!(sparse_maps,nsparse) + sparse_maps = sizehint!(SparseExpansionMap[],nsparse) for cone in cones if @conedispatch is_sparse_expandable(cone) diff --git a/src/kktsolvers/direct-ldl/directldl_kkt_assembly.jl b/src/kktsolvers/direct-ldl/directldl_kkt_assembly.jl index d8d85fee..a24b3fbc 100644 --- a/src/kktsolvers/direct-ldl/directldl_kkt_assembly.jl +++ b/src/kktsolvers/direct-ldl/directldl_kkt_assembly.jl @@ -1,21 +1,14 @@ using SparseArrays, StaticArrays -function _allocate_kkt_Hsblocks(type::Type{T}, cones) where{T <: Real} +function _allocate_kkt_Hsblocks( + ::Type{Z}, + cones::CompositeCone{T} +) where{T <: AbstractFloat, Z <: Real} - ncones = length(cones) - Hsblocks = Vector{Vector{T}}(undef,ncones) + rng_blocks = cones.rng_blocks + nnz = length(rng_blocks) == 0 ? 0 : last(rng_blocks[end]) + zeros(Z,nnz) - for (i, cone) in enumerate(cones) - nvars = numel(cone) - if Hs_is_diagonal(cone) - numelblock = nvars - else #dense triangle - numelblock = triangular_number(nvars) #must be Int - end - Hsblocks[i] = Vector{T}(undef,numelblock) - end - - return Hsblocks end @@ -84,11 +77,11 @@ function _kkt_assemble_colcounts( sparse_map_iter = Iterators.Stateful(map.sparse_maps) for (i,cone) = enumerate(cones) - row = cones.headidx[i] + n + row = first(cones.rng_cones[i]) + n #add the the Hs blocks in the lower right - blockdim = numel(cone) - if Hs_is_diagonal(cone) + @conedispatch blockdim = numel(cone) + if @conedispatch Hs_is_diagonal(cone) _csc_colcount_diag(K,row,blockdim) else _csc_colcount_dense_triangle(K,row,blockdim,shape) @@ -137,15 +130,16 @@ function _kkt_assemble_fill( sparse_map_iter = Iterators.Stateful(map.sparse_maps) for (i,cone) = enumerate(cones) - row = cones.headidx[i] + n + row = first(cones.rng_cones[i]) + n #add the the Hs blocks in the lower right - blockdim = numel(cone) + @conedispatch blockdim = numel(cone) + block = view(map.Hsblocks,cones.rng_blocks[i]) - if Hs_is_diagonal(cone) - _csc_fill_diag(K,map.Hsblocks[i],row,blockdim) + if @conedispatch Hs_is_diagonal(cone) + _csc_fill_diag(K,block,row,blockdim) else - _csc_fill_dense_triangle(K,map.Hsblocks[i],row,blockdim,shape) + _csc_fill_dense_triangle(K,block,row,blockdim,shape) end #add sparse expansions columns for sparse cones diff --git a/src/kktsolvers/kktsolver_defaults.jl b/src/kktsolvers/kktsolver_defaults.jl index b9982b4e..98964e6e 100644 --- a/src/kktsolvers/kktsolver_defaults.jl +++ b/src/kktsolvers/kktsolver_defaults.jl @@ -17,8 +17,8 @@ end #solve and assign LHS function kktsolver_solve!( kktsolver::AbstractKKTSolver{T}, - x::Union{Nothing,AbstractVector{T}}, - z::Union{Nothing,AbstractVector{T}} + x::Option{AbstractVector{T}}, + z::Option{AbstractVector{T}} ) where{T} error("function not implemented") end diff --git a/src/kktsolvers/kktsolver_directldl.jl b/src/kktsolvers/kktsolver_directldl.jl index d5236b7b..9cf23459 100644 --- a/src/kktsolvers/kktsolver_directldl.jl +++ b/src/kktsolvers/kktsolver_directldl.jl @@ -24,7 +24,7 @@ mutable struct DirectLDLKKTSolver{T} <: AbstractKKTSolver{T} # a vector for storing the Hs blocks # on the in the KKT matrix block diagonal - Hsblocks::Vector{Vector{T}} + Hsblocks::Vector{T} #unpermuted KKT matrix KKT::SparseMatrixCSC{T,Int} @@ -43,7 +43,13 @@ mutable struct DirectLDLKKTSolver{T} <: AbstractKKTSolver{T} diagonal_regularizer::T - function DirectLDLKKTSolver{T}(P,A,cones,m,n,settings) where {T} + function DirectLDLKKTSolver{T}( + P::SparseMatrixCSC{T}, + A::SparseMatrixCSC{T}, + cones::CompositeCone{T}, + m::Int,n::Int, + settings::Settings{T} + ) where {T} # get a constructor for the LDL solver we should use, # and also the matrix shape it requires @@ -56,13 +62,13 @@ mutable struct DirectLDLKKTSolver{T} <: AbstractKKTSolver{T} p = pdim(map.sparse_maps) #LHS/RHS/work for iterative refinement - x = Vector{T}(undef,n+m+p) - b = Vector{T}(undef,n+m+p) - work_e = Vector{T}(undef,n+m+p) - work_dx = Vector{T}(undef,n+m+p) + x = zeros(T,n+m+p) + b = zeros(T,n+m+p) + work_e = zeros(T,n+m+p) + work_dx = zeros(T,n+m+p) #the expected signs of D in LDL - Dsigns = Vector{Int}(undef,n+m+p) + Dsigns = zeros(Int,n+m+p) _fill_Dsigns!(Dsigns,m,n,map) #updates to the diagonal of KKT will be @@ -216,11 +222,10 @@ function _kktsolver_update_inner!( #Set the elements the W^tW blocks in the KKT matrix. get_Hs!(cones,kktsolver.Hsblocks) - for (index, values) in zip(map.Hsblocks,kktsolver.Hsblocks) - #change signs to get -W^TW - @. values *= -one(T) - _update_values!(ldlsolver,KKT,index,values) - end + (values, index) = (kktsolver.Hsblocks, map.Hsblocks) + #change signs to get -W^TW + @. values *= -one(T) + _update_values!(ldlsolver,KKT,index,values) sparse_map_iter = Iterators.Stateful(map.sparse_maps) @@ -324,8 +329,8 @@ end function kktsolver_getlhs!( kktsolver::DirectLDLKKTSolver{T}, - lhsx::Union{Nothing,AbstractVector{T}}, - lhsz::Union{Nothing,AbstractVector{T}} + lhsx::Option{AbstractVector{T}}, + lhsz::Option{AbstractVector{T}} ) where {T} x = kktsolver.x @@ -340,8 +345,8 @@ end function kktsolver_solve!( kktsolver::DirectLDLKKTSolver{T}, - lhsx::Union{Nothing,AbstractVector{T}}, - lhsz::Union{Nothing,AbstractVector{T}} + lhsx::Option{AbstractVector{T}}, + lhsz::Option{AbstractVector{T}} ) where {T} (x,b) = (kktsolver.x,kktsolver.b) diff --git a/src/kktsystem.jl b/src/kktsystem.jl index 9b58abea..76d5f133 100644 --- a/src/kktsystem.jl +++ b/src/kktsystem.jl @@ -17,8 +17,8 @@ mutable struct DefaultKKTSystem{T} <: AbstractKKTSystem{T} #work vectors for assembling/disassembling vectors workx::Vector{T} - workz::ConicVector{T} - work_conic::ConicVector{T} + workz::Vector{T} + work_conic::Vector{T} function DefaultKKTSystem{T}( data::DefaultProblemData{T}, @@ -33,19 +33,19 @@ mutable struct DefaultKKTSystem{T} <: AbstractKKTSystem{T} kktsolver = DirectLDLKKTSolver{T}(data.P,data.A,cones,m,n,settings) #the LHS constant part of the reduced solve - x1 = Vector{T}(undef,n) - z1 = Vector{T}(undef,m) + x1 = zeros(T,n) + z1 = zeros(T,m) #the LHS for other solves - x2 = Vector{T}(undef,n) - z2 = Vector{T}(undef,m) + x2 = zeros(T,n) + z2 = zeros(T,m) #workspace compatible with (x,z) - workx = Vector{T}(undef,n) - workz = ConicVector{T}(cones) + workx = zeros(T,n) + workz = zeros(T,m) #additional conic workspace vector compatible with s and z - work_conic = ConicVector{T}(cones) + work_conic = zeros(T,m) return new(kktsolver,x1,z1,x2,z2,workx,workz,work_conic) diff --git a/src/presolver.jl b/src/presolver.jl index 85ed1d83..c12473dc 100644 --- a/src/presolver.jl +++ b/src/presolver.jl @@ -1,119 +1,148 @@ -struct PresolverRowReductionIndex +function Presolver{T}( + A::AbstractMatrix{T}, + b::Vector{T}, + cones::Vector{SupportedCone}, + settings::Settings{T} +) where {T} - # vector of length = original RHS. Entries are false - # for those rows that should be eliminated before solve - keep_logical::Vector{Bool} + infbound = Clarabel.get_infinity() - # vector of length = reduced RHS, taking values - # that map reduced b back to their original index - # This is just findall(keep_logical) and is held for - # efficient solution repopulation - keep_index::Vector{Int64} + # make copy of cones to protect from user interference + init_cones = Vector{SupportedCone}(cones) -end + mfull = length(b) -struct Presolver{T} - - # possibly reduced internal copy of user cone specification - cone_specs::Vector{SupportedCone} - - # record of reduced constraints for NN cones with inf bounds - reduce_map::Union{Nothing,PresolverRowReductionIndex} - - # size of original and reduced RHS, respectively - mfull::Int64 - mreduced::Int64 - - # inf bound that was taken from the module level - # and should be applied throughout. Held here so - # that any subsequent change to the module's state - # won't mess up our solver mid-solve - infbound::Float64 - - function Presolver{T}( - A::AbstractMatrix{T}, - b::Vector{T}, - cone_specs::Vector{<:SupportedCone}, - settings::Settings{T} - ) where {T} - - infbound = Clarabel.get_infinity() - - # make copy of cone_specs to protect from user interference after - # setup and explicitly cast as the abstract type. This also prevents - # errors arising from input vectors that are all the same cone, and - # therefore more concretely typed than we want - cone_specs = Vector{SupportedCone}(cone_specs) - mfull = length(b) - - (reduce_map, mreduced) = let - if settings.presolve_enable - reduce_cones!(cone_specs,b,T(infbound)) - else - (nothing,mfull) - end - end - - return new(cone_specs,reduce_map,mfull, mreduced, infbound) + (reduce_map, mreduced) = make_reduction_map(cones,b,T(infbound)) - end + return Presolver{T}(init_cones, reduce_map, mfull, mreduced, infbound) end -Presolver(args...) = Presolver{DefaultFloat}(args...) is_reduced(ps::Presolver{T}) where {T} = !isnothing(ps.reduce_map) count_reduced(ps::Presolver{T}) where {T} = ps.mfull - ps.mreduced -function reduce_cones!( - cone_specs::Vector{<:SupportedCone}, +function presolve( + presolver::Presolver{T}, + A::SparseMatrixCSC{T}, + b::Vector{T}, + cones::Vector{SupportedCone} +) where {T} + + A_new, b_new = reduce_A_b(presolver,A,b) + cones_new = reduce_cones(presolver, cones) + + return A_new, b_new, cones_new + +end + +function reduce_A_b( + presolver::Presolver{T}, + A::SparseMatrixCSC{T}, + b::Vector{T} +) where{T} + + @assert !isnothing(presolver.reduce_map) + map = presolver.reduce_map + A = A[map.keep_logical,:] + b = b[map.keep_logical] + + A, b + +end + +function reduce_cones( + presolver::Presolver{T}, + cones::Vector{SupportedCone}, +) where {T} + + @assert !isnothing(presolver.reduce_map) + map = presolver.reduce_map + + cones_new = sizehint!(SupportedCone[],length(cones)) + keep_iter = Iterators.Stateful(map.keep_logical) + + for cone in cones + numel_cone = nvars(cone) + + if isa(cone, NonnegativeConeT) + nkeep = count(Iterators.take(keep_iter,numel_cone)) + if nkeep > 0 + push!(cones_new, NonnegativeConeT(nkeep)) + end + else + keep_iter = Iterators.drop(keep_iter,numel_cone) + push!(cones_new, deepcopy(cone)) + end + end + + return cones_new + +end + +function reverse_presolve!( + presolver::Presolver{T}, + solution::DefaultSolution{T}, + variables::DefaultVariables{T} +) where {T} + + @. solution.x = variables.x + + map = presolver.reduce_map + ctr = 1 + + for (idx,keep) in enumerate(map.keep_logical) + if keep + solution.s[idx] = variables.s[ctr] + solution.z[idx] = variables.z[ctr] + ctr += 1 + else + solution.s[idx] = T(presolver.infbound) + solution.z[idx] = zero(T) + end + end + +end + + +function make_reduction_map( + cones::Vector{SupportedCone}, b::Vector{T}, - infbound::T) where {T} + infbound::T +) where {T} keep_logical = trues(length(b)) mreduced = length(b) - # we loop through b and remove any entries that are both infinite - # and in a nonnegative cone + # only try to reduce nn cones. Make a slight contraction + # so that we are firmly "less than" here + infbound *= (1-10*eps(T)) - is_reduced = false - bptr = 1 # index into the b vector - - for (cidx,cone) in enumerate(cone_specs) + idx = 1 # index into the b vector + for cone in cones numel_cone = nvars(cone) - # only try to reduce nn cones. Make a slight contraction - # so that we are firmly "less than" here - infbound *= (1-10*eps(T)) - if isa(cone, NonnegativeConeT) - num_finite = 0 - for i in bptr:(bptr + numel_cone - 1) - if b[i] < infbound - num_finite += 1 - else - keep_logical[i] = false + for _ in 1:numel_cone + if b[idx] > infbound + keep_logical[idx] = false mreduced -= 1 end + idx += 1 end - if num_finite < numel_cone - # contract the cone to a smaller size - cone_specs[cidx] = NonnegativeConeT(num_finite) - is_reduced = true - end + else + # skip this cone + idx += numel_cone end - - bptr += numel_cone end outoption = let - if is_reduced - keep_index = findall(keep_logical) - PresolverRowReductionIndex(keep_logical, keep_index) + if mreduced < length(b) + PresolverRowReductionIndex(keep_logical) else nothing end diff --git a/src/problemdata.jl b/src/problemdata.jl index c5ccb6fc..fa5fb8a6 100644 --- a/src/problemdata.jl +++ b/src/problemdata.jl @@ -1,5 +1,74 @@ import LinearAlgebra +function DefaultProblemData{T}( + P::AbstractMatrix{T}, + q::AbstractVector{T}, + A::AbstractMatrix{T}, + b::AbstractVector{T}, + cones::Vector{SupportedCone}, + settings::Settings{T} +) where {T} + + # some caution is required to ensure we take a minimal, + # but nonzero, number of data copies during presolve steps + P_new, q_new, A_new, b_new, cones_new = + (nothing,nothing,nothing,nothing,nothing) + + # make a copy if badly formatted. istriu is very fast + if !istriu(P) + P_new = triu(P) #copies + P = P_new #rust borrow + end + + # presolve : return nothing if disabled or no reduction + # -------------------------------------- + presolver = try_presolver(A,b,cones,settings) + + if !isnothing(presolver) + (A_new, b_new, cones_new) = presolve(presolver, A, b, cones) + (A, b, cones) = (A_new, b_new, cones_new) + end + + # chordal decomposition : return nothing if disabled or no decomp + # -------------------------------------- + chordal_info = try_chordal_info(A,b,cones,settings) + + if !isnothing(chordal_info) + (P_new, q_new, A_new, b_new, cones_new) = + decomp_augment!(chordal_info, P, q, A, b, settings) + end + + # now make sure we have a clean copy of everything if we + #haven't made one already. Necessary since we will scale + # scale the internal copy and don't want to step on the user + copy_if_nothing(x,y) = isnothing(x) ? deepcopy(y) : x + P_new = copy_if_nothing(P_new,P) + q_new = copy_if_nothing(q_new,q) + A_new = copy_if_nothing(A_new,A) + b_new = copy_if_nothing(b_new,b) + cones_new = copy_if_nothing(cones_new,cones) + + #cap entries in b at INFINITY. This is important + #for inf values that were not in a reduced cone + #this is not considered part of the "presolve", so + #can always happen regardless of user settings + @. b_new .= min(b_new,T(Clarabel.get_infinity())) + + #this ensures m is the *reduced* size m + (m,n) = size(A_new) + + equilibration = DefaultEquilibration{T}(n,m) + + normq = norm(q_new, Inf) + normb = norm(b_new, Inf) + + DefaultProblemData{T}( + P_new,q_new,A_new,b_new,cones_new, + n,m,equilibration,normq,normb, + presolver,chordal_info) + +end + function data_get_normq!(data::DefaultProblemData{T}) where {T} if isnothing(data.normq) @@ -125,7 +194,7 @@ function scale_data!( A::AbstractMatrix{T}, q::AbstractVector{T}, b::AbstractVector{T}, - d::Union{Nothing,AbstractVector{T}}, + d::Option{AbstractVector{T}}, e::AbstractVector{T} ) where {T <: AbstractFloat} @@ -141,3 +210,50 @@ function scale_data!( return nothing end + +function try_chordal_info( + A::SparseMatrixCSC{T}, + b::Vector{T}, + cones::Vector{SupportedCone}, + settings::Settings{T} +) where {T} + + if !settings.chordal_decomposition_enable + return nothing + end + + # nothing to do if there are no PSD cones + if !any(c -> isa(c,PSDTriangleConeT), cones) + return nothing + end + + chordal_info = ChordalInfo(A, b, cones, settings) + + # no decomposition possible + if !is_decomposed(chordal_info) + return nothing + end + + return chordal_info +end + + +function try_presolver( + A::AbstractMatrix{T}, + b::Vector{T}, + cones::Vector{SupportedCone}, + settings::Settings{T} +) where {T} + + if(!settings.presolve_enable) + return nothing + end + + presolver = Presolver{T}(A,b,cones,settings) + + if !is_reduced(presolver) + return nothing + end + + return presolver +end \ No newline at end of file diff --git a/src/settings.jl b/src/settings.jl index ed530776..b98721b9 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -30,7 +30,6 @@ reduced\\_tol\\_infeas_rel | 5e-5 | reduced relative infeasibility reduced\\_tol\\_ktratio | 1e-4 | reduced κ/τ tolerance || __Data Equilibration Settings__|| -|| equilibrate\\_enable | true | enable data equilibration pre-scaling equilibrate\\_max\\_iter | 10 | maximum equilibration scaling iterations equilibrate\\_min\\_scaling | 1e-5 | minimum equilibration scaling allowed @@ -42,7 +41,6 @@ min\\_switch\\_step\\_length | 1e-1 | minimum step size allowed min\\_terminate\\_step\\_length | 1e-4 | minimum step size allowed for symmetric cones && asymmetric cones with Dual scaling || __Linear Solver Settings__|| -|| direct\\_kkt\\_solver | true | use a direct linear solver method (required true) direct\\_solve\\_method | :qdldl | direct linear solver (:qdldl, :mkl or :cholmod) static\\_regularization\\_enable | true | enable KKT static regularization @@ -56,8 +54,15 @@ iterative\\_refinement\\_reltol | 1e-12 | iterative refinement relat iterative\\_refinement\\_abstol | 1e-12 | iterative refinement absolute tolerance iterative\\_refinement\\_max\\_iter | 10 | iterative refinement maximum iterations iterative\\_refinement\\_stop\\_ratio | 5.0 | iterative refinement stalling tolerance +|| __Preprocessing Settings presolve_enable | true | enable presolve constraint reduction +|| +__Chordal Decomposition Settings__|| +chordal\\_decomposition\\_enable | true | enable chordal decomposition +chordal\\_decomposition\\_merge_method | :clique_graph | chordal decomposition merge method (:none, :parent_child or :clique_graph) +chordal\\_decomposition\\_compact | true | assemble decomposed system in "compact" form +chordal\\_decomposition\\_complete_dual | false | complete PSD dual variables after decomposition """ Base.@kwdef mutable struct Settings{T <: AbstractFloat} @@ -120,6 +125,12 @@ Base.@kwdef mutable struct Settings{T <: AbstractFloat} #preprocessing presolve_enable::Bool = true + #chordal decomposition + chordal_decomposition_enable::Bool = true + chordal_decomposition_merge_method::Symbol = :clique_graph + chordal_decomposition_compact::Bool = true + chordal_decomposition_complete_dual::Bool = false + end Settings(args...) = Settings{DefaultFloat}(args...) @@ -183,9 +194,9 @@ function Base.show(io::IO, settings::Clarabel.Settings{T}) where {T} end println(io) # and the settings - for row in 1:size(table,1) + for row in axes(table,1) @printf(io, " ") - for col in 1:size(table,2) + for col in axes(table,2) @printf(io, " %s ", table[row,col]) end println(io) diff --git a/src/solution.jl b/src/solution.jl index b471f745..90879ae4 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -1,5 +1,5 @@ -function solution_finalize!( +function solution_post_process!( solution::DefaultSolution{T}, data::DefaultProblemData{T}, variables::DefaultVariables{T}, @@ -7,54 +7,52 @@ function solution_finalize!( settings::Settings{T} ) where {T} - solution.status = info.status - solution.obj_val = info.cost_primal - solution.obj_val_dual = info.cost_dual + solution.status = info.status + is_infeasible = status_is_infeasible(info.status) - #if we have an infeasible problem, normalize - #using κ to get an infeasibility certificate. - #Otherwise use τ to get a solution. - if status_is_infeasible(info.status) - scaleinv = one(T) / variables.κ - solution.obj_val = NaN + if is_infeasible + solution.obj_val = NaN solution.obj_val_dual = NaN - else - scaleinv = one(T) / variables.τ - end - - #also undo the equilibration - d = data.equilibration.d; dinv = data.equilibration.dinv - e = data.equilibration.e; einv = data.equilibration.einv - cscale = data.equilibration.c[] + else + solution.obj_val = info.cost_primal + solution.obj_val_dual = info.cost_dual + end - @. solution.x = variables.x * d * scaleinv + solution.iterations = info.iterations + solution.r_prim = info.res_primal + solution.r_dual = info.res_dual - map = data.presolver.reduce_map - if !isnothing(map) - map = data.presolver.reduce_map - @. solution.z[map.keep_index] = variables.z * e * (scaleinv / cscale) - @. solution.s[map.keep_index] = variables.s * einv * scaleinv + # unscale the variables to get a solution + # to the internal problem as we solved it + variables_unscale!(variables,data,is_infeasible) - #eliminated constraints get huge slacks - #and are assumed to be nonbinding - @. solution.s[!map.keep_logical] = T(data.presolver.infbound) - @. solution.z[!map.keep_logical] = zero(T) + # unwind the chordal decomp and presolve, in the + # reverse of the order in which they were applied + if !isnothing(data.chordal_info) + variables = decomp_reverse!( + data.chordal_info, variables, data.cones, settings) + end + if !isnothing(data.presolver) + reverse_presolve!(data.presolver, solution, variables) else - @. solution.z = variables.z * e * (scaleinv / cscale) - @. solution.s = variables.s * einv * scaleinv + @. solution.x = variables.x + @. solution.z = variables.z + @. solution.s = variables.s end - - - solution.iterations = info.iterations - solution.solve_time = info.solve_time - solution.r_prim = info.res_primal - solution.r_dual = info.res_dual return nothing end +function solution_finalize!( + solution::DefaultSolution{T}, + info::DefaultInfo{T}, +) where {T} + + solution.solve_time = info.solve_time + +end function Base.show(io::IO, solution::DefaultSolution) diff --git a/src/solver.jl b/src/solver.jl index 6fa057ec..43a326ad 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -82,6 +82,10 @@ function setup!( cones::Vector{<:SupportedCone}, ) where{T} + # project against cones with overly specific type, e.g. + # when all of the cones are NonnegativeConeT + cones = convert(Vector{SupportedCone},cones) + #sanity check problem dimensions _check_dimensions(P,q,A,b,cones) @@ -90,16 +94,19 @@ function setup!( @timeit s.timers "setup!" begin - #reduce the cone sizes. (A,b) will be reduced - #within the problem data constructor. Also makes - #an internal copy of the user cone specification - presolver = Presolver{T}(A,b,cones,s.settings) + # user facing results go here + s.solution = DefaultSolution{T}(A.n,A.m) + + # presolve / chordal decomposition if needed, + # then take an internal copy of the problem data + @timeit s.timers "presolve" begin + s.data = DefaultProblemData{T}(P,q,A,b,cones,s.settings) + end - s.cones = CompositeCone{T}(presolver.cone_specs) - s.data = DefaultProblemData{T}(P,q,A,b,s.cones,presolver) + s.cones = CompositeCone{T}(s.data.cones) s.data.m == s.cones.numel || throw(DimensionMismatch()) - s.variables = DefaultVariables{T}(s.data.n,s.cones) + s.variables = DefaultVariables{T}(s.data.n,s.data.m) s.residuals = DefaultResiduals{T}(s.data.n,s.data.m) #equilibrate problem data immediately on setup. @@ -114,14 +121,11 @@ function setup!( end # work variables for assembling step direction LHS/RHS - s.step_rhs = DefaultVariables{T}(s.data.n,s.cones) - s.step_lhs = DefaultVariables{T}(s.data.n,s.cones) + s.step_rhs = DefaultVariables{T}(s.data.n,s.data.m) + s.step_lhs = DefaultVariables{T}(s.data.n,s.data.m) # a saved copy of the previous iterate - s.prev_vars = DefaultVariables{T}(s.data.n,s.cones) - - # user facing results go here - s.solution = DefaultSolution{T}(s.data.presolver.mfull,s.data.n) + s.prev_vars = DefaultVariables{T}(s.data.n,s.data.m) end @@ -231,7 +235,10 @@ function solve!( #update the scalings #-------------- + @timeit s.timers "scale cones" begin is_scaling_success = variables_scale_cones!(s.variables,s.cones,μ,scaling) + end + # check whether variables are interior points (action,scaling) = _strategy_checkpoint_is_scaling_success(s,is_scaling_success,scaling) if action === Fail; break; @@ -336,8 +343,16 @@ function solve!( @notimeit info_print_status(s.info,s.settings) end - info_finalize!(s.info,s.residuals,s.settings,s.timers) #halts timers - solution_finalize!(s.solution,s.data,s.variables,s.info,s.settings) + @timeit s.timers "post-process" begin + #check for "almost" convergence checks and then extract solution + info_post_process!(s.info,s.residuals,s.settings) + solution_post_process!(s.solution,s.data,s.variables,s.info,s.settings) + end + + #halt timers + info_finalize!(s.info,s.timers) + solution_finalize!(s.solution,s.info) + @notimeit info_print_footer(s.info,s.settings) diff --git a/src/types.jl b/src/types.jl index e27c3d17..f7b716a4 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,21 +1,41 @@ using TimerOutputs # ------------------------------------- -# abstract type defs +# default solver component types # ------------------------------------- -abstract type AbstractVariables{T <: AbstractFloat} end -abstract type AbstractEquilibration{T <: AbstractFloat} end -abstract type AbstractResiduals{T <: AbstractFloat} end -abstract type AbstractProblemData{T <: AbstractFloat} end -abstract type AbstractKKTSystem{T <: AbstractFloat} end -abstract type AbstractKKTSolver{T <: AbstractFloat} end -abstract type AbstractInfo{T <: AbstractFloat} end -abstract type AbstractSolution{T <: AbstractFloat} end -abstract type AbstractSolver{T <: AbstractFloat} end -# ------------------------------------- -# default solver subcomponent implementations -# ------------------------------------- +# --------------------- +# presolver and internals +# --------------------- + +struct PresolverRowReductionIndex + + # vector of length = original RHS. Entries are false + # for those rows that should be eliminated before solve + keep_logical::Vector{Bool} + +end +struct Presolver{T} + + # original cones of the problem + init_cones::Vector{SupportedCone} + + # record of reduced constraints for NN cones with inf bounds + reduce_map::Option{PresolverRowReductionIndex} + + # size of original and reduced RHS, respectively + mfull::Int64 + mreduced::Int64 + + # inf bound that was taken from the module level + # and should be applied throughout. Held here so + # that any subsequent change to the module's state + # won't mess up our solver mid-solve + infbound::Float64 + +end + +Presolver(args...) = Presolver{DefaultFloat}(args...) # --------------- # variables @@ -24,19 +44,19 @@ abstract type AbstractSolver{T <: AbstractFloat} end mutable struct DefaultVariables{T} <: AbstractVariables{T} x::Vector{T} - s::ConicVector{T} - z::ConicVector{T} + s::Vector{T} + z::Vector{T} τ::T κ::T function DefaultVariables{T}( n::Integer, - cones::CompositeCone{T} + m::Integer, ) where {T} - x = Vector{T}(undef,n) - s = ConicVector{T}(cones) - z = ConicVector{T}(cones) + x = zeros(T,n) + s = zeros(T,m) + z = zeros(T,m) τ = T(1) κ = T(1) @@ -84,14 +104,14 @@ mutable struct DefaultResiduals{T} <: AbstractResiduals{T} function DefaultResiduals{T}(n::Integer, m::Integer) where {T} - rx = Vector{T}(undef,n) - rz = Vector{T}(undef,m) + rx = zeros(T,n) + rz = zeros(T,m) rτ = T(1) - rx_inf = Vector{T}(undef,n) - rz_inf = Vector{T}(undef,m) + rx_inf = zeros(T,n) + rz_inf = zeros(T,m) - Px = Vector{T}(undef,n) + Px = zeros(T,n) new(rx,rz,rτ,rx_inf,rz_inf,zero(T),zero(T),zero(T),zero(T),Px) end @@ -112,22 +132,22 @@ struct DefaultEquilibration{T} <: AbstractEquilibration{T} #to be treated as diagonal scaling data d::Vector{T} dinv::Vector{T} - e::ConicVector{T} - einv::ConicVector{T} + e::Vector{T} + einv::Vector{T} #overall scaling for objective function c::Base.RefValue{T} function DefaultEquilibration{T}( - nvars::Int, - cones::CompositeCone{T}, + n::Int64, + m::Int64, ) where {T} #Left/Right diagonal scaling for problem data - d = ones(T,nvars) - dinv = ones(T,nvars) - e = ConicVector{T}(cones); e .= one(T) - einv = ConicVector{T}(cones); einv .= one(T) + d = ones(T,n) + dinv = ones(T,n) + e = ones(T,m) + einv = ones(T,m) c = Ref(T(1.)) @@ -149,6 +169,7 @@ mutable struct DefaultProblemData{T} <: AbstractProblemData{T} q::Vector{T} A::AbstractMatrix{T} b::Vector{T} + cones::Vector{SupportedCone} n::DefaultInt m::DefaultInt equilibration::DefaultEquilibration{T} @@ -157,55 +178,11 @@ mutable struct DefaultProblemData{T} <: AbstractProblemData{T} # during data updating to allow for multiple updates, and # then recalculated during solve if needed - normq::Union{T,Nothing} #unscaled inf norm of q - normb::Union{T,Nothing} #unscaled inf norm of b - - presolver::Presolver{T} - - function DefaultProblemData{T}( - P::AbstractMatrix{T}, - q::AbstractVector{T}, - A::AbstractMatrix{T}, - b::AbstractVector{T}, - cones::CompositeCone{T}, - presolver::Presolver{T}, - ) where {T} - - # dimension checks will have already been - # performed during problem setup, so skip here - - #take an internal copy of all problem - #data, since we are going to scale it - P = triu(P) - q = deepcopy(q) - - (A,b) = let - map = presolver.reduce_map; - if !isnothing(map) - ( - A[map.keep_logical,:], - b[map.keep_logical] - ) - else - (deepcopy(A),deepcopy(b)) - end - end + normq::Option{T} #unscaled inf norm of q + normb::Option{T} #unscaled inf norm of b - #cap entries in b at INFINITY. This is important - #for inf values that were not in a reduced cone - @. b .= min(b,T(presolver.infbound)) - - #this ensures m is the *reduced* size m - (m,n) = size(A) - - equilibration = DefaultEquilibration{T}(n,cones) - - normq = norm(q, Inf) - normb = norm(b, Inf) - - new(P,q,A,b,n,m,equilibration,normq,normb,presolver) - - end + presolver::Option{Presolver{T}} + chordal_info::Option{ChordalInfo{T}} end @@ -295,11 +272,11 @@ mutable struct DefaultSolution{T} <: AbstractSolution{T} end -function DefaultSolution{T}(m,n) where {T <: AbstractFloat} +function DefaultSolution{T}(n,m) where {T <: AbstractFloat} - x = Vector{T}(undef,n) - z = Vector{T}(undef,m) - s = Vector{T}(undef,m) + x = zeros(T,n) + z = zeros(T,m) + s = zeros(T,m) # seemingly reasonable defaults status = UNSOLVED @@ -329,16 +306,16 @@ Initializes an empty Clarabel solver that can be filled with problem data using: """ mutable struct Solver{T <: AbstractFloat} <: AbstractSolver{T} - data::Union{AbstractProblemData{T},Nothing} - variables::Union{AbstractVariables{T},Nothing} - cones::Union{CompositeCone{T},Nothing} - residuals::Union{AbstractResiduals{T},Nothing} - kktsystem::Union{AbstractKKTSystem{T},Nothing} - info::Union{AbstractInfo{T},Nothing} - step_lhs::Union{AbstractVariables{T},Nothing} - step_rhs::Union{AbstractVariables{T},Nothing} - prev_vars::Union{AbstractVariables{T},Nothing} - solution::Union{AbstractSolution{T},Nothing} + data::Option{AbstractProblemData{T}} + variables::Option{AbstractVariables{T}} + cones::Option{CompositeCone{T}} + residuals::Option{AbstractResiduals{T}} + kktsystem::Option{AbstractKKTSystem{T}} + info::Option{AbstractInfo{T}} + step_lhs::Option{AbstractVariables{T}} + step_rhs::Option{AbstractVariables{T}} + prev_vars::Option{AbstractVariables{T}} + solution::Option{AbstractSolution{T}} settings::Settings{T} timers::TimerOutput @@ -370,3 +347,4 @@ function Solver{T}(d::Dict) where {T} end Solver(args...; kwargs...) = Solver{DefaultFloat}(args...; kwargs...) + diff --git a/src/utils/mathutils.jl b/src/utils/mathutils.jl index 564fb09b..b9c061d1 100644 --- a/src/utils/mathutils.jl +++ b/src/utils/mathutils.jl @@ -329,6 +329,27 @@ function triangular_index(k::Int) triangular_number(k) end +# given an index into the upper triangular part of a matrix, return +# its row and column position + +function upper_triangular_index_to_coord(linearidx::Int) + col = (isqrt(8 * linearidx) + 1) >> 1 + row = linearidx - triangular_number(col - 1) + (row,col) +end + +# given a row and column position, return the index into the upper +# triangular part of the matrix + +function coord_to_upper_triangular_index(coord::Tuple{Int, Int}) + (i,j) = coord + if i <= j + return triangular_number(j-1) + i + else + return triangular_number(i-1) + j + end +end + # Put elements from a vector of length n*(n+1)/2 into the upper triangle # of an nxn matrix. It does NOT perform any scaling of the vector entries. diff --git a/src/variables.jl b/src/variables.jl index 60d0284b..53cfb735 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -231,9 +231,51 @@ function variables_rescale!(variables) invscale = 1/scale; variables.x .*= invscale - variables.z.vec .*= invscale - variables.s.vec .*= invscale - variables.τ *= invscale - variables.κ *= invscale + variables.z .*= invscale + variables.s .*= invscale + variables.τ *= invscale + variables.κ *= invscale +end + + +# ---------------------- +# remaining functions are specific to DefaultVariables types +# only and not part of the general AbstractVariables interface + + +function variables_unscale!( + variables::DefaultVariables{T}, + data::DefaultProblemData{T}, + is_infeasible::Bool +) where {T} + + #if we have an infeasible problem, normalize + #using κ to get an infeasibility certificate. + #Otherwise use τ to get an unscaled solution. + if is_infeasible + scaleinv = one(T) / variables.κ + else + scaleinv = one(T) / variables.τ + end + + #also undo the equilibration + d = data.equilibration.d + e = data.equilibration.e + einv = data.equilibration.einv + cscale = data.equilibration.c[] + + @. variables.x *= d * scaleinv + @. variables.z *= e * (scaleinv / cscale) + @. variables.s *= einv * scaleinv + + variables.τ *= scaleinv + variables.κ *= scaleinv + +end + + +# return (n_variables, n_duals) +function variables_dims(variables::DefaultVariables{T}) where {T} + return (length(variables.x), length(variables.s)) end \ No newline at end of file diff --git a/test/OptTests/data_updating.jl b/test/OptTests/data_updating.jl index d347d6a7..a73a6b60 100644 --- a/test/OptTests/data_updating.jl +++ b/test/OptTests/data_updating.jl @@ -17,6 +17,7 @@ function updating_test_data(Type::Type{T}) where {T <: AbstractFloat} settings = Clarabel.Settings{T}() settings.presolve_enable = false + settings.chordal_decomposition_enable = false return (P,q,A,b,cones,settings) end diff --git a/test/UnitTests/test_coneops_psdtrianglecone.jl b/test/UnitTests/test_coneops_psdtrianglecone.jl index ae9ce183..f30858fe 100644 --- a/test/UnitTests/test_coneops_psdtrianglecone.jl +++ b/test/UnitTests/test_coneops_psdtrianglecone.jl @@ -28,12 +28,12 @@ FloatT = Float64 y = zeros(Clarabel.triangular_number(n)) # check inner product identity - map((v,M)->Clarabel._mat_to_svec!(v,M), (x,y), (X,Y)) + map((v,M)->Clarabel.mat_to_svec!(v,M), (x,y), (X,Y)) @test x'y - tr(X'Y)≈ 0 atol = 1e-12 # check round trip - Clarabel._mat_to_svec!(x,X) - Clarabel._svec_to_mat!(Z,x) + Clarabel.mat_to_svec!(x,X) + Clarabel.svec_to_mat!(Z,x) @test norm(X-Z) ≈ 0 atol = 1e-12 @@ -48,13 +48,13 @@ FloatT = Float64 Y = randsym(rng, n) Z = randsym(rng, n) (x,y,z) = map(v->zeros(Clarabel.triangular_number(n)),1:3) - map((v,M)->Clarabel._mat_to_svec!(v,M), (y,z), (Y,Z)) + map((v,M)->Clarabel.mat_to_svec!(v,M), (y,z), (Y,Z)) K = Clarabel.PSDTriangleCone(n) X1 .= 0.5*(Y*Z + Z*Y) Clarabel.circ_op!(K,x,y,z) - Clarabel._svec_to_mat!(X2,x) + Clarabel.svec_to_mat!(X2,x) @test tr(X1) ≈ dot(y,z) @test norm(X2-X1) ≈ 0 atol = 1e-12 @@ -72,7 +72,7 @@ FloatT = Float64 Z = zeros(n,n) W = zeros(n,n) (x,z,λ,w) = map(m->zeros(K.numel), 1:4) - map((v,M)->Clarabel._mat_to_svec!(v,M),(x,z,λ,w),(X,Z,Λ,W)) + map((v,M)->Clarabel.mat_to_svec!(v,M),(x,z,λ,w),(X,Z,Λ,W)) #Z = 1/2(ΛX + XΛ) Clarabel.circ_op!(K,z,λ,x) @@ -80,7 +80,7 @@ FloatT = Float64 #W should now be the solution to 1/2(ΛW + WΛ) = Z K.data.λ .= λdiag #diagonal internal scaling Clarabel.λ_inv_circ_op!(K,w,z) - Clarabel._svec_to_mat!(W,w) + Clarabel.svec_to_mat!(W,w) #now we should have x = w @test norm(x - w) ≈ 0 atol = 100*eps(FloatT) @@ -101,9 +101,9 @@ FloatT = Float64 x = zeros(K.numel) XplusaI = X + a*I(n) - Clarabel._mat_to_svec!(x,X) + Clarabel.mat_to_svec!(x,X) Clarabel.scaled_unit_shift!(K,x,a,Clarabel.PrimalCone) - Clarabel._svec_to_mat!(X,x) + Clarabel.svec_to_mat!(X,x) @test norm(X - XplusaI) ≈ 0 @@ -120,7 +120,7 @@ FloatT = Float64 Z = randpsd(rng,n) (s,z) = map(m->zeros(K.numel), 1:2) - map((v,M)->Clarabel._mat_to_svec!(v,M),(s,z),(S,Z)) + map((v,M)->Clarabel.mat_to_svec!(v,M),(s,z),(S,Z)) μ = 0.0 #placeholder value, not used strategy = Clarabel.PrimalDual @@ -149,7 +149,7 @@ FloatT = Float64 S = randpsd(rng,n); dS = randsym(rng,n) (s,z,ds,dz) = map(m->zeros(K.numel), 1:4) - map((v,M)->Clarabel._mat_to_svec!(v,M),(s,z,ds,dz),(S,Z,dS,dZ)) + map((v,M)->Clarabel.mat_to_svec!(v,M),(s,z,ds,dz),(S,Z,dS,dZ)) #compute internal scaling required for step calc μ = 0.0 #placeholder value, not used @@ -172,7 +172,7 @@ FloatT = Float64 #should reach maximum step dS .= randpsd(rng,n); dZ .= randpsd(rng,n) - map((v,M)->Clarabel._mat_to_svec!(v,M),(ds,dz),(dS,dZ)) + map((v,M)->Clarabel.mat_to_svec!(v,M),(ds,dz),(dS,dZ)) (αz,αs) = Clarabel.step_length(K,dz,ds,z,s,settings,1.0) @test min(αz,αs) ≈ 1.0 rtol = 10*eps(FloatT) @@ -185,7 +185,7 @@ FloatT = Float64 (Z,S,V1,V2) = map(m->randpsd(rng,n), 1:4) (s,z,v1,v2) = map(m->zeros(K.numel), 1:4) - map((v,M)->Clarabel._mat_to_svec!(v,M),(s,z,v1,v2),(S,Z,V1,V2)) + map((v,M)->Clarabel.mat_to_svec!(v,M),(s,z,v1,v2),(S,Z,V1,V2)) #compute internal scaling required for step calc μ = 0.0 #placeholder value, not used @@ -218,7 +218,7 @@ FloatT = Float64 (Z,S,V1,V2,V3) = map(m->randpsd(rng,n), 1:5) (z,s,v1,v2,v3) = map(m->zeros(K.numel), 1:5) - map((v,M)->Clarabel._mat_to_svec!(v,M),(s,z,v1,v2,v3),(S,Z,V1,V2,V3)) + map((v,M)->Clarabel.mat_to_svec!(v,M),(s,z,v1,v2,v3),(S,Z,V1,V2,V3)) #compute internal scaling required for step calc μ = 0.0 #placeholder value, not used