Avoid zero division in decompose

docs/src/optimization.md

@@ -27,10 +27,10 @@ For instance, consider the 1-simplex:
 julia> using Polyhedra
 
 julia> simplex = HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0) ∩ HyperPlane([1, 1], 1)
-H-representation Polyhedra.Intersection{Int64,Array{Int64,1},Int64}:
-1-element iterator of HyperPlane{Int64,Array{Int64,1}}:
+H-representation Polyhedra.Intersection{Int64, Vector{Int64}, Int64}:
+1-element iterator of HyperPlane{Int64, Vector{Int64}}:
  HyperPlane([1, 1], 1),
-2-element iterator of HalfSpace{Int64,Array{Int64,1}}:
+2-element iterator of HalfSpace{Int64, Vector{Int64}}:
  HalfSpace([-1, 0], 0)
  HalfSpace([0, -1], 0)
 ```
@@ -47,26 +47,26 @@ CachingOptimizer state: NO_OPTIMIZER
 Solver name: No optimizer attached.
 
 julia> @variable(model, λ[1:2])
-2-element Array{VariableRef,1}:
+2-element Vector{VariableRef}:
  λ[1]
  λ[2]
 ```
 
 The variables can be constrained to belong to the simplex as follows:
 ```jldoctest jump-in-hrep
 julia> @constraint(model, λ in simplex)
-[λ[1], λ[2]] ∈ Polyhedra.PolyhedraOptSet{Int64,Polyhedra.Intersection{Int64,Array{Int64,1},Int64}}(HyperPlane([1, 1], 1) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0))
+[λ[1], λ[2]] ∈ Polyhedra.PolyhedraOptSet{Int64, Polyhedra.Intersection{Int64, Vector{Int64}, Int64}}(HyperPlane([1, 1], 1) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0))
 ```
 but a vector of affine or quadratic expressions can also be constrained to belong to the simplex:
 ```jldoctest jump-in-hrep
 julia> A = [1  1
             1 -1]
-2×2 Array{Int64,2}:
+2×2 Matrix{Int64}:
  1   1
  1  -1
 
 julia> @constraint(model, A * λ in simplex)
-[λ[1] + λ[2], λ[1] - λ[2]] ∈ Polyhedra.PolyhedraOptSet{Int64,Polyhedra.Intersection{Int64,Array{Int64,1},Int64}}(HyperPlane([1, 1], 1) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0))
+[λ[1] + λ[2], λ[1] - λ[2]] ∈ Polyhedra.PolyhedraOptSet{Int64, Polyhedra.Intersection{Int64, Vector{Int64}, Int64}}(HyperPlane([1, 1], 1) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0))
 ```
 We can verify that the model contains both constraints:
 ```julia

docs/src/plot.md

@@ -36,10 +36,10 @@ Suppose for instance that we want to visualize the polyhedron having the followi
 julia> using Polyhedra
 
 julia> v = convexhull([0, 0, 0]) + conichull([1, 0, 0], [0, 1, 0], [0, 0, 1])
-V-representation Polyhedra.Hull{Int64,Array{Int64,1},Int64}:
-1-element iterator of Array{Int64,1}:
+V-representation Polyhedra.Hull{Int64, Vector{Int64}, Int64}:
+1-element iterator of Vector{Int64}:
  [0, 0, 0],
-3-element iterator of Ray{Int64,Array{Int64,1}}:
+3-element iterator of Ray{Int64, Vector{Int64}}:
  Ray([1, 0, 0])
  Ray([0, 1, 0])
  Ray([0, 0, 1])
@@ -48,10 +48,10 @@ V-representation Polyhedra.Hull{Int64,Array{Int64,1},Int64}:
 The V-representation cannot be given to [MeshCat](https://github.com/rdeits/MeshCat.jl) or [Makie](https://github.com/JuliaPlots/Makie.jl) directly, it first need to be transformed into a polyhedron:
 ```jldoctest plots3
 julia> p = polyhedron(v)
-Polyhedron DefaultPolyhedron{Rational{BigInt},Polyhedra.Intersection{Rational{BigInt},Array{Rational{BigInt},1},Int64},Polyhedra.Hull{Rational{BigInt},Array{Rational{BigInt},1},Int64}}:
-1-element iterator of Array{Rational{BigInt},1}:
+Polyhedron DefaultPolyhedron{Rational{BigInt}, Polyhedra.Intersection{Rational{BigInt}, Vector{Rational{BigInt}}, Int64}, Polyhedra.Hull{Rational{BigInt}, Vector{Rational{BigInt}}, Int64}}:
+1-element iterator of Vector{Rational{BigInt}}:
  Rational{BigInt}[0//1, 0//1, 0//1],
-3-element iterator of Ray{Rational{BigInt},Array{Rational{BigInt},1}}:
+3-element iterator of Ray{Rational{BigInt}, Vector{Rational{BigInt}}}:
  Ray(Rational{BigInt}[1//1, 0//1, 0//1])
  Ray(Rational{BigInt}[0//1, 1//1, 0//1])
  Ray(Rational{BigInt}[0//1, 0//1, 1//1])
@@ -60,7 +60,7 @@ Polyhedron DefaultPolyhedron{Rational{BigInt},Polyhedra.Intersection{Rational{Bi
 Then, we need to create a mess from the polyhedron:
 ```jldoctest plots3
 julia> m = Polyhedra.Mesh(p)
-Polyhedra.Mesh{3,Rational{BigInt},DefaultPolyhedron{Rational{BigInt},Polyhedra.Intersection{Rational{BigInt},Array{Rational{BigInt},1},Int64},Polyhedra.Hull{Rational{BigInt},Array{Rational{BigInt},1},Int64}}}(convexhull([0//1, 0//1, 0//1]) + convexhull(Ray(Rational{BigInt}[1//1, 0//1, 0//1]), Ray(Rational{BigInt}[0//1, 1//1, 0//1]), Ray(Rational{BigInt}[0//1, 0//1, 1//1])), nothing, nothing, nothing)
+Polyhedra.Mesh{3, Rational{BigInt}, DefaultPolyhedron{Rational{BigInt}, Polyhedra.Intersection{Rational{BigInt}, Vector{Rational{BigInt}}, Int64}, Polyhedra.Hull{Rational{BigInt}, Vector{Rational{BigInt}}, Int64}}}(convexhull([0//1, 0//1, 0//1]) + convexhull(Ray(Rational{BigInt}[1//1, 0//1, 0//1]), Ray(Rational{BigInt}[0//1, 1//1, 0//1]), Ray(Rational{BigInt}[0//1, 0//1, 1//1])), nothing, nothing, nothing)
 ```
 
 ```@docs

src/decompose.jl

@@ -57,10 +57,19 @@ function scene(vr::VRep, ::Type{T}) where T
     # Intersection of rays with the limits of the scene
     (start, r) -> begin
         ray = coord(r)
-        times = max.((Vector(minimum(scene))-start) ./ ray, (Vector(maximum(scene))-start) ./ ray)
-        times[ray .== 0] .= Inf # To avoid -Inf with .../(-0)
-        time = minimum(times)
-        start + time * ray
+        λ = nothing
+        min_scene = minimum(scene)
+        max_scene = maximum(scene)
+        for i in 1:3
+            r = ray[i]
+            if !iszero(r)
+                cur = max((min_scene[i] - start[i]) / r, (max_scene[i] - start[i]) / r)
+                if λ === nothing || cur < λ
+                    λ = cur
+                end
+            end
+        end
+        start + λ * ray
     end
 end
 

src/vecrep.jl

@@ -23,10 +23,10 @@ For instance, the simplex
 can be created as follows:
 ```jldoctest
 julia> hrep([HyperPlane([1, 1], 1)], [HalfSpace([0, -1], 0), HalfSpace([-1, 0], 0)])
-H-representation Polyhedra.Intersection{Int64,Array{Int64,1},Int64}:
-1-element iterator of HyperPlane{Int64,Array{Int64,1}}:
+H-representation Polyhedra.Intersection{Int64, Vector{Int64}, Int64}:
+1-element iterator of HyperPlane{Int64, Vector{Int64}}:
  HyperPlane([1, 1], 1),
-2-element iterator of HalfSpace{Int64,Array{Int64,1}}:
+2-element iterator of HalfSpace{Int64, Vector{Int64}}:
  HalfSpace([0, -1], 0)
  HalfSpace([-1, 0], 0)
 ```