Get all changes of sign of some vectors

General Usage Performance
1

Hello,

Given some vectors of numbers, I need to get all vectors obtained from the given ones by assigning the signs + or - to the numbers. For example:

vctrs = [[0.4; 0.3], [0.5; 0.0]]
julia> changesOfSign(vctrs)
6-element Vector{Vector{Float64}}:
 [-0.4, -0.3]
 [0.4, -0.3]
 [-0.4, 0.3]
 [0.4, 0.3]
 [-0.5, 0.0]
 [0.5, 0.0]

Here is my implementation:

function expandgrid(vectors) # Cartesian product
    cp = Iterators.product(ntuple(i -> vectors[i], length(vectors))...)
    return hcat(map(collect, collect(cp))...)
end

function changesOfSign(vectors)
    vvv = [[ifelse(x==0, [0.0], [-x, x]) for x=v] for v=vectors]
    mats = [expandgrid(vv) for vv=vvv]
    mat = hcat(mats...)
    return [mat[:, i] for i=1:size(mat, 2)]
end

Is it good? I’d like to receive your comments or suggestions, if you have ones.

2

Hi @stla,

from the first look at it there should be quite some room for improvement. But why don’t you post a complete M(inimal)W(orking)E(xample) so we know the size of the data set you want to apply your implementation too?

Bonus tip: best would be something involving @btime from BenchmarkTools.jl:wink:

3

Oh, that depends… For example I use it with 5 vectors of dimension 4. That’s rather a small dataset.

4

This is a version that is faster, but not very pretty either:

julia> vctrs = [[0.4; 0.3], [0.5; 0.0]]
2-element Vector{Vector{Float64}}:
 [0.4, 0.3]
 [0.5, 0.0]

julia> function signs(vecs,::Val{N}) where N
           list = eltype(vecs)[]
           signs = Iterators.product(ntuple( i -> -1:2:1, N)...)
           for sign in signs
               for v in vecs
                   push!(list, [ v[i] == 0 ? 0. : sign[i] * v[i] for i in 1:length(v) ])
               end
           end
           return unique(list)
       end
signs (generic function with 3 methods)

julia> @btime signs($vctrs,Val(length($vctrs[1])))
  934.520 ns (16 allocations: 1.25 KiB)
6-element Vector{Vector{Float64}}:
 [-0.4, -0.3]
 [-0.5, 0.0]
 [0.4, -0.3]
 [0.5, 0.0]
 [-0.4, 0.3]
 [0.4, 0.3]

It is sort of ugly for two reasons: first, because to avoid a lot of allocations, you need to specialize the signs iterator at compile time by providing the length of the vectors. Second, because of the [ v[i] == 0 ? 0. : sign[i] * v[i] for i in 1:length(v) ], which is needed here only because unique does not eliminate the cases where -0.0 results not to be considered identical to 0.0:

julia> unique([ [-0.0], [0.0] ])
2-element Vector{Vector{Float64}}:
 [-0.0]
 [0.0]

julia> unique([ [0.0], [0.0] ])
1-element Vector{Vector{Float64}}:
 [0.0]

If it wasn’t for this, that line would be just push!(list,sign .* v) and the function would be prettier.

1 Like
5

Here is a version that is faster than both, and does not require the Val argument. But, it was a lot of effort and is actually not much faster on bigger inputs.

Removing the push! and preallocating with similar takes off about 50ns, using unique! another 50ns, the new algorithm takes off a further 100ns or so. Making this branchless did not help.

using BenchmarkTools
vctrs = [[0.4; 0.3], [0.5; 0.0]]
bigger_vctrs = [[0.4; 0.3; 0.1; 0.0], [0.5; 0.0; -0.2; 0.3], [0.4; 0.3; 0.1; 0.0], [0.0; 0.2; 0.3; -0.4]]


function signs(vecs, ::Val{N}) where {N}
    list = eltype(vecs)[]
    signs = Iterators.product(ntuple(i -> -1:2:1, N)...)
    for sign in signs
        for v in vecs
            push!(list, [v[i] == 0 ? 0.0 : sign[i] * v[i] for i in 1:length(v)])
        end
    end
    return unique(list)
end

@info "baseline"
@btime signs($vctrs, Val(length($vctrs[1])))
@btime signs($bigger_vctrs, Val(length($bigger_vctrs[1])))



function signs_new(vecs)
    M = length(vecs)
    N = length(vecs[1])
    nsigns = 2^N
    list = [similar(vecs[1]) for i in 1:nsigns*M]
    for (i, v) in enumerate(vecs)
        for j in 1:nsigns
            ctr = copy(j - 1)
            ind = j + (i - 1) * nsigns
            for k in 1:N
                if v[k] != 0.0 && ctr & 0b1 == 1
                    list[ind][k] = v[k]
                else
                    list[ind][k] = -v[k]
                end
                ctr = ctr >> 1
            end
        end
    end
    unique!(list)
return list

@assert all(sort(signs_new(vctrs)) .== sort(signs(vctrs, Val(2))))
@assert all(sort(signs_new(bigger_vctrs)) .== sort(signs(bigger_vctrs, Val(4))))

@btime signs_new($vctrs)
@btime signs_new($bigger_vctrs)

julia> include("bench.jl")
[ Info: baseline
  632.316 ns (16 allocations: 1.25 KiB)
  3.923 μs (78 allocations: 9.47 KiB)
[ Info: new version
  419.445 ns (13 allocations: 1.11 KiB)
  3.354 μs (72 allocations: 7.67 KiB)

if OP is able to use StaticArray, then that is probably where the real speedups lie:

julia> @btime signs_new(bigger_vctrs)
  3.145 μs (72 allocations: 7.67 KiB)

julia> @btime signs_new($(SVector{4}.(bigger_vctrs)))
  1.853 μs (72 allocations: 4.67 KiB)
1 Like
6

I erred. Based on @lmiq and @pjentsch0 proposals I have

using Random, BenchmarkTools

function expandgrid(vectors) # Cartesian product
    cp = Iterators.product(ntuple(i -> vectors[i], length(vectors))...)
    return hcat(map(collect, collect(cp))...)
end

function changesOfSign(vectors)
    vvv = [[ifelse(x==0, [0.0], [-x, x]) for x=v] for v=vectors]
    mats = [expandgrid(vv) for vv=vvv]
    mat = hcat(mats...)
    return [mat[:, i] for i=1:size(mat, 2)]
end

function signs1(vecs, ::Val{N}) where {N}
    list = eltype(vecs)[]
    signs = Iterators.product(ntuple(i -> -1:2:1, N)...)
    for v in vecs
        for sign in signs
            push!(list, [v[i] == 0 ? 0.0 : sign[i] * v[i] for i in 1:length(v)])
        end
    end
    return unique(list)
end

function signs2(vecs)
    M = length(vecs)
    N = length(vecs[1])
    nsigns = 2^N
    list = [similar(vecs[1]) for i in 1:nsigns*M]
    for (i, v) in enumerate(vecs)
        for j in 1:nsigns
            ctr = copy(j - 1)
            ind = j + (i - 1) * nsigns
            for k in 1:N
                if v[k] != 0.0 && ctr & 0b1 == 1
                    list[ind][k] = v[k]
                else
                    list[ind][k] = -v[k]
                end
                ctr = ctr >> 1
            end
        end
    end
    unique!(list)
    return list
end

function signs3(vecs) 
    N = length(vecs[1])
    signs = Iterators.product(ntuple( i -> (-1, 1), N)...)
    return [[sign[i] * vec[i] for i in 1:N] 
        for vec in vecs for sign in signs 
            if all(vec[i] != 0 || sign[i] == 1 for i in eachindex(vec))]
end 


vctrs = [[0.4; 0.3; 0.2], [0.5; 0.0; -0.0]]
result1 = changesOfSign(vctrs)
@show result1
result2 = signs1(vctrs, Val(length(vctrs[1]))) 
@show result2
@assert result1 == result2
result2 = signs2(vctrs) 
@show result2
@assert result2 == result1
result2 = signs3(vctrs) 
@show result2
@assert result2 == result1

@btime changesOfSign(vctrs) setup=(Random.seed!(42); vctrs = [rand(10) for i in 1:10])
@btime signs1(vctrs, Val(length(vctrs[1]))) setup=(Random.seed!(42); vctrs = [rand(10) for i in 1:10])
@btime signs2(vctrs) setup=(Random.seed!(42); vctrs = [rand(10) for i in 1:10])
@btime signs3(vctrs) setup=(Random.seed!(42); vctrs = [rand(10) for i in 1:10])
println()

yielding

result1 = [[-0.4, -0.3, -0.2], [0.4, -0.3, -0.2], [-0.4, 0.3, -0.2], [0.4, 0.3, -0.2], [-0.4, -0.3, 0.2], [0.4, -0.3, 0.2], [-0.4, 0.3, 0.2], [0.4, 0.3, 0.2], [-0.5, 0.0, 0.0], [0.5, 
0.0, 0.0]]
result2 = [[-0.4, -0.3, -0.2], [0.4, -0.3, -0.2], [-0.4, 0.3, -0.2], [0.4, 0.3, -0.2], [-0.4, -0.3, 0.2], [0.4, -0.3, 0.2], [-0.4, 0.3, 0.2], [0.4, 0.3, 0.2], [-0.5, 0.0, 0.0], [0.5, 
0.0, 0.0]]
result2 = [[-0.4, -0.3, -0.2], [0.4, -0.3, -0.2], [-0.4, 0.3, -0.2], [0.4, 0.3, -0.2], [-0.4, -0.3, 0.2], [0.4, -0.3, 0.2], [-0.4, 0.3, 0.2], [0.4, 0.3, 0.2], [-0.5, -0.0, 0.0], [0.5, -0.0, 0.0]]
result2 = [[-0.4, -0.3, -0.2], [0.4, -0.3, -0.2], [-0.4, 0.3, -0.2], [0.4, 0.3, -0.2], [-0.4, -0.3, 0.2], [0.4, -0.3, 0.2], [-0.4, 0.3, 0.2], [0.4, 0.3, 0.2], [-0.5, 0.0, -0.0], [0.5, 0.0, -0.0]]
  1.146 ms (20864 allocations: 5.65 MiB)
  3.055 ms (10280 allocations: 2.23 MiB)
  3.159 ms (10264 allocations: 1.67 MiB)
  803.500 μs (10265 allocations: 1.73 MiB)
1 Like
7

Thank you guys. Sorry for the delay in replying, my Ubuntu laptop is broken.

8 
julia> v = [[0.4; 0.3; 0.2], [0.5; 0.0; -0.0]]
2-element Vector{Vector{Float64}}:
 [0.4, 0.3, 0.2]
 [0.5, 0.0, -0.0]

julia> @btime [[i,j,k].*v[e] for i in [-1,1],j in [-1,1], k in [-1,1], e in [1,2]][:]
  5.200 μs (74 allocations: 4.39 KiB)
16-element Vector{Vector{Float64}}:
 [-0.4, -0.3, -0.2]
 [0.4, -0.3, -0.2]
 [-0.4, 0.3, -0.2]
 [0.4, 0.3, -0.2]
 [-0.4, -0.3, 0.2]
 [0.4, -0.3, 0.2]
 [-0.4, 0.3, 0.2]
 [0.4, 0.3, 0.2]
 [-0.5, -0.0, 0.0]
 [0.5, -0.0, 0.0]
 [-0.5, 0.0, 0.0]
 [0.5, 0.0, 0.0]
 [-0.5, -0.0, -0.0]
 [0.5, -0.0, -0.0]
 [-0.5, 0.0, -0.0]
 [0.5, 0.0, -0.0]