Zonotope order reduction (Girard's algorithm) #211
Conversation
src/Zonotope.jl
Outdated
| This function implements the algorithm described in A. Girard's | ||
| *Reachability of Uncertain Linear Systems Using Zonotopes*, HSCC. Vol. 5. 2005. | ||
| """ | ||
| function reduce_order(Z::Zonotope, r::Int) |
src/Zonotope.jl
Outdated
| function reduce_order(Z::Zonotope, r::Int) | ||
| c, G = Z.center, Z.generators | ||
| p, d = ngens(Z), dim(Z) | ||
| h = zeros(p) |
There was a problem hiding this comment.
zeros(N, p) (you need to define N first), assuming that vecnorm is not always Float64 (I did not check)?
There was a problem hiding this comment.
I checked it: For Rational it returns Float64, but for Float32 it returns Float32. So let us keep it as is (Float32 is automatically promoted to Float64, I checked it).
src/Zonotope.jl
Outdated
| end | ||
| ind = sortperm(h) | ||
|
|
||
| if r * d >= p |
There was a problem hiding this comment.
Pull up as far as possible? No need to compute h in that case!
src/Zonotope.jl
Outdated
| rg = G[:, ind[1:m]] | ||
|
|
||
| # interval hull computation of reduced generators | ||
| d = sum(abs.(rg), 2) |
There was a problem hiding this comment.
What does this line do? d is never used again...
There was a problem hiding this comment.
the sum is used just below, to build Gbox
There was a problem hiding this comment.
Yes, but this line can be removed I think (or use d below instead).
src/Zonotope.jl
Outdated
| Gred = Gbox | ||
| end | ||
| Zred = Zonotope(c, Gred) | ||
| return Zred |
src/Zonotope.jl
Outdated
| d = sum(abs.(rg), 2) | ||
| Gbox = diagm(sum(abs.(rg), 2)[:]) | ||
| if m > p | ||
| Gred = [G[:, ind[m+1]:end] Gbox] |
There was a problem hiding this comment.
Should it be ind[m+1:end] (not sure)? Maybe add a comment what happens here...
There was a problem hiding this comment.
ind[m+1] points to a number (index), so this chooses from that index to the last one
There was a problem hiding this comment.
but ind is a vector of length p thus the check is nonsense (we need m+1 <= p), updated.
There was a problem hiding this comment.
equivalent to m < p for integers ;)
There was a problem hiding this comment.
actually, ind[m+1:end] because these is the set of non reduced generators, well done.
src/Zonotope.jl
Outdated
| p, d = ngens(Z), dim(Z) | ||
| h = zeros(p) | ||
| for i in 1:p | ||
| h[i] = vecnorm(G[:, i], 1) - vecnorm(G[:, i], Inf) |
There was a problem hiding this comment.
For vectors norm and vecnorm should be equivalent. Any reason why you used vecnorm here?
There was a problem hiding this comment.
not particular, no. i tried to use vecnorm column-wise in the matrix G with a broadcast but it didn't work, thus the for loop.
Closes #195.