`unique()` with `isapprox()` instead of `isequal()`

[wsshin]

I hope to have the version of unique() that uses isapprox() instead of isequal(). Such a version will be useful, for example, to pick unique floating point numbers. MATLAB has uniquetol() for that purpose.

Instead of creating a separate function, we may generalize unique() such that it accepts an "equality" function as an extra argument. Then we could provide isapprox() or other equality functions as the extra argument.

[stevengj]

This would be a nice feature to have. Unfortunately, it requires a completely different implementation than the ones we have now (for arrays and generic iterables), since those are based on hashing. Changing the equality function requires changing the hash function, but I'm not sure how you would make a hash function compatible with a fuzzy comparison like isapprox.

Developing an efficient unique algorithm that works with an arbitrary equivalence relation seems like an interesting problem. I wonder if there is a good known algorithm for this?

[stevengj]

For equivalence relations, there is the disjoint-set data structure. But the problem is more difficult here because isapprox is not even an equivalence relation because it is not transitive: 1 ≈ 1+1e-8 and 1+1e-8 ≈ 1+2e-8, but !(1 ≈ 1+2e-8).

The lack of transitivity makes it unclear to me how uniquetol is even defined. Can you give a precise definition?

[KristofferC]

I dont want to copy matlabs documentation. https://se.mathworks.com/help/matlab/ref/uniquetol.html says explicitly how its defined.

[stevengj]

@KristofferC that documentation does not precisely define it, because it doesn't define how the non-transitivity is handled.

For example, what is uniquetol([1, 1+1e-8, 1+2e-8, 1+3e-8], 1.5e-8)? Valid answers would apparently include [1], [1, 1+2e-8], [1,1+3e-8], and [1+1e-8,1+3e-8]. Its documentation does not seem to distinguish between the validity of these three options. Of course, Matlab must give some answer, but "whatever Matlab does" is not a good definition.

[stevengj]

For an ordered set, like an array, I guess you could define it like:

function uniquetol{T}(A::AbstractArray{T}; kws...)
     S = T[]
     for a in A
          if !any(s -> isapprox(s, a; kws...), S)
              push!(S, a)
          end
     end
     return S
end

which is an O(N^2) algorithm. But the answer depends on the ordering of A, and so you wouldn't want to use this definition for an unordered collection like a Set (unless you don't mind having implementation-defined results).

For a real-valued array, you could devise an O(N log N) algorithm by sorting, but doing so for something like a complex array seems harder (one could use a kd tree, I guess, or some other data structure for nearest-neighbor search).

[standarddeviant]

I'm late to the party, but a naive implementation is

UPDATE: fixed the function to actually work...

function uniquetol_map(itr,tol=1e-6)
    # utVals = Array{Number}(0)
    ut_vals_map = Dict()
    ut_idxs_map = Dict()
    if length(itr) == 0
        return utMap
    end

    _itr = copy(itr)
    while(length(_itr) > 0)
        idxs = find(abs.(itr - _itr[1]) .< tol)
        utval = mean(itr[idxs])
        ut_vals_map[ utval ] = itr[idxs]
        ut_idxs_map[ utval ] =     idxs
        deleteat!(_itr,find(abs.(_itr - _itr[1]) .< tol))
    end
    
    return ut_vals_map, ut_idxs_map
end

I agree that there's probably a pathological test-case that will change behavior of almost any uniquetol implementation, but working 99.9% of the time is better than 0%. I'm using this to find repeated roots of a polynomial.

[stevengj]

@standarddeviant, that implementation is O(N^2), allocates lots of unnecessary temporary arrays, and still has the problem of producing results that depend in a confusing way on input order as I mentioned above. It is surely fine for some purposes, but I'd be reluctant to include something like this in Base.

[StefanKarpinski]

I'm also against having this in base. Fundamentally, isapprox does not define an equivalence relation, so any behavior here is massively implementation-dependent. The unique function works by hashing values in a way that respects isequal so unique with an arbitrary comparison function would have to work completely differently. What would make sense is unique with a transformation function which returns a representative for each unique value of the transformation (i.e. transformation implicitly induces equivalence classes on its input via preimage).

I agree that there's probably a pathological test-case that will change behavior of almost any uniquetol implementation, but working 99.9% of the time is better than 0%. I'm using this to find repeated roots of a polynomial.

I've written this before, but "this mostly works" is a luxury that user code and lightly used packages can afford, but which a language's standard library cannot. Anything that goes in the standard library has to work in all cases, not just the common ones. If you want to have a uniquetol function, you could make a PR to the MatlabCompat package.

[standarddeviant]

I wanted to update this thread just in terms of the algorithm complexity to do this, not to reopen the issue.

I made a version that is much more performant than the above naive approach here :

julia> a = randn(10) + im*randn(10);

julia> @time uniquetol_list(a);
  0.000102 seconds (211 allocations: 12.563 KiB)

julia> @time uniquetol_dict(a);
  0.033212 seconds (227 allocations: 98.828 KiB)

Re: performance, I'm not a CS person, but the heavy lifting is all done by sort and diff. I believe the default sort in julia (merge sort) has a worst case of n log(n). The real parts get sorted once, and the imaginary parts get sorted in smaller chunks - this seems good because you're incurring n log(n) for multiple smaller values of n. I tried to re-use the buffers I already had, although I didn't see a diff! function.

Also, because of the sorting, this approach yields the same answer independent of input order.

[jamblejoe]

As @StefanKarpinski mentioned, it's better to transform your collection and call unique on the transformed collection. For comparing floating point numbers, one could specify a rounding scheme like so:

A = [1.123, 1.124, 1.0]
unique(A) == A
unique(x->round(x, digits=2), A) == [1.123, 1.0]