forked from JuliaPolyhedra/Polyhedra.jl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
aff.jl
228 lines (202 loc) · 7.59 KB
/
aff.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
export HyperPlanesIntersection, LinesHull, affinehull, linespace
# It is easy to go from H-rep of affine space to V-rep of affine space by computing the kernel of a matrix using RowEchelon
# However, it is really worth it since Base.in for an HRepElement in HyperPlanesIntersection and Base.in for an VRepElement in LinesHull are false already.
function remproj(x::RepElement{<:Integer}, l::RepElement{<:Integer})
# generates large numbers but keeps the type integer
x * dot(coord(l), coord(l)) - l * dot(coord(x), coord(l))
end
function remproj(x, l)
ll = dot(coord(l), coord(l))
if isapproxzero(ll)
x
else
simplify(x - l * ((dot(coord(x), coord(l)) / ll)))
end
end
# H-representation
"""
hrep(hyperplanes::HyperPlaneIt; d::FullDim)
Creates an affine space of full dimension `d` from the list of hyperplanes
`hyperplanes`.
### Examples
```julia
hrep([HyperPlane([0, 1, 0], 1), HyperPlane([0, 0, 1], 0)])
```
creates the 1-dimensional affine subspace containing all the points ``(x_1, 0, 0)``, i.e. the ``x_1``-axis.
```julia
hrep([HyperPlane([1, 1], 1), HyperPlane([1, 0], 0)])
```
creates the 0-dimensional affine subspace only containing the point ``(0, 1)``.
"""
function hrep(hyperplanes::HyperPlaneIt; d::FullDim = FullDim_rec(hyperplanes))
return HyperPlanesIntersection(d, hyperplanes)
end
# Representation of an affine space as the intersection of hyperplanes.
# Also called affine set, affine manifold, affine variety, linear variety or flat.
# An affine space L satisfies:
# λx + (1-λ)y ∈ L, ∀x, y ∈ L, ∀ λ ∈ R
# Note that λ is not rhyperplaneuired to be between 0 and 1 as in convex sets.
mutable struct HyperPlanesIntersection{T, AT, D<:FullDim} <: HAffineSpace{T}
d::D
# HyperPlanes whose intersection is the affine space
hyperplanes::Vector{HyperPlane{T, AT}}
orthogonalized::Bool
function HyperPlanesIntersection{T, AT, D}(d::FullDim,
hps::HyperPlaneIt{T}) where {T, AT, D}
new{T, AT, D}(FullDim_convert(D, d), lazy_collect(hps), false)
end
end
function HyperPlanesIntersection(d::FullDim,
it::ElemIt{HyperPlane{T, AT}}) where {T, AT}
return HyperPlanesIntersection{T, AT, typeof(d)}(d, it)
end
HyperPlanesIntersection{T, AT}(d::FullDim) where {T, AT} = HyperPlanesIntersection(d, HyperPlane{T, AT}[])
FullDim(h::HyperPlanesIntersection) = h.d
hvectortype(L::Type{<:HyperPlanesIntersection{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:HyperPlanesIntersection}, d::FullDim, ::Type{T}) where T = HyperPlanesIntersection{T, similar_type(hvectortype(PT), d, T), typeof(d)}
function Base.empty!(L::HyperPlanesIntersection{T, AT}) where {T, AT}
# We want the dimension to be `-1` and the hyperplanes to make
# an infeasible system so we build the system `x = 0 & 0 = 1`.
empty!(L.hyperplanes)
for i in 1:fulldim(L)
push!(L.hyperplanes, HyperPlane(basis(AT, L.d, i), zero(T)))
end
push!(L.hyperplanes, HyperPlane(origin(AT, fulldim(L)), one(T)))
return L
end
function Base.intersect!(L::HyperPlanesIntersection, h::HyperPlane)
if L.orthogonalized
h_p = remproj(h, L)
if isapproxzero(h_p.a)
if !isapproxzero(h_p.β)
# `L` is empty
empty!(L)
end
else
push!(L.hyperplanes, h_p)
end
else
push!(L.hyperplanes, h)
end
return L
end
@vecrepelem HyperPlanesIntersection HyperPlane hyperplanes
hreptype(::Type{<:HyperPlanesIntersection{T, AT}}) where {T, AT} = Intersection{T, AT}
# Returns an HyperPlanesIntersection representing the affine hull of p.
# The affine hull is defined as
# {λx + (1-λ)y | x, y ∈ p, λ ∈ R}
function affinehull(h::HRep, current=false)
if !current
detecthlinearity!(h)
end
HyperPlanesIntersection(FullDim(h), hyperplanes(h))
end
function detecthlinearity!(L::HyperPlanesIntersection, solver = nothing)
if !L.orthogonalized
V = detecthlinearity(L)
L.hyperplanes = V.hyperplanes
L.orthogonalized = true
end
return L
end
function detecthlinearity(L::HyperPlanesIntersection{T, AT}, solver=nothing) where {T, AT}
H = HyperPlanesIntersection{T, AT}(FullDim(L))
for h in hyperplanes(L)
hp = remproj(h, H)
if !(h in H)
intersect!(H, h)
end
end
H
end
# V-representation
struct VEmptySpace{T, AT, D <: FullDim} <: VLinearSpace{T}
d::D
end
FullDim(v::VEmptySpace) = v.d
function emptyspace(v::VRep{T}) where {T}
d = FullDim(v)
return VEmptySpace{T, vvectortype(typeof(v)), typeof(d)}(d)
end
Base.in(v::VRepElement, L::VEmptySpace) = isapproxzero(v)
"""
vrep(lines::LineIt; d::FullDim)
Creates an affine space of full dimension `d` from the list of lines `lines`.
### Examples
```julia
vrep([Line([1, 0, 0]), Line([0, 1, 0])])
```
creates the 2-dimensional affine subspace containing all the points ``(x_1, x_2, 0)``, i.e. the ``x_1````x_2``-plane.
"""
vrep(lines::LineIt; d::FullDim = FullDim_rec(lines)) = LinesHull(d, lines)
# Representation of an affine space containing the origin by the minkowsky sum of lines
mutable struct LinesHull{T, AT, D<:FullDim} <: VLinearSpace{T}
d::D
lines::Vector{Line{T, AT}}
orthogonalized::Bool
function LinesHull{T, AT, D}(d::FullDim, lines::LineIt{T}) where {T, AT, D}
new{T, AT, D}(FullDim_convert(D, d), lazy_collect(lines), false)
end
end
function LinesHull(d::FullDim, it::ElemIt{Line{T, AT}}) where {T, AT}
return LinesHull{T, AT, typeof(d)}(d, it)
end
LinesHull{T, AT}(d::FullDim) where {T, AT} = LinesHull(d, Line{T, AT}[])
FullDim(v::LinesHull) = v.d
vvectortype(::Type{<:LinesHull{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:LinesHull}, d::FullDim, ::Type{T}) where T = LinesHull{T, similar_type(vvectortype(PT), d, T), typeof(d)}
function convexhull!(L::LinesHull, l::Line)
if L.orthogonalized
l_p = remproj(l, L)
if !isapproxzero(l_p)
push!(L.lines, l_p)
end
else
push!(L.lines, l)
end
return L
end
@vecrepelem LinesHull Line lines
conetype(::Type{LinesHull{T, AT, D}}) where {T, AT, D} = RaysHull{T, AT, D}
vreptype(::Type{LinesHull{T, AT, D}}) where {T, AT, D} = Hull{T, AT, D}
# Returns a LinesHull representing the following set (TODO does it have a name?, does someone has a reference talking about it ?)
# {x | y in p => x + y in p}
function linespace(v::VRep, current=false)
if !current
detectvlinearity!(v)
end
LinesHull(FullDim(v), lines(v))
end
function detectvlinearity!(L::LinesHull, solver = nothing)
if !L.orthogonalized
V = detectvlinearity(L)
L.lines = V.lines
L.orthogonalized = true
end
return L
end
function detectvlinearity(L::LinesHull{T, AT}, solver = nothing) where {T, AT}
V = LinesHull{T, AT}(FullDim(L))
V.orthogonalized = true
for l in lines(L)
convexhull!(V, l)
end
return V
end
detectlinearity!(h::HRepresentation) = detecthlinearity!(h)
detectlinearity!(v::VRepresentation) = detectvlinearity!(v)
_linearity(h::HRepresentation) = hyperplanes(h)
_linearity(v::VRepresentation) = lines(v)
function remproj(el::RepElement, L::Union{HyperPlanesIntersection, LinesHull})
detectlinearity!(L)
# Modified Gram-Schmidt process, it is more stable than
# the classical Gram-Schmidt process.
# see https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process#Numerical_stability
for lin in _linearity(L)
el = remproj(el, lin)
end
return el
end
function Base.in(el::RepElement, L::Union{HyperPlanesIntersection, LinesHull})
return isapproxzero(remproj(el, L))
end