Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor + small improvements in order reduction functions #155

Merged
merged 2 commits into from
May 5, 2020
Merged

Refactor + small improvements in order reduction functions #155

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
85be5b3
refactor + update order reduction functions
mforets May 2, 2020
d424fc5
reorder args
mforets May 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
157 changes: 81 additions & 76 deletions src/Flowpipes/setops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -433,59 +433,94 @@ end

abstract type AbstractReductionMethod end

# These methods split the given zonotope Z into two zonotopes, K and L, where
# K contains the most "representative" generators and L contains the generators
# that are reduced, Lred using a box overapproximation
struct GIR05 <: AbstractReductionMethod end
struct COMB03 <: AbstractReductionMethod end

const _COMB03 = COMB03()
const _GIR05 = GIR05()

# algorithm selection
_reduce_order(Z::Zonotope, r::Number) = _reduce_order_GIR05(Z, r) # default
_reduce_order(Z::Zonotope, r::Number, ::GIR05) = _reduce_order_GIR05(Z, r)
_reduce_order(Z::Zonotope, r::Number, ::COMB03) = _reduce_order_COMB03(Z, r)

# Implements zonotope order reduction method from [COMB03]
# We follow the notation from [YS18]
function _reduce_order_COMB03(Z::Zonotope{N}, r::Number) where {N}
r >= 1 || throw(ArgumentError("the target order should be at least 1, but it is $r"))
c = Z.center
G = Z.generators
n, p = size(G)

# r is bigger than the order of Z => don't reduce
(r * n >= p) && return Z

# split Z into K and L, where K contains the most "representative" generators
# and L contains the generators that are reduced, Lred

# this algorithm sort generators by decreasing 2-norm
weights = zeros(N, p)
@inbounds for i in 1:p
weights[i] = norm(view(G, :, i), 2)
# return the indices of the generators in G (= columns) sorted according to the COMB03 method
# the generator index with highest score goes first
function _weighted_gens!(indices, G::AbstractMatrix{N}, ::COMB03) where {N}
p = size(G, 2)
weights = Vector{N}(undef, p)
@inbounds for j in 1:p
v = view(G, :, j)
weights[j] = norm(v, 2)
end
indices = Vector{Int}(undef, p)
sortperm!(indices, weights, rev=true, initialized=false)
return indices
end

# the first m generators with greatest weight
m = floor(Int, n * (r - 1))
# return the indices of the generators in G (= columns) sorted according to the GIR05 method
# the generator index with highest score goes first
function _weighted_gens!(indices, G::AbstractMatrix{N}, ::GIR05) where {N}
n, p = size(G)
weights = Vector{N}(undef, p)
@inbounds for j in 1:p
aux_norm_1 = zero(N)
aux_norm_inf = zero(N)
for i in 1:n
abs_Gij = abs(G[i, j])
aux_norm_1 += abs_Gij
if aux_norm_inf < abs_Gij
aux_norm_inf = abs_Gij
end
end
weights[j] = aux_norm_1 - aux_norm_inf
end
sortperm!(indices, weights, rev=true, initialized=false)
return indices
end

# compute interval hull of L
# compute interval hull of the generators of G (= columns) corresponding to `indices`
function _interval_hull(G::AbstractMatrix{N}, indices) where {N}
n, p = size(G)
Lred = zeros(N, n, n)
@inbounds for i in 1:n
for j in indices[m+1:end]
for j in indices
Lred[i, i] += abs(G[i, j])
end
end
return Lred
end

if isone(r)
return Zonotope(c, Lred)
# implementation for static arrays
function _interval_hull(G::SMatrix{n, p, N, L}, indices) where {n, p, N, L}
Lred = zeros(MMatrix{n, n, N})
@inbounds for i in 1:n
for j in indices
Lred[i, i] += abs(G[i, j])
end
end
return SMatrix{n, n}(Lred)
end

K = G[:, indices[1:m]]
Gred = hcat(K, Lred)
return Zonotope(c, Gred)
# given an n x p matrix G and a vector of m integer indices with m <= p,
# concatenate the columns of G given by `indices` with the matrix Lred
function _hcat_KLred(G::AbstractMatrix, indices, Lred::AbstractMatrix)
K = view(G, :, indices)
return hcat(K, Lred)
end

# Implements zonotope order reduction method from [GIR05]
# implementation for static arrays
function _hcat_KLred(G::SMatrix{n, p, N, L1}, indices, Lred::SMatrix{n, n, N, L2}) where {n, p, N, L1, L2}
m = length(indices)
K = SMatrix{n, m}(view(G, :, indices))
return hcat(K, Lred)
end

# Implements zonotope order reduction method from [COMB03]
# We follow the notation from [YS18]
function _reduce_order_GIR05(Z::Zonotope{N}, r::Number) where {N}
function _reduce_order_COMB03(Z::Zonotope{N}, r::Number) where {N}
r >= 1 || throw(ArgumentError("the target order should be at least 1, but it is $r"))
c = Z.center
G = Z.generators
Expand All @@ -494,76 +529,46 @@ function _reduce_order_GIR05(Z::Zonotope{N}, r::Number) where {N}
# r is bigger than the order of Z => don't reduce
(r * n >= p) && return Z

# split Z into K and L, where K contains the most "representative" generators
# and L contains the generators that are reduced, Lred

# this algorithm sort generators by ||⋅||₁ - ||⋅||∞ difference
weights = zeros(N, p)
@inbounds for i in 1:p
v = view(G, :, i)
weights[i] = norm(v, 1) - norm(v, Inf)
end
# this algorithm sort generators by decreasing 2-norm
indices = Vector{Int}(undef, p)
sortperm!(indices, weights, rev=true, initialized=false)
_weighted_gens!(indices, G, _COMB03)

# the first m generators with greatest weight
# the first m generators have greatest weight
m = floor(Int, n * (r - 1))

# compute interval hull of L
Lred = zeros(N, n, n)
@inbounds for i in 1:n
for j in indices[m+1:end]
Lred[i, i] += abs(G[i, j])
end
end
Lred = _interval_hull(G, view(indices, (m+1):p))

if isone(r)
return Zonotope(c, Lred)
end
isone(r) && return Zonotope(c, Lred)

K = G[:, indices[1:m]]
Gred = hcat(K, Lred)
Gred = _hcat_KLred(G, view(indices, 1:m), Lred)
return Zonotope(c, Gred)
end

# implementation for zonotopes using static arrays
function _reduce_order_GIR05(Z::Zonotope{N, SVector{n, N}, SMatrix{n, p, N, LM}}, r::Number) where {N, n, p, LM}
# Implements zonotope order reduction method from [GIR05]
# We follow the notation from [YS18]
function _reduce_order_GIR05(Z::Zonotope{N}, r::Number) where {N}
r >= 1 || throw(ArgumentError("the target order should be at least 1, but it is $r"))
c = Z.center
G = Z.generators
n, p = size(G)

# r is bigger than the order of Z => don't reduce
(r * n >= p) && return Z

# split Z into K and L, where K contains the most "representative" generators
# and L contains the generators that are reduced, Lred

# this algorithm sort generators by ||⋅||₁ - ||⋅||∞ difference
weights = zeros(N, p)
@inbounds for i in 1:p
v = view(G, :, i)
weights[i] = norm(v, 1) - norm(v, Inf)
end
# this algorithm sorts generators by ||⋅||₁ - ||⋅||∞ difference
indices = Vector{Int}(undef, p)
sortperm!(indices, weights, rev=true, initialized=false)
_weighted_gens!(indices, G, _GIR05)

# the first m generators with greatest weight
# the first m generators have greatest weight
m = floor(Int, n * (r - 1))

# compute interval hull of L
Lred = zeros(MMatrix{n, n, N})
@inbounds for i in 1:n
for j in indices[m+1:end]
Lred[i, i] += abs(G[i, j])
end
end
Lred = _interval_hull(G, view(indices, (m+1):p))

if isone(r)
return Zonotope(c, SMatrix{n, n}(Lred))
end
isone(r) && return Zonotope(c, Lred)

K = SMatrix{n, m}(view(G, :, indices[1:m]))
Gred = hcat(K, Lred)
Gred = _hcat_KLred(G, view(indices, 1:m), Lred)
return Zonotope(c, Gred)
end

Expand Down