Count occurances for matrix rows (where column order does not matter)

General Usage
,
21

One fix is to track deleted items to prevent them from being re-added.

function onlyoneof2(F)
    Ft = sort.(SVector{4}.(eachrow(F)))
    d = Dict{eltype(Ft),Int}()
    del = Vector{eltype(Ft)}()
    for i in eachindex(Ft)
        if haskey(d, Ft[i])
            delete!(d, Ft[i])
            push!(del, Ft[i])
        elseif Ft[i] ∈ del
            continue
        else
            d[Ft[i]] = i
        end
    end
    return collect(values(d))
end

Now if we modify the short F you gave to have an odd number of occurrences of a row: Fβ€²=[1 2 3 4; 5 6 7 8; 4 3 2 1; 7 8 6 5; 1 2 5 6; 8 7 6 5; 5 6 7 8;1 2 3 4;].

julia> onlyoneof(F)
2-element Vector{Int64}:
 5
 8 # wrong

julia> onlyoneof2(F)
1-element Vector{Int64}:
 5

I’m surprised (and somewhat confused) to find that this new version is faster!

julia> @benchmark onlyoneof(F4) evals=5
BenchmarkTools.Trial: 6531 samples with 5 evaluations.
 Range (min … max):  136.330 ΞΌs … 470.626 ΞΌs  β”Š GC (min … max): 0.00% … 59.64%
 Time  (median):     141.526 ΞΌs               β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   152.676 ΞΌs Β±  50.194 ΞΌs  β”Š GC (mean Β± Οƒ):  6.92% Β± 12.69%

  β–†β–ˆβ–‚β–‚                                                       ▁▂ ▁
  β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–„β–…β–ƒβ–„β–β–„β–ƒβ–„β–ƒβ–β–ƒβ–β–ƒβ–β–β–ƒβ–ƒβ–β–ƒβ–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–ˆβ–ˆ β–ˆ
  136 ΞΌs        Histogram: log(frequency) by time        384 ΞΌs <

 Memory estimate: 405.27 KiB, allocs estimate: 210.

julia> @benchmark onlyoneof2(F4) evals=5
BenchmarkTools.Trial: 9865 samples with 5 evaluations.
 Range (min … max):   89.992 ΞΌs … 308.097 ΞΌs  β”Š GC (min … max): 0.00% … 64.94%
 Time  (median):      95.344 ΞΌs               β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   100.918 ΞΌs Β±  31.735 ΞΌs  β”Š GC (mean Β± Οƒ):  5.52% Β± 11.16%

  β–†β–ˆβ–‡β–                                                    ▁▁  ▁ β–‚
  β–ˆβ–ˆβ–ˆβ–ˆβ–†β–†β–ˆβ–…β–„β–„β–„β–„β–ƒβ–β–„β–…β–ƒβ–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–„β–ˆβ–ˆβ–‡β–†β–ˆ β–ˆ
  90 ΞΌs         Histogram: log(frequency) by time        280 ΞΌs <

 Memory estimate: 286.78 KiB, allocs estimate: 15.
22

this doesn’t scale well to large F. Try benchmarking with 10_000 rows (at least)

23

That does scale poorly compared to your function with increasing duplicates. Without a better picture of real F data, I don’t know whether that would be a problem in practice.

24

Here is the MATLAB code for it, which on my machine can process a (25143270,4) array in about 6 seconds:

[~,ind1,ind2]=unique(sort(F,2),'rows'); %use unique
c=accumarray(ind2,1,[length(ind1) 1]); %Counts for unique set
L=c(ind2)==1;  %Expand counts to match input, and check for occurrence = 1

@Dan’s onlyoneofX code, and also my version onlyoneof3 below (based on @halleysfifthinc’s code above :point_up: but using a Dict to track the deleted entries) both take about 7.5-8 seconds (I just used @elapsed a couple times) on my machine for a (25143270,4) array.

function onlyoneof3(F)
    Ft = sort.(SVector{4}.(eachrow(F)))
    d = Dict{eltype(Ft),Int}()
    del = Dict{eltype(Ft),Int}()
    for i in eachindex(Ft)
        if haskey(d, Ft[i]) & ~haskey(del, Ft[i])
            delete!(d, Ft[i]) #Remove from d             
            del[Ft[i]] = i #Add to deleted
        else
            d[Ft[i]] = i
        end
    end
    return collect(values(d))
end

Here is the code I am using to create that β€œgiant ish” face array:

function grid3D(x,y,z)

    X = [ i for i ∈ x          , j ∈ 1:length(y), k ∈ 1:length(z) ]
    Y = [ j for i ∈ 1:length(x), j ∈ y          , k ∈ 1:length(z) ]
    Z = [ k for i ∈ 1:length(x), j ∈ 1:length(y), k ∈ z           ]

    return X, Y, Z
end


function getElements(ijk)
    
    #Cartesian index offsets
    iStep = CartesianIndex(1,0,0)
    jStep = CartesianIndex(0,1,0)
    kStep = CartesianIndex(0,0,1)

    F=zeros(Int64,6*length(ijk),4)
    for q=1:1:length(ijk)
        #Build 8-noded hex element
        e=[LinearIndices(sizV)[ijk[q]] LinearIndices(sizV)[ijk[q]+iStep] LinearIndices(sizV)[ijk[q]+iStep+jStep] LinearIndices(sizV)[ijk[q]+jStep]  LinearIndices(sizV)[ijk[q]+kStep] LinearIndices(sizV)[ijk[q]+iStep+kStep] LinearIndices(sizV)[ijk[q]+iStep+jStep+kStep] LinearIndices(sizV)[ijk[q]+jStep+kStep]]                

        F[1+(q-1)*6:q*6,:]=[e[[4 3 2 1]]; 
                            e[[5 6 7 8]]; 
                            e[[1 2 6 5]];
                            e[[2 3 7 6]];
                            e[[3 4 8 7]];
                            e[[4 1 5 8]];
                            ]
    end

    return F
end

s=0.01
r=1;
X,Y,Z= grid3D(-r:s:r,-r:s:r,-r:s:r)

M = sqrt.(X.^2 .+ Y.^2 .+ Z.^2)

siz=size(M)
sizV=siz.+1

L=M.<=(r+s/100) #Bool defining segmentation
ijk=findall(L) #Cartesian indices
F = getElements(ijk)
25

Yes it matters :slight_smile: I think I might have fixed your code, but performance is similar to that shorter one by Dan. I now also provided an example code to create big F arrays too.

1 Like
26

Is it still too slow for your needs? b/c it can be made faster (>2x) than MATLAB. But I feel, correct me if wrong, this is decent performance.

27

It would be good to find that β€œconvincingly faster than MATLAB” solution, as I am trying to convince myself and my colleagues that Julia is the way to go for the future. I also have ever larger meshes at times so any time saved would be great.

Do you think that using additional outputs from the unique function (like for the MATLAB version), would perhaps lead to the fastest approach?

28

The dictionary in Julia is essentially simlar to unique+accumarray from MATLAB and this is why the solutions take similar times. Julia is not magically faster in such small non-custom tasks (just as matrix-multiplication is the same speed in MATLAB and Julia).
But, Julia is fast overall, and open and easily extendible, so I would recommend it.
Any specific task, if scrutinized carefuly admits more performance squeeze. Do you really believe this problem will clinch it for the team?
(BTW same performance squeeze available in other β€œuncrippled” languages such as C/C++, javascript, java, rust).

30

countmap seems like a good idea. Combined with SVectors, the following is a little faster than the fastest onlyoneof above.

using StatsBase

function onlyoneof4(F)
    S = [sort(SVector{4}(r)) for r in eachrow(F)]
    countmap(S)
end
2 Likes
31

in order get fair comparison the sorted array of indices β€œmust” be returned :slight_smile:

1 Like
32

i found that a for large F, the sort based approach the β€œmost effective”,
here is the r=1.0, s=0.01 case:

julia-1.8> include("facestest.jl")
  Activating project at `~/Asztal/git/discourse.julialang.org/faces`
size of F: 25143270
β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”
β”‚                            fun β”‚  time(s) β”‚  mem(MB) β”‚ allocs(#) β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚                           sort β”‚ 3.78e+00 β”‚ 2.42e+03 β”‚        10 β”‚
β”‚               countmap by @dan β”‚ 4.37e+00 β”‚ 2.64e+03 β”‚        83 β”‚
β”‚ 2 dicts by @halleysfifthinc+OP β”‚ 5.25e+00 β”‚ 3.27e+03 β”‚       432 β”‚
β”‚                         1 dict β”‚ 3.86e+00 β”‚ 2.62e+03 β”‚        94 β”‚
β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜

perhaps tweaking in this direction would lead to more acceptable results.