As mentioned more briefly in Tests Not Giving Good Error Info, I’m encountering a strange phantom error which only occurs when running tests or stepping thru via the debugger – when I just run the code normally, no error arises. This strikes me as something that should basically never happen, so I’m quite baffled.
I haven’t been able to construct toy examples which replicate the issue, so a decent amount of my code will be copied below to enable reproducing it. (This is almost the same code that I eventually copied in the other thread as well.)
The functions at the very end are the relevant ones. no_tests(50) runs fine, while with_tests(50) has an error… despite the only difference between them being the @test macros.
Although with my actual code I was also seeing the phantom error when stepping through via the debugger, with the code below I’m not able to replicate that for some reason – now I’m seeing the error when using @test, but stepping into no_tests(50) and then telling the debugger to continue appears to hang (whereas the function runs in a reasonable amount of time when not debugging).
As I mentioned in the other thread, I’m in the process of porting this code from an older version of Julia, so if something looks weird that’s my excuse. Any advice about the code itself is appreciated, but mainly I’m worried about this really weird phantom-error behavior.
Thanks for any insights!
using Test
using Random
srand(x) = Random.seed!(x)
using JuMP
using Clp
import Base.==
abstract type Split end
struct LinSplit <: Split
coefs::Vector{Float64}
con::Float64
end
function ==(l1::LinSplit,l2::LinSplit)
if l1.con!=l2.con
return false
else
for i in 1:max(length(l1.coefs), length(l2.coefs))
if get(l1.coefs, i, 0.0) != get(l2.coefs, i, 0.0)
return false
end
end
return true
end
end
LinSplit(coef,con) =
LinSplit(convert(Vector{Float64}, coef), convert(Float64, con))
abstract type BSPNode end
struct None <: BSPNode
end
struct SplitNode <: BSPNode
split::Split # split definition
neg::BSPNode # negative side
pos::BSPNode # positive side
end
==(s1::SplitNode, s2::SplitNode) =
s1.split == s2.split && s1.neg == s2.neg && s1.pos == s2.pos
struct LeafNode <: BSPNode
value::Any
end
==(l1::LeafNode,l2::LeafNode) =
l1.value == l2.value
abstract type Region end
struct LinBound
split::LinSplit
positive::Bool
end
==(l1::LinBound,l2::LinBound) =
l1.split == l2.split && l1.positive == l2.positive
import Base.isless
function isless(l1::LinBound, l2::LinBound) # defining so they can be sorted
if l1.positive < l2.positive
return true
elseif l2.positive < l1.positive
return false
elseif l1.split.con < l2.split.con
return true
elseif l2.split.con < l1.split.con
return false
else
for i in 1:length(l1.split.coefs)
if l1.split.coefs[i] < l2.split.coefs[i]
return true
elseif l2.split.coefs[i] < l1.split.coefs[i]
return false
end
end
return false
end
end
struct LinRegion <: Region
bounds::Vector{LinBound}
dimensionality::Int
end
function regularize_bound(b::LinBound)
z = sqrt(sum(b.split.coefs .^ 2))
if z == 0.0
return b
else
return LinBound(LinSplit(b.split.coefs / z, b.split.con / z), b.positive)
end
end
function ==(l1::LinRegion,l2::LinRegion)
b1 = map(regularize_bound, trim_region(l1).bounds)
b2 = map(regularize_bound, trim_region(l2).bounds)
return sort(b1) == sort(b2)
end
boundless_reg(d::Int) = LinRegion(Vector{LinBound}(),d)
function rect_split(dim,loc,numdims)
z = zeros(numdims)
z[dim]=1.0
return LinSplit(z,loc)
end
function rect_region(mins, maxs)
max_d = length(mins)
s = Vector{LinBound}()
for d in 1:max_d
z = zeros(max_d)
z[d] = 1.0
s = union(s,[LinBound(LinSplit(z,mins[d]), true)])
end
for d in 1:max_d
z = zeros(max_d)
z[d] = 1.0
s = union(s,[LinBound(LinSplit(z,maxs[d]), false)])
end
return LinRegion(s,max_d)
end
import Base.isequal
isequal(t1::SplitNode, t2::SplitNode) =
isequal(t1.split, t2.split) &&
isequal(t1.pos, t2.pos) &&
isequal(t1.neg, t2.neg)
isequal(t1::LinSplit, t2::LinSplit) =
t1.coefs == t2.coefs &&
t1.con == t2.con
struct BSPTree
boundary::Region
root::BSPNode
end
==(t1::BSPTree, t2::BSPTree) = (t1.boundary == t2.boundary) && isequal(t1.root, t2.root)
isequal(t1::BSPTree, t2::BSPTree) = t1 == t2
function randomrect(n)
mins = rand(n)
maxs = rand(n)+mins
return rect_region(mins,maxs)
end
function randomrect()
randomrect(2)
end
function randomsplit(n)
LinSplit(rand(n),rand())
end
function random_rect_split(n)
coefs = zeros(n)
coefs[rand(1:n)] = 1.0
LinSplit(coefs,rand())
end
# creates a balanced tree of the given depth
function randomtree(dims, depth)
if depth == 0
return LeafNode(rand())
else
return SplitNode(randomsplit(dims),
randomtree(dims,depth-1),
randomtree(dims,depth-1))
end
end
function random_rect_tree(dims, depth)
if depth == 0
return LeafNode(rand())
else
return SplitNode(random_rect_split(dims),
random_rect_tree(dims,depth-1),
random_rect_tree(dims,depth-1))
end
end
function epsilonequal(t1::BSPTree, t2::BSPTree, epsilon::Float64=0.0000001)
if !(regvalid(intersection(t1.boundary, t2.boundary), epsilon))
return true # If they have no intersection, they trivially agree about it!
elseif typeof(t1.root) == SplitNode
return epsilonequal(takepos(t1), t2, epsilon) && epsilonequal(takeneg(t1), t2, epsilon)
elseif typeof(t2.root) == SplitNode
return epsilonequal(takepos(t2), t1, epsilon) && epsilonequal(takeneg(t2), t1, epsilon)
elseif t1.root == None() && t2.root == None()
return true
elseif t1.root == None() || t2.root == None()
return false
else
return epsilonequal(t1.root.value, t2.root.value, epsilon)
end
end
epsilonequal(a::Float64, b::Float64, epsilon::Float64=0.0000001) = (a == b) || (abs(a - b) < epsilon)
epsilonequal(a::Number, b::Number, epsilon::Float64=0.0000001) = epsilonequal(convert(Float64,a), convert(Float64,b), epsilon)
epsilonequal(a::LinBound, b::LinBound, epsilon::Float64=0.0000001) =
a.positive == b.positive && epsilonequal(a.split, b.split, epsilon)
epsilonequal(a::LinSplit, b::LinSplit, epsilon::Float64=0.0000001) =
epsilonequal(a.con, b.con, epsilon) && all(map(((x,y) -> epsilonequal(x,y, epsilon)), a.coefs, b.coefs)) # shouldn't really be assuming the vectors have same length here
function intersection(r1::LinRegion, r2::LinRegion)
return LinRegion(union(r1.bounds, r2.bounds), max(r1.dimensionality, r2.dimensionality))
end
function rect_valid(r::Region)
max_d = r.dimensionality
maxs = fill(Inf, max_d)
maxs_closed = trues(max_d) # whether the maximum is closed
mins = fill(-Inf, max_d)
mins_closed = trues(max_d) # whether the minimum is closed
for bound in r.bounds # Go through bounds checking number of nonzero coefs; if >1, return 0; if 1, record the new min/max implied.
bound_coefs = bound.split.coefs
bound_length = length(bound_coefs)
if bound_length > max_d
throw("rect_valid handed a region with more coefficients than dimensions: $r")
end
nonzero_count = 0
nonzero_loc = 0
for dim in 1:bound_length
if !(bound_coefs[dim] == 0.0)
nonzero_count = nonzero_count + 1
if nonzero_count > 1
return 0
end
nonzero_loc = dim
end
end
if nonzero_loc > 0 #if the bound is doing anything at all
# Now that we know the bound is uni-dimensional, add it to the list of dimension bounds.
adjusted_con = bound.split.con/bound_coefs[nonzero_loc] # The actual location for this bound.
if bound.positive # closed bound
if bound_coefs[nonzero_loc] > 0.0 # positive coef, so bound.positive=true means the bound is expressing a minimum
if adjusted_con > mins[nonzero_loc] # lower bound would be increased
mins[nonzero_loc] = adjusted_con
mins_closed[nonzero_loc] = true # inherits the closed status, since we increased past whatever the bound previously might have been
end # if the bound is less than or equal to one already established, we have no need to update the open or closed status in this case, since it could only be loosening the constraint, which we don't want to do
elseif bound_coefs[nonzero_loc] < 0.0 # negative coef, so bound is actually a maximum, since -x>=c means x<=c
if adjusted_con < maxs[nonzero_loc] # existing upper bound would be decreased
maxs[nonzero_loc] = adjusted_con
maxs_closed[nonzero_loc] = true # inherits the open status, since we decreased past whatever the previous bound would have been
end # no need to update open/closed status otherwise
end
else # open bound
if bound_coefs[nonzero_loc] > 0.0 # positive coef, so bound is saying x<c, IE giving a max
if adjusted_con < maxs[nonzero_loc] # new max would be lower than old max
maxs[nonzero_loc] = adjusted_con # lower the max
maxs_closed[nonzero_loc] = false # inherit the open status
elseif adjusted_con == maxs[nonzero_loc] # if they're equal, we don't need to update the max, but we may need to switch it from closed to open since an open constraint is stricter
maxs_closed[nonzero_loc] = false # let's just assign; no point checking
end # otherwise, nothing to update; we were already maximally strict
elseif bound_coefs[nonzero_loc] < 0.0 && adjusted_con >= mins[nonzero_loc] # negative coef, so we're saying -x<c, IE x>c, so we are giving a min
if adjusted_con > mins[nonzero_loc] # new min would be higher
mins[nonzero_loc] = adjusted_con
mins_closed[nonzero_loc] = false # inherit the open status
elseif adjusted_con == mins[nonzero_loc] # if they're equal, we still need to enforce openness
mins_closed[nonzero_loc] = false
end
end
end
else # no nonzero locs in this bound, because we already checked whether there's 1, and returned for the case where there's more than one
if bound.positive && bound.split.con > 0 # no way to get > 0 with no coefficients; this region is invalid
return -1
elseif !bound.positive && bound.split.con <= 0 # if the constant is 0 or below, there's no way we can get strictly under it
return -1
end
end
end
# Now that we've added all the bounds, check for validity.
for dim in 1:max_d
if maxs[dim] < mins[dim] || (maxs[dim]==mins[dim] && !(maxs_closed[dim] && mins_closed[dim]))
return -1
end
end
return +1
end
function zero_bound_check(b::LinBound) # using the same idea as the zero-coef bound case in rect_valid, returns +1 if the bound rules out nothing, -1 if everything, otherwise zero
if all((x->x==0.0), b.split.coefs)
if b.positive
if b.split.con > 0
return -1
else
return +1
end
else
if b.split.con <= 0
return -1
else
return +1
end
end
else
return 0
end
end
function regvalid(r::LinRegion, epsilon::Float64=0.0000001)
r_v = rect_valid(r)
if r_v==1
return true
elseif r_v==-1
return false
end
m = makemodel(r)
for bound in r.bounds
if !(bound.positive)
d = -bound.split.coefs
@objective(m[1], Max, sum(d[i]*(m[2])[i] for i=1:length(d))) # optimize away from each face
optimize!(m[1])
s = termination_status(m[1])
if s == :Infeasible # not even the region's closure is valid
return false
elseif s == :Optimal
optimize_away = getobjectivevalue(m[1])
difference = optimize_away + bound.split.con # Since things are negative, I want the optimized value minus the **negative of** the constant; so really I want to add them.
if difference < epsilon
return false
end
end
end
end
return true
end
regvalid(r::BSPTree) =
regvalid(r.boundary)
takepos(reg::LinRegion, split::LinSplit) = LinRegion(union(reg.bounds, [LinBound(split, true)]), reg.dimensionality)
takeneg(reg::LinRegion, split::LinSplit) = LinRegion(union(reg.bounds, [LinBound(split, false)]), reg.dimensionality)
takepos(tree::BSPTree) =
BSPTree(takepos(tree.boundary, tree.root.split), tree.root.pos)
takeneg(tree::BSPTree) =
BSPTree(takeneg(tree.boundary, tree.root.split), tree.root.neg)
function makemodel(r::LinRegion)
m = Model(optimizer_with_attributes(Clp.Optimizer, "InfeasibleReturn" => 1, "LogLevel" => 0))
d = r.dimensionality
@variable(m,x[1:d])
for bound in r.bounds
if bound.positive
@constraint(m, sum(x[i]*bound.split.coefs[i] for i=1:length(bound.split.coefs)) >= bound.split.con)
else
@constraint(m, sum(x[i]*bound.split.coefs[i] for i=1:length(bound.split.coefs)) <= bound.split.con)
end
end
return (m, x)
end
function cheap_treearithmetic(operator, base_type)
eval(
quote
function $operator(n1::BSPNode, n2::$base_type) # Only work on raw nodes if we're combining with a number; otherwise we want to track regions to avoid computing things which would be trimmed, since that can make an exponential difference.
if typeof(n1) == LeafNode
return LeafNode( $operator(n1.value, n2))
elseif n1 == None()
# return None()
return LeafNode(n2) # None() combined with anything else acts like an identity element.
elseif typeof(n1) == SplitNode
return SplitNode(n1.split, $operator(n1.neg, n2),
$operator(n1.pos, n2))
end
end
function $operator(n1::$base_type, n2::BSPNode)
if typeof(n2) == LeafNode
return LeafNode( $operator(n1, n2.value))
elseif n2 == None()
#return None()
return LeafNode(n1) # None() combined with anything else acts like an identity element.
elseif typeof(n2) == SplitNode
return SplitNode(n2.split, $operator(n1, n2.neg),
$operator(n1, n2.pos))
end
end
function $operator(t1::BSPTree, t2::BSPTree)
i = intersection(t1.boundary, t2.boundary)
if typeof(t1.root) == LeafNode
if !regvalid(i)
return BSPTree(i, None())
end
return merge_bottom_up(BSPTree(i, ($operator(t1.root.value, cheap_trim(t2, i).root))))
elseif typeof(t1.root) == SplitNode
pos = takepos(t1)
neg = takeneg(t1)
return BSPTree(i,
SplitNode(t1.root.split,
identity($operator(neg, t2)).root, # I have to use "identity" here to prevent the order of operations from being changed when $operator is expanded.
identity($operator(pos, t2)).root)) #
elseif t1.root == None()
return t2
elseif t2.root == None()
return t1
end
end
$operator(t1::BSPTree, t2::$base_type) = BSPTree(t1.boundary, $operator(t1.root, t2))
$operator(t1::$base_type, t2::BSPTree) = BSPTree(t2.boundary, $operator(t1, t2.root))
end
)
end
cheap_treearithmetic(op) = cheap_treearithmetic(op, Number)
import Base.+
cheap_treearithmetic(:+)
import Base.*
cheap_treearithmetic(:*)
import Base.max
cheap_treearithmetic(:max)
import Base.min
cheap_treearithmetic(:min)
import Base./
cheap_treearithmetic(:/)
import Base.-
cheap_treearithmetic(:-)
import Base.^
cheap_treearithmetic(:^,Integer) # adding the integer definitions to eliminate some compiler warnings about ambiguous matches
cheap_treearithmetic(:^)
function with_tests(n)
for i=1:n
srand(i)
dims = rand(1:5)
a = BSPTree(randomrect(dims), randomtree(dims, rand(1:4)))
@test epsilonequal(a,a)
@test epsilonequal(a, a + 0.000001, 0.01)
@test !epsilonequal(a, a + 0.1, 0.01)
end
end
function no_tests(n)
for i=1:n
srand(i)
dims = rand(1:5)
a = BSPTree(randomrect(dims), randomtree(dims, rand(1:4)))
epsilonequal(a,a)
epsilonequal(a, a + 0.000001, 0.01)
!epsilonequal(a, a + 0.1, 0.01)
end
end
no_tests(50)
with_tests(50)