{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "struct CutoffPointIndex\n",
    "    cutoff::Int\n",
    "    index::Int\n",
    "end\n",
    "Base.show(io::IO, p::CutoffPointIndex) = print(io, \"p[$(p.cutoff), $(p.index)]\")\n",
    "struct CutoffRayIndex\n",
    "    cutoff::Int\n",
    "    index::Int\n",
    "end\n",
    "Base.show(io::IO, r::CutoffRayIndex) = print(io, \"r[$(r.cutoff), $(r.index)]\")\n",
    "struct DoubleDescriptionData{PointT, RayT, LineT, HST}\n",
    "    fulldim::Int\n",
    "    halfspaces::Vector{HST}\n",
    "    # Elements ordered by first halfspace cutting it off\n",
    "    points::Vector{PointT}\n",
    "    pz::Vector{BitSet}\n",
    "    cutpoints::Vector{Vector{PointT}}\n",
    "    cutpz::Vector{Vector{BitSet}}\n",
    "    pin::Vector{Vector{CutoffPointIndex}}\n",
    "    rays::Vector{RayT}\n",
    "    rz::Vector{BitSet}\n",
    "    cutrays::Vector{Vector{RayT}}\n",
    "    cutrz::Vector{Vector{BitSet}}\n",
    "    rin::Vector{Vector{CutoffRayIndex}}\n",
    "    lines::Vector{LineT}\n",
    "    cutline::Vector{Union{Nothing, LineT}}\n",
    "    lineray::Vector{Union{Nothing, CutoffRayIndex}}\n",
    "    nlines::Vector{Int}\n",
    "end\n",
    "function Base.show(io::IO, data::DoubleDescriptionData)\n",
    "    println(io, \"DoubleDescriptionData in $(data.fulldim) dimension:\")\n",
    "    println(io, data.points)\n",
    "    println(io, data.rays)\n",
    "    println(io, data.lines)\n",
    "    for i in reverse(eachindex(data.cutpoints))\n",
    "        println(io, \" Halfspace $i: $(data.halfspaces[i]):\")\n",
    "        if !isempty(data.cutpoints[i])\n",
    "            println(io, \"  Cut points:\")\n",
    "            for j in eachindex(data.cutpoints[i])\n",
    "                println(io, \"  $j: \", data.cutpoints[i][j], \" zero at: \", data.cutpz[i][j])\n",
    "            end\n",
    "        end\n",
    "        if !isempty(data.pin[i])\n",
    "            println(io, \"  In: \", data.pin[i])\n",
    "        end\n",
    "        if !isempty(data.cutrays[i])\n",
    "            println(io, \"  Cut rays:\")\n",
    "            for j in eachindex(data.cutrays[i])\n",
    "                println(io, \"  $j: \", data.cutrays[i][j], \" zero at: \", data.cutrz[i][j])\n",
    "            end\n",
    "        end\n",
    "        if !isempty(data.rin[i])\n",
    "            println(io, \"  In: \", data.rin[i])\n",
    "        end\n",
    "        if data.cutline[i] !== nothing\n",
    "            println(io, \"  Cut line: \", data.cutline[i])\n",
    "            if data.lineray[i] !== nothing\n",
    "                println(io, \"  Line ray: \", data.lineray[i])\n",
    "            end\n",
    "        end\n",
    "        if !iszero(data.nlines[i])\n",
    "            println(io, \"  $(data.nlines[i]) uncut lines left\")\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "function DoubleDescriptionData{PointT, RayT, LineT}(fulldim::Integer, halfspaces) where {PointT, RayT, LineT}\n",
    "    n = length(halfspaces)\n",
    "    return DoubleDescriptionData{PointT, RayT, LineT, eltype(halfspaces)}(\n",
    "        fulldim,\n",
    "        halfspaces,\n",
    "        PointT[],\n",
    "        BitSet[],\n",
    "        [PointT[] for i in 1:n],\n",
    "        [BitSet[] for i in 1:n],\n",
    "        [CutoffPointIndex[] for i in 1:n],\n",
    "        RayT[],\n",
    "        BitSet[],\n",
    "        [RayT[] for i in 1:n],\n",
    "        [BitSet[] for i in 1:n],\n",
    "        [CutoffRayIndex[] for i in 1:n],\n",
    "        LineT[],\n",
    "        Union{Nothing, LineT}[nothing for i in 1:n],\n",
    "        Union{Nothing, CutoffRayIndex}[nothing for i in 1:n],\n",
    "        zeros(Int, n)\n",
    "    )\n",
    "end\n",
    "function Base.getindex(data::DoubleDescriptionData, p::CutoffPointIndex)\n",
    "    if p.cutoff == 0\n",
    "        return data.points[p.index]\n",
    "    else\n",
    "        return data.cutpoints[p.cutoff][p.index]\n",
    "    end\n",
    "end\n",
    "function Base.getindex(data::DoubleDescriptionData, r::CutoffRayIndex)\n",
    "    if r.cutoff == 0\n",
    "        return data.rays[r.index]\n",
    "    else\n",
    "        return data.cutrays[r.cutoff][r.index]\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "isin (generic function with 2 methods)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function _bitdot_range(b1::BitSet, b2::BitSet, i, n)\n",
    "    count = 1 # They share the hyperplance `i`\n",
    "    for j in (i+1):n\n",
    "        if j in b1 && j in b2\n",
    "            count += 1\n",
    "        end\n",
    "    end\n",
    "    return count\n",
    "end\n",
    "function isadjacent(data, i::Integer, p1::CutoffPointIndex, p2::CutoffPointIndex)\n",
    "    pz1 = data.cutpz[p1.cutoff][p1.index]\n",
    "    pz2 = data.cutpz[p2.cutoff][p2.index]\n",
    "    n = length(data.halfspaces)\n",
    "    return _bitdot_range(pz1, pz2, i, n) + data.lines[i] + 1 == data.fulldim\n",
    "end\n",
    "function isadjacent(data, i::Integer, p::CutoffPointIndex, r::CutoffRayIndex)\n",
    "    pz = data.cutpz[p.cutoff][p.index]\n",
    "    rz = data.cutrz[r.cutoff][r.index]\n",
    "    n = length(data.halfspaces)\n",
    "    return _bitdot_range(pz, rz, i, n) + data.lines[i] + 1 == data.fulldim\n",
    "end\n",
    "function isadjacent(data, i::Integer, r::CutoffRayIndex, p::CutoffPointIndex)\n",
    "    return isadjacent(data, i, p, r)\n",
    "end\n",
    "function isadjacent(data, i::Integer, r1::CutoffPointIndex, r2::CutoffPointIndex)\n",
    "    rz1 = data.cutrz[r1.cutoff][r1.index]\n",
    "    rz2 = data.cutrz[r2.cutoff][r2.index]\n",
    "    n = length(data.halfspaces)\n",
    "    return _bitdot_range(rz1, rz2, i, n) + data.lines[i] + 2 == data.fulldim\n",
    "end\n",
    "function isin(data, i, p::CutoffPointIndex)\n",
    "    if p.cutoff == 0\n",
    "        return i in data.pz[p.index]\n",
    "    else\n",
    "        return i in data.cutpz[p.cutoff][p.index]\n",
    "    end\n",
    "end\n",
    "function isin(data, i, r::CutoffRayIndex)\n",
    "    if r.cutoff == 0\n",
    "        return i in data.rz[r.index]\n",
    "    else\n",
    "        return i in data.cutrz[r.cutoff][r.index]\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "add_index! (generic function with 4 methods)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resized_bitset(data) = sizehint!(BitSet(), length(data.halfspaces))\n",
    "function add_index!(data, cutoff::Nothing, p::AbstractVector)\n",
    "    push!(data.points, p)\n",
    "    push!(data.pz, resized_bitset(data))\n",
    "    return CutoffPointIndex(0, length(data.points))\n",
    "end\n",
    "function add_index!(data, cutoff::Integer, p::AbstractVector)\n",
    "    push!(data.cutpoints[cutoff], p)\n",
    "    push!(data.cutpz[cutoff], resized_bitset(data))\n",
    "    return CutoffPointIndex(cutoff, length(data.cutpoints[cutoff]))\n",
    "end\n",
    "function add_index!(data, cutoff::Nothing, r::Polyhedra.Ray)\n",
    "    push!(data.rays, r)\n",
    "    push!(data.rz, resized_bitset(data))\n",
    "    return CutoffRayIndex(0, length(data.rays))\n",
    "end\n",
    "function add_index!(data, cutoff::Integer, r::Polyhedra.Ray)\n",
    "    push!(data.cutrays[cutoff], r)\n",
    "    push!(data.cutrz[cutoff], resized_bitset(data))\n",
    "    return CutoffRayIndex(cutoff, length(data.cutrays[cutoff]))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "add_adjacent_element! (generic function with 1 method)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function add_in!(data, i, index::CutoffPointIndex)\n",
    "    push!(data.pin[i], index)\n",
    "    if index.cutoff == 0\n",
    "        push!(data.pz[index.index], i)\n",
    "    else\n",
    "        push!(data.cutpz[index.cutoff][index.index], i)\n",
    "    end\n",
    "end\n",
    "function add_in!(data, i, index::CutoffRayIndex)\n",
    "    push!(data.rin[i], index)\n",
    "    if index.cutoff == 0\n",
    "        push!(data.rz[index.index], i)\n",
    "    else\n",
    "        push!(data.cutrz[index.cutoff][index.index], i)\n",
    "    end\n",
    "end\n",
    "function set_in!(data, I, el, index)\n",
    "    for i in I\n",
    "        if el in Polyhedra.hyperplane(data.halfspaces[i])\n",
    "            add_in!(data, i, index)\n",
    "        end\n",
    "    end\n",
    "end\n",
    "function add_element!(data, k, el)\n",
    "    cutoff = nothing\n",
    "    for i in reverse(1:k)\n",
    "        if data.cutline[i] !== nothing\n",
    "            el = line_project(el, data.cutline[i], data.halfspaces[i])\n",
    "            index = add_adjacent_element!(data, i - 1, el, data.lineray[i])\n",
    "            set_in!(data, i:k, el, index)\n",
    "            return index\n",
    "        end\n",
    "        if !(el in data.halfspaces[i])\n",
    "            cutoff = i\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    index = add_index!(data, cutoff, el)\n",
    "    set_in!(data, (index.cutoff+1):k, el, index)\n",
    "    return index\n",
    "end\n",
    "function add_adjacent_element!(data, k, el, parent)\n",
    "    index = add_element!(data, k, el)\n",
    "    # Condition (c_k) in [FP96]\n",
    "    if index.cutoff != parent.cutoff\n",
    "        addintersection!(data, index, parent)\n",
    "    end\n",
    "    return index\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "combine (generic function with 4 methods)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "function combine(β, p1::AbstractVector, value1, p2::AbstractVector, value2)\n",
    "    λ = (value2 - β) / (value2 - value1)\n",
    "    return λ * p1 + (1 - λ) * p2\n",
    "end\n",
    "function combine(β, p::AbstractVector, pvalue, r::Polyhedra.Ray, rvalue)\n",
    "    λ = (β - pvalue) / rvalue\n",
    "    return p + λ * r\n",
    "end\n",
    "combine(β, r::Polyhedra.Ray, rvalue, p::AbstractVector, pvalue) = combine(β, p, pvalue, r, rvalue)\n",
    "function combine(r1::Polyhedra.Ray, value1, r2::Polyhedra.Ray, value2)\n",
    "    # should take\n",
    "    # λ = value2 / (value2 - value1)\n",
    "    @assert 0 <= value2 / (value2 - value1) <= 1\n",
    "    # By homogeneity we can avoid the division and do\n",
    "    #newr = value2 * r1 - value1 * r2\n",
    "    # but this can generate very large numbers (see JuliaPolyhedra/Polyhedra.jl#48)\n",
    "    # so we still divide\n",
    "    newr = (value2 * r1 - value1 * r2) / (r2 - r1)\n",
    "    # In CDD, it does value2 * r1 - value1 * r2 but then it normalize the ray\n",
    "    # by dividing it by its smallest nonzero entry (see dd_CreateNewRay)\n",
    "    return Polyhedra.simplify(newr)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "add_if_adjacent! (generic function with 1 method)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function addintersection!(data, idx1, idx2)\n",
    "    if idx1.cutoff > idx2.cutoff\n",
    "        return addintersection!(data, idx2, idx1)\n",
    "    end\n",
    "    i = idx2.cutoff\n",
    "    if isin(data, i, idx1)\n",
    "        return\n",
    "    end\n",
    "    el1 = data[idx1]\n",
    "    el2 = data[idx2]\n",
    "    h = data.halfspaces[i]\n",
    "    newel = combine(h.β, el1, h.a ⋅ el1, el2, h.a ⋅ el2)\n",
    "    add_adjacent_element!(data, i - 1, newel, idx1)\n",
    "end\n",
    "function add_if_adjacent!(data, i::Integer, el1, el2)\n",
    "    # Condition (c_k) in [FP96]\n",
    "    if el1.cutoff != el2.cutoff\n",
    "        if isadjacent(data, i, el1, el2)\n",
    "            addintersection!(data, el1, el2)\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "double_description (generic function with 1 method)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_shift(el::AbstractVector, line::Line) = el + Polyhedra.coord(line)\n",
    "_shift(el::Line, line::Line) = el + line\n",
    "_shift(el::Ray, line::Line) = el + Polyhedra.Ray(Polyhedra.coord(line))\n",
    "function line_project(el, line, h)\n",
    "    # (line + λ * cutline) ⋅ h.a == h.β\n",
    "    # λ = (h.β - line ⋅ h.a) / (cutline ⋅ h.a)\n",
    "    λ = (h.β - el ⋅ h.a) / (line ⋅ h.a)\n",
    "    return Polyhedra.simplify(_shift(el, λ * line))\n",
    "end\n",
    "function hline(data, line::Line, i, h)\n",
    "    value = h.a ⋅ line\n",
    "    if !Polyhedra.isapproxzero(value)\n",
    "        if data.cutline[i] === nothing\n",
    "            if value > 0\n",
    "                line = -line # Make `lineray` point inward\n",
    "            end\n",
    "            data.cutline[i] = line\n",
    "            cut = true\n",
    "            return true\n",
    "        else\n",
    "            line = line_project(line, data.cutline[i], hs)\n",
    "        end\n",
    "    end\n",
    "    data.nlines[i] += 1\n",
    "    return false\n",
    "end\n",
    "function double_description(hr::HRepresentation)\n",
    "    v = Polyhedra.dualfullspace(hr)\n",
    "    hps = Polyhedra.lazy_collect(hyperplanes(hr))\n",
    "    hss = Polyhedra.lazy_collect(halfspaces(hr))\n",
    "    data = DoubleDescriptionData{pointtype(v), raytype(v), linetype(v)}(fulldim(hr), hps, hss)\n",
    "    for line in lines(v)\n",
    "        cut = false\n",
    "        for i in reverse(eachindex(hps))\n",
    "            cut = hline(data, line, nhalfspaces(hr) + i, hss[i])\n",
    "            if cut\n",
    "                break\n",
    "            end\n",
    "        end\n",
    "        if !cut\n",
    "            for i in reverse(eachindex(hss))\n",
    "                cut = hline(data, line, i, hss[i])\n",
    "            end\n",
    "            if cut\n",
    "                break\n",
    "            end\n",
    "        end\n",
    "        if !cut\n",
    "            push!(data.lines, line)\n",
    "        end\n",
    "    end\n",
    "    # Add line rays after all lines are added so that the rays can be `line_project`ed.\n",
    "    # We only do that for halfspaces, hyperplanes do not create rays from cutoff lines.\n",
    "    # We use increasing index order since higher index may need the `lineray` of lower index.\n",
    "    for i in eachindex(hss)\n",
    "        line = data.cutline[i]\n",
    "        if line !== nothing\n",
    "            ray = Polyhedra.Ray(Polyhedra.coord(line))\n",
    "            data.lineray[i] = add_element!(data, i - 1, ray)\n",
    "        end\n",
    "    end\n",
    "    @assert isone(npoints(v))\n",
    "    add_element!(data, nhalfspaces(hr), first(points(v))) # Add the origin\n",
    "    for i in reverse(eachindex(hss))\n",
    "        if isempty(data.cutpoints[i]) && isempty(data.cutrays[i])\n",
    "            # Redundant, remove its contribution to avoid incorrect `isadjacent`\n",
    "            for p in data.pin\n",
    "                if p.cutoff == 0\n",
    "                    delete!(data.pz, i)\n",
    "                else\n",
    "                    delete!(data.cutpz[p.cutoff], i)\n",
    "                end\n",
    "            end\n",
    "            for r in data.rin\n",
    "                if r.cutoff == 0\n",
    "                    delete!(data.rz, i)\n",
    "                else\n",
    "                    delete!(data.cutrz[pr.cutoff], i)\n",
    "                end\n",
    "            end\n",
    "            continue\n",
    "        end\n",
    "        if i > 1\n",
    "            for p1 in data.pin[i], p2 in data.pin[i]\n",
    "                add_if_adjacent!(data, i, p1, p2)\n",
    "            end\n",
    "            for p in data.pin[i], r in data.rin[i]\n",
    "                add_if_adjacent!(data, i, p, r)\n",
    "            end\n",
    "        end\n",
    "        deleteat!(data.cutpoints, i)\n",
    "        deleteat!(data.cutpz, i)\n",
    "        if i > 1\n",
    "            for r1 in data.rin[i], r2 in data.rin[i]\n",
    "                add_if_adjacent!(data, i, r1, r2)\n",
    "            end\n",
    "        end\n",
    "        deleteat!(data.cutrays, i)\n",
    "        deleteat!(data.cutrz, i)\n",
    "        deleteat!(data.pin, i)\n",
    "        deleteat!(data.rin, i)\n",
    "    end\n",
    "    if isempty(data.points)\n",
    "        # Empty polyhedron, there may be rays left,\n",
    "        # Example 1: for 0x_1 + x_2 = -1 ∩ 0x_1 + x_2 = 1, the line (0, 1) is detected as correct\n",
    "        # Example 2: for 0x_1 + 0x_2 = 1, the lines (1, 0) and (0, 1) are detected as correct\n",
    "        # but since there is no point, the polyhedron is empty and we should drop all rays/lines\n",
    "        empty!(data.lines)\n",
    "        empty!(data.rays)\n",
    "    end\n",
    "    similar(v, data.points, data.lines, data.rays)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples\n",
    "\n",
    "## Square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}:\n",
       "4-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
       " HalfSpace([1.0, 0.0], 1.0)\n",
       " HalfSpace([-1.0, 0.0], 1.0)\n",
       " HalfSpace([0.0, 1.0], 1.0)\n",
       " HalfSpace([0.0, -1.0], 1.0)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Polyhedra\n",
    "h = HalfSpace([1.0, 0.0], 1.0) ∩ HalfSpace([-1.0, 0.0], 1.0) ∩\n",
    "    HalfSpace([0.0, 1.0], 1.0) ∩ HalfSpace([0.0, -1.0], 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000197 seconds (155 allocations: 9.531 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Float64,Array{Float64,1},Int64}:\n",
       "4-element iterator of Array{Float64,1}:\n",
       " [-1.0, -1.0]\n",
       " [1.0, -1.0]\n",
       " [-1.0, 1.0]\n",
       " [1.0, 1.0]"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time double_description(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000984 seconds (1.51 k allocations: 70.188 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Float64,Array{Float64,1},Int64}:\n",
       "4-element iterator of Array{Float64,1}:\n",
       " [1.0, 1.0]\n",
       " [-1.0, 1.0]\n",
       " [1.0, -1.0]\n",
       " [-1.0, -1.0]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time doubledescription(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Triangle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "H-representation Polyhedra.Intersection{Float64,Array{Float64,1},Int64}:\n",
       "3-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
       " HalfSpace([1.0, 1.0], 1.0)\n",
       " HalfSpace([1.0, -1.0], 0.0)\n",
       " HalfSpace([-1.0, 0.0], 0.0)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = hrep([HalfSpace([ 1.,  1], 1),\n",
    "          HalfSpace([ 1., -1], 0),\n",
    "          HalfSpace([-1.,  0], 0)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000179 seconds (141 allocations: 8.641 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Float64,Array{Float64,1},Int64}:\n",
       "3-element iterator of Array{Float64,1}:\n",
       " [0.0, 0.0]\n",
       " [0.0, 1.0]\n",
       " [0.5, 0.5]"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time double_description(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000830 seconds (1.26 k allocations: 61.578 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Float64,Array{Float64,1},Int64}:\n",
       "3-element iterator of Array{Float64,1}:\n",
       " [0.5, 0.5]\n",
       " [0.0, 0.0]\n",
       " [0.0, 1.0]"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time doubledescription(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "H-representation Polyhedra.Intersection{Float64,SArray{Tuple{2},Float64,1,2},Size{(2,)}}:\n",
       "3-element iterator of HalfSpace{Float64,SArray{Tuple{2},Float64,1,2}}:\n",
       " HalfSpace([1.0, 1.0], 1.0)\n",
       " HalfSpace([1.0, -1.0], 0.0)\n",
       " HalfSpace([-1.0, 0.0], 0.0)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using StaticArrays\n",
    "sh = HalfSpace((@SVector [ 1,  1]), 1) ∩\n",
    "     HalfSpace((@SVector [ 1, -1]), 0) ∩\n",
    "     HalfSpace((@SVector [-1,  0]), 0)\n",
    "shf = HalfSpace((@SVector [ 1,  1]), 1.0) ∩\n",
    "      HalfSpace((@SVector [ 1, -1]), 0.0) ∩\n",
    "      HalfSpace((@SVector [-1,  0]), 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000688 seconds (8.06 k allocations: 165.719 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Rational{BigInt},SArray{Tuple{2},Rational{BigInt},1,2},Size{(2,)}}:\n",
       "3-element iterator of SArray{Tuple{2},Rational{BigInt},1,2}:\n",
       " Rational{BigInt}[0//1, 0//1]\n",
       " Rational{BigInt}[0//1, 1//1]\n",
       " Rational{BigInt}[1//2, 1//2]"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time double_description(sh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.009481 seconds (43.65 k allocations: 923.055 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Rational{BigInt},SArray{Tuple{2},Rational{BigInt},1,2},Size{(2,)}}:\n",
       "3-element iterator of SArray{Tuple{2},Rational{BigInt},1,2}:\n",
       " Rational{BigInt}[1//2, 1//2]\n",
       " Rational{BigInt}[0//1, 0//1]\n",
       " Rational{BigInt}[0//1, 1//1]"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time doubledescription(sh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000447 seconds (250 allocations: 12.750 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Float64,SArray{Tuple{2},Float64,1,2},Size{(2,)}}:\n",
       "3-element iterator of SArray{Tuple{2},Float64,1,2}:\n",
       " [0.0, 0.0]\n",
       " [0.0, 1.0]\n",
       " [0.5, 0.5]"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time double_description(shf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.002401 seconds (3.48 k allocations: 142.328 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation Polyhedra.Hull{Float64,SArray{Tuple{2},Float64,1,2},Size{(2,)}}:\n",
       "3-element iterator of SArray{Tuple{2},Float64,1,2}:\n",
       " [0.5, 0.5]\n",
       " [0.0, 0.0]\n",
       " [0.0, 1.0]"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time doubledescription(shf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "H-representation MixedMatHRep{Float64,Array{Float64,2}}:\n",
       "3-element iterator of HalfSpace{Float64,Array{Float64,1}}:\n",
       " HalfSpace([1.0, 1.0], 1.0)\n",
       " HalfSpace([1.0, -1.0], 0.0)\n",
       " HalfSpace([-1.0, 0.0], 0.0)"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hmat = hrep([ 1  1\n",
    "              1 -1\n",
    "             -1  0],\n",
    "            [1., 0, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000121 seconds (130 allocations: 7.969 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation MixedMatVRep{Float64,Array{Float64,2}}:\n",
       "3-element iterator of Array{Float64,1}:\n",
       " [0.0, 0.0]\n",
       " [0.0, 1.0]\n",
       " [0.5, 0.5]"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time double_description(hmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.000910 seconds (1.67 k allocations: 84.406 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "V-representation MixedMatVRep{Float64,Array{Float64,2}}:\n",
       "3-element iterator of Array{Float64,1}:\n",
       " [0.5, 0.5]\n",
       " [0.0, 0.0]\n",
       " [0.0, 1.0]"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time doubledescription(hmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}