[WIP] affine

examples/Controlled Invariant Set.ipynb

@@ -39,7 +39,7 @@
     "0 & 0\n",
     "\\end{bmatrix}.\n",
     "$$\n",
-    "As shown in [LTJ18], a set $S$ is controlled invariant if\n",
+    "As shown in [LTJ18], a set $S$ is controlled invariant iff\n",
     "$$\n",
     "\\begin{bmatrix}\n",
     "1 & \\Delta t\n",
@@ -59,7 +59,9 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true
+   },
    "outputs": [],
    "source": [
     "using Polyhedra\n",
@@ -77,7 +79,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true
+   },
    "outputs": [],
    "source": [
     "using SCS\n",
@@ -88,7 +92,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true
+   },
    "outputs": [],
    "source": [
     "using CSDP    # CSDP is less appropriate than SCS and Mosek because it does not natively support\n",
@@ -99,7 +105,9 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true
+   },
    "outputs": [],
    "source": [
     "using MosekTools\n",
@@ -1942,7 +1950,7 @@
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.1.0"
+   "version": "1.1.1"
   }
  },
  "nbformat": 4,

src/Sets/ellipsoids.jl

@@ -92,7 +92,7 @@ function _HPH(D, d, δ)
          d D]
 end
 
-_HPH(ell::EllipsoidAtOrigin) = _HPH(ell.Q, zeros(size(q.Q, 1)), -1.0)
+_HPH(ell::EllipsoidAtOrigin) = _HPH(ell.Q, zeros(size(ell.Q, 1)), -1.0)
 
 """
     struct ShiftedEllipsoid{T}

src/Sets/transformations.jl

@@ -7,6 +7,16 @@ perspective_variable(li::LinearImage) = perspective_variable(li.set)
 space_variables(li::LinearImage) = nothing
 dimension(li::LinearImage) = size(li.A, 1)
 
+# A^{-1} * S + b
+struct AffineImage{S, T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractSet{T}
+    set::S
+    A::MT
+    b::VT
+end
+perspective_variable(ai::AffineImage) = perspective_variable(ai.set)
+space_variables(ai::AffineImage) = nothing
+dimension(ai::AffineImage) = size(ai.A, 1)
+
 # A^{-1} * S
 struct LinearPreImage{S, T, MT <: AbstractMatrix{T}} <: AbstractSet{T}
     set::S

src/affine.jl

@@ -8,6 +8,10 @@ struct AffineExpression{T, F <: AbstractScalarFunction} <: AbstractScalarFunctio
     constant::T
 end
 
+function Base.copy(aff::AffineExpression)
+    return AffineExpression(copy(aff.terms), aff.constant)
+end
+
 function Base.:(+)(f::F, g::F) where F <: AbstractScalarFunction
     return AffineExpression([AffineTerm(1.0, f), AffineTerm(1.0, g)], 0.0)
 end

src/inclusion.jl

@@ -50,6 +50,10 @@ function set_space(space::Space, ::InclusionConstraint{<:Sets.LinearImage,
                                                        <:Sets.LinearImage})
     return set_space(space, DualSpace)
 end
+function set_space(space::Space, ::InclusionConstraint{<:Sets.AffineImage,
+                                                       <:Sets.LinearImage})
+    return set_space(space, DualSpace)
+end
 # AS ⊆ S <=> S ⊆ A^{-1}S so PrimalSpace work
 # AS ⊆ S <=> A^{-T}S∘ ⊆ S∘ so DualSpace work
 function set_space(space::Space, ::InclusionConstraint{<:Sets.LinearImage,
@@ -140,14 +144,23 @@ function JuMP.build_constraint(_error::Function,
 end
 
 # See [LTJ18]
+function _apply_maps(_error, subset, supset; kws...)
+    dim = Sets.dimension(subset)
+    @polyvar x[1:dim]
+    JuMP.build_constraint(_error, apply_map(subset, x),
+                          PowerSet(apply_map(supset, x)); kws...)
+end
 function JuMP.build_constraint(_error::Function,
                                subset::Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}},
                                sup_powerset::PowerSet{<:Sets.LinearImage{<:Union{Sets.Polar, Sets.PerspectiveDual}}};
                                kws...)
-    dim = Sets.dimension(subset)
-    @polyvar x[1:dim]
-    JuMP.build_constraint(_error, apply_map(subset, x),
-                          PowerSet(apply_map(sup_powerset.set, x)); kws...)
+    _apply_maps(_error, subset, sup_powerset.set; kws...)
+end
+function JuMP.build_constraint(_error::Function,
+                               subset::Sets.AffineImage{<:Sets.PerspectiveDual},
+                               sup_powerset::PowerSet{<:Sets.LinearImage{<:Sets.PerspectiveDual}};
+                               kws...)
+    _apply_maps(_error, subset, sup_powerset.set; kws...)
 end
 function JuMP.build_constraint(_error::Function,
                                subset::Sets.AbstractSet,

src/map.jl

@@ -2,6 +2,11 @@ need_variablify(lm::Sets.LinearImage) = need_variablify(lm.set)
 function variablify(lm::Sets.LinearImage)
     return Sets.LinearImage(variablify(lm.set), lm.A)
 end
+need_variablify(lm::Sets.AffineImage) = need_variablify(lm.set)
+function variablify(lm::Sets.AffineImage)
+    return Sets.AffineImage(variablify(lm.set), lm.A, lm.b)
+end
+
 
 """
     apply_map(li::Sets.LinearImage{<:Sets.PolarOf{<:Sets.EllipsoidAtOrigin}}, new_vars)
@@ -64,7 +69,28 @@ function apply_map(li::Sets.LinearImage{<:Sets.PerspectiveDualOf{<:Sets.Househol
     return Sets.perspective_dual(dual)
 end
 
+# [x'] = [A b] [x]
+# [z'] = [0 1] [z]
+# The transpose is
+# [A' 0]
+# [b' 1]
+# So we substitute x for A' x' and then z for b' x' + z'
+function apply_map(li::Sets.AffineImage{<:Sets.PerspectiveDualOf{<:Sets.Householder{T}}},
+                   new_vars) where T
+    z = Sets.perspective_variable(li.set)
+    old_vars = Sets.space_variables(li.set)
+    q = subs(Sets.perspective_gauge0(li.set.set), old_vars => li.A' * new_vars)
+    # We do this substitution after in case `new_vars` and `old_vars` have an intersection
+    p = subs(q, z => z + dot(li.b, new_vars))
+    dual = Sets.Householder(Sets.UnknownSet{T}(), p, li.A * li.set.set.h + li.b,
+                            z, new_vars)
+    return Sets.perspective_dual(dual)
+end
+
 # FIXME, for Sets.AbstractSet, we should apply it directly
 function Base.:(*)(A::AbstractMatrix, set::Union{SetVariableRef, Sets.AbstractSet})
     return Sets.LinearImage(set, A)
 end
+function Base.:(+)(li::Sets.LinearImage, b::AbstractVector)
+    return Sets.AffineImage(li.set, li.A, b)
+end

src/variables.jl

@@ -1,3 +1,5 @@
+export InteriorPoint, CenterPoint
+
 abstract type AbstractVariable <: JuMP.AbstractVariable end
 
 abstract type HintPoint end