Julia: add modphase keyword for isapprox, deprecate `@test_approx_eq_modphase`

[StefanKarpinski]

Since it seems that both Base and some packages use the horrifically named @test_approx_eq_modphase macro, perhaps it would be useful for a modphase keyword to be introduced for the isapprox function at which point @test_approx_eq_modphase a b could be cleanly deprecated in favor of @test a ≈ b modphase=true instead of simply deleting the macro or copy-pasta-ing everywhere. I, however, do not feel that I have the expertise to correctly implement such a keyword argument for the function, so if someone who understands what this macro is for could implement it, that would be greatly appreciated. cc: @andreasnoack, @jiahao, @stevengj

[stevengj]

It should probably be just isapprox(a, b * cis(angle(b ⋅ a))), since this is the angle ϕ that minimizes the L2 norm |a - b eⁱᵠ|.

[StefanKarpinski]

So would it make sense to implement this as something like

function isapprox(...; modphase::Bool=false)
    modphase && (b *= cis(angle(b⋅a)))
    # rest of isapprox definition
end

Or do you think it simply shouldn't be part of isapprox at all?

[stevengj]

Something like that seems fine to me. I would tend to modphase && return isapprox(a,b*fixphase(b⋅a); kws...) rather than *= to ensure type stability of b, and I would probably define:

fixphase(bdota::Real) = sign(bdota)
fixphase(bdota) = cis(angle(bdota))

in order to avoid unnecessarily complexifying real vectors.

[stevengj]

I assume you mean the Base.Test.test_approx_eq_modphase function, since there is no such macro. One complication here is that it tests equality columnwise (e.g. for an array of eigenvectors, each column of which may have a different phase), so it is a bit different than isapprox.

[stevengj]

In particular, I think the scalar and vector cases should look like this:

# given bdota = dot(b,a), return the phase factor eⁱᵠ that minimizes
# the L2 norm |b - a|.   For real vectors, this is just a ± sign,
# and we special-case bdota::Real to avoid complexifying real vectors.
fixphase(bdota::Real) = bdota < 0 ? -one(bdota) : +one(bdota)
fixphase(bdota) = isfinite(bdota) ? cis(angle(bdota)) : complex(float(fixphase(real(bdota))))

function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x,y), atol::Real=0,
                                        nans::Bool=false, ignorephase::Bool=false)
    ignorephase && return isapprox(a,b*fixphase(b⋅a); rtol=rtol, atol=atol, nans=nans)
    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y))) || (nans && isnan(x) && isnan(y))
end

function isapprox{T,S}(x::AbstractArray{T}, y::AbstractArray{S};
                  rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y)),
                  atol::Real=0, nans::Bool=false, ignorephase::Bool=false,
                  norm::Function=vecnorm)
    if ignorephase
        ydotx = y⋅x
        if isnan(ydotx) && nans
            # Re-compute dot products ignoring nans, since NaN*phase is a NaN.
            # (We can use reduce(op, itr) because x & y must be nonempty here.)
            ydotx = reduce(ab -> let adotb = ab[1]⋅ab[2]
                                     isnan(adotb) ? zero(adotb) : adotb
                                 end, zip(y, x))
        end
        return isapprox(x,y*Base.fixphase(ydotx); rtol=rtol, atol=atol, nans=nans, norm=norm)
    end
    d = norm(x - y)
    if isfinite(d)
        return d <= atol + rtol*max(norm(x), norm(y))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
    end
end

function test_approx_eq_modphase{S,T}(a::AbstractVector{S}, b::AbstractVector{T},
                                      err = length(indices(a,1))^3*(eps(S)+eps(T)))
    @test a ≈ b atol=err ignorephase=true
end

function test_approx_eq_modphase{S,T}(A::AbstractMatrix{S}, B::AbstractMatrix{T},
                                      err = length(indices(a,1))^3*(eps(S)+eps(T)))
    @test all(i -> isapprox(A[:,i], B[:,i], atol=err, ignorephase=true), i in size(A,2))
end

However, the all(...) construction is verbose enough that I think we would end up still wanting something like the test_approx_eq_modphase for our tests, so I'm not sure what the deprecation path is.

[StefanKarpinski]

Couldn't the all construction be the default behavior? I'm a bit unclear why that's necessary.

[stevengj]

It could, but then the modphase keyword is rather odd, because its behavior would depend on the dimensionality of the array. What would it do for 3d arrays?

[StefanKarpinski]

I have no idea, I don't really understand what this function does or why it exists, which is why I want to deleted it. I feel like you (@stevengj), @andreasnoack and @jiahao need to figure this one out.

[stevengj]

The reason it exists is for testing things like eigensolvers, which return a matrix whose columns are the eigenvectors, but each eigenvector has a "random" phase.

[StefanKarpinski]

Would a function to "standardize" the phase of certain slices of an array allow this to be expressed more concisely?

[andreasnoack]

I think this is too specific to linear algebra to be part of isapprox. We could just have a small convenience function in LinAlg for this instead of making isapprox more complex. E.g. something like

_make_first_row_positive(A) = A*Diagonal(conj.(sign.(A[1,:])))
@test _make_first_row_positive(V1) ≈ _make_first_row_positive(V2)

or simply

@test V1 ≈ V2 .* cis.(angle.(sum(conj.(V2) .* V1, 1)))

in the few cases where it is needed.

[fredrikekre]

I removed this in secret in #25571