intersection of CPA and unbdd set

intersectionCPA.md

function intersection(X::CartesianProductArray{N},
                      Y::AbstractPolyhedron{N}) where {N}

    # preallocate the resulting CartesianProductArray with size hint
    result = CartesianProductArray(length(X), N)

    if isbounded(Y)
        # no free variables
        free_variables = []
    else
        free_variables = setdiff(1:dim(Y), constrained_dimensions(Y))
    end
    free_blocks = block_indices(X, free_variables)

    for bi in 1:length(X.array)
        if bi in free_blocks
            # if this block is not constrained, just push the set
            push!(result.array, X.array[i])
        else
            # otherwise, make the intersection with the projection of the halfspace
            push!(result.array, intersection(X.array[i], project(variable_indices(X, bi), Y))
        end
    end
    return result
end

I would sort the free_blocks and then keep track of the current index. If you have n blocks, you need to search n times with expected runtime of n/2, which gives you expected runtime of n²/2. On the other hand, sorting can be done in n log n and then you do not need to search.