[Bug]: Result value is Supposition.Fail even though the test function returns true

[sztal]

What happened?

Hi, first of all, thanks for the great work on Supposition.jl! And let me also start with apologies, since this will not be a proper bug report - I first want to ask on a more general level, since providing a reproducible example for my case may be quite difficult.

So, what happens is that in some tests run using @check macro I get Supposition.Fail result even though the result produced by the tested function is correct, and the function returns true when executed on the failing example.

I want to ask about any known nuances that may lead to such an outcome before providing a reproducible example, as this problem occurs for me when testing a somewhat complex code including a custom type inheriting from LazyArray on examples generated from a custom Possibility producing matrices.

The crux is that I am running the testing function on the examples that fail according to Supposition and I see that they are correct and the testing function returns true, so the situation is really confusing. More concrtely:

# Let this be 'Supposition.SuppositionReport' where `matrices` is a possibility generating values of type `Matrix`
report = @check testfunction(X=matrices) = begin
    # something is computed here and the output is `Bool`
end

# And the problem is that there are tests for which:
report.result.value isa Supposition.Fail
# while
testfunction(report.result.value.example...) === true

Is there any piece of knowledge about Supposition that could be helpful for solving this problem?

What did you expect to happen?

I expect examples for which the testing function returns true to produce Supposition.Pass results. More concretely:

testfunction(report.result.value.example...) === true
# should imply
report.result.value isa Supposition.Pass

Julia Version

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 9 5900HX with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_PROJECT = @.
  JULIA_PKG_PRESERVE_TIERED_INSTALLED = true
  JULIA_REVISE = manual

Package Environment

[7d9fca2a] Arpack v0.5.4
⌃ [86223c79] Graphs v1.11.2
  [5078a376] LazyArrays v2.2.1
  [5a0628fe] Supposition v0.3.5
  [37e2e46d] LinearAlgebra v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [f489334b] StyledStrings v1.11.0
  [8dfed614] Test v1.11.0
Info Packages marked with ⌃ have new versions available and may be upgradable.

[Seelengrab]

Hi! Thanks for the kind words :) If you're unsure whether your problem is a bug or not, you can also open a discussion topic. Apologies for only getting to it today!

Onto your issue - the only situation I can think of that this could happen in is if you reuse some object in your custom Possibility for generating your matrices, while also modifying that same object in your property 🤔 Supposition.jl (currently) assumes that the inputs to a property are isolated from each other. If that's not the case, it can happen that the shrinking process modifies a failing example into a "working" one, thereby hiding the failure.

Without a reproducible example (I don't mind if it's somewhat large, as long as it's reproducible!), your options are checking for the above and/or storing the first failure in a global Ref with a flag that it's already been filled. This should work, since you mention that the result is a Supposition.Fail and not a Supposition.Error, so you should be able to store the first failure just before returning from your property.

If that also leads to a failing result in spite of being true, I'll definitely need a reproducer!

[sztal]

Thanks! I will look into this and report here if that is okay.

Btw. could you give a simple example of a situation in which an object is incorrectly modified in a property that could lead to such an unexpected behavior?

I am still new to both Julia and Supposition, so I want to make sure I understand what you mean exactly.

[Seelengrab]

Inside of a property you should be able to modify the argument however much you want :) The input is recreated from scratch for new executions. What I assumed was that you defined a custom type subtyping Possibility, with a custom produce! method.

If that's not the case and you've only combined existing Possibility, it should work properly. Are you perhaps sharing some global variable in your custom Possibilitiy? An example showcasing the problem would be very helpful here!

[sztal]

Yes, I did define some custom Possibility subtypes. Below are the exact definitions, each code (sorry for the code dump!). Let me know if anything is unclear or stupid.

Here are possiblity types for generating real and complex scalars.

# numbers.jl
module _Numbers
export RealNumbers, ComplexNumbers

import Supposition: example, TestCase, reject
import Supposition.Data: Possibility, Floats, produce!, postype
import StyledStrings: @styled_str
import Senet.Testing: typename


struct RealNumbers{
    T <: AbstractFloat
} <: Possibility{T}
    floats::Floats{T}
    rtol::T

    function RealNumbers{T}(minimum=-T(Inf), maximum=+T(Inf); rtol=nothing) where {T}
        rtol = isnothing(rtol) ? sqrt(eps(T)) : rtol
        rtol > eps(T) ||
            throw(ArgumentError("'rtol' must be greater 'eps($T)'"))
        minimum < maximum ||
            throw(ArgumentError("'minimum' must be lower than 'maximum'"))
        floats = Floats{T}(; minimum=minimum, maximum=maximum, infs=false, nans=false)
        return new{T}(floats, rtol)
    end
end

RealNumbers(args...; kwargs...) = RealNumbers{Float64}(args...; kwargs...)

function Base.propertynames(r::RealNumbers, private::Bool=false)
    return (fieldnames(RealNumbers)..., :minimum, :maximum)
end

function Base.getproperty(r::RealNumbers, name::Symbol)
    if name == :minimum
        return r.floats.minimum
    elseif name == :maximum
        return r.floats.maximum
    else
        return getfield(r, name)
    end
end

Base.:(==)(r1::RealNumbers, r2::RealNumbers) = r1.floats == r2.floats

function Base.show(io::IO, r::RealNumbers)
    print(io, typename(r), "(")
    print(io, "; ")
    print(io, "minimum=", string(r.minimum), ", ")
    print(io, "maximum=", string(r.maximum), ")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", r::RealNumbers)
    left = isfinite(r.minimum) ? "[$(r.minimum)" : "(-∞"
    right = isfinite(r.maximum) ? "$(r.maximum)]" : "∞)"
    print(io, styled"""
    {code,underline:$RealNumbers}:

        Produce real number in {code:$(left), $(right)}.
    """)
end

function produce!(tc::TestCase, r::RealNumbers{T}) where {T}
    x = produce!(tc, r.floats)
    !issubnormal(x) || reject(tc)
    x = abs(x) < eps(T) ? zero(T) : x
    if !isnothing(r.rtol)
        abs(x) >= r.rtol || reject(tc)
    end
    return x
end


struct ComplexNumbers{T <: AbstractFloat} <: Possibility{Complex{T}}
    magnitudes::RealNumbers{T}

    function ComplexNumbers{T}(magnitudes::RealNumbers{T}) where {T}
        magnitudes.minimum >= 0 ||
            throw(ArgumentError("'magnitudes' must have non-negative minimum"))
        return new{T}(magnitudes)
    end
end

function ComplexNumbers{T}(; minabs=zero(T), maxabs=T(Inf), kwargs...) where {T<:Real}
    magnitudes = RealNumbers{T}(minabs, maxabs; kwargs...)
    return ComplexNumbers{T}(magnitudes)
end

function ComplexNumbers(; minabs=0.0, maxabs=Inf, kwargs...)
    minabs = convert(AbstractFloat, minabs)
    maxabs = convert(AbstractFloat, maxabs)
    T = promote_type(typeof(minabs), typeof(maxabs))
    return ComplexNumbers{T}(; minabs=minabs, maxabs=maxabs, kwargs...)
end

function Base.:(==)(cn1::ComplexNumbers, cn2::ComplexNumbers)
    return cn1.magnitudes == cn2.magnitudes
end

function Base.show(io::IO, cn::ComplexNumbers)
    print(io, typename(cn), "(")
    show(io, cn.magnitudes)
    print(io, ")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", cn::ComplexNumbers)
    min = cn.magnitudes.minimum
    max = cn.magnitudes.maximum
    print(io, styled"""
    {code,underline:$ComplexNumbers}:

        Produce complex number with magnitude from {code:$min} to {code:$max}.
    """)
end

function produce!(tc::TestCase, cn::ComplexNumbers{T}) where {T}
    z = randn(tc.rng, Complex{T})
    z /= abs(z)
    m = produce!(tc, cn.magnitudes)
    return z * m
end

end

And here for generating square matrices from matrices of (eigen)basis vectors and vectors of eigenvalues.

# matrices.jl
module _Matrices
export BasisVectors, SquareMatrices

import Supposition: example, TestCase
import Supposition.Data: Possibility, Integers, produce!, postype, reject
import StyledStrings: @styled_str

using LinearAlgebra:
    norm, rank, qr, diag, pinv, eigvecs, eigvals,
    Diagonal, Symmetric, Hermitian
using Senet.Testing: typename
using Senet.Testing._Numbers: RealNumbers, ComplexNumbers


struct BasisVectors{
    T <: Union{AbstractFloat, Complex},
} <: Possibility{Matrix{T}}
    mindim::Int
    maxdim::Int
    sizes::Possibility{Int}
    unitary::Bool
    deficiency::Int

    function BasisVectors{T}(mindim, maxdim; unitary=false, deficiency=0) where {T}
        mindim >= 0 ||
            throw(ArgumentError("'mindim' cannot be negative"))
        maxdim >= mindim ||
            throw(ArgumentError("'maxdim' cannot be lower than 'mindim'"))
        deficiency > mindim &&
            throw(ArgumentError("rank 'deficiency' cannot be greater than 'mindim'"))
        unitary && deficiency > 0 &&
            throw(ArgumentError("unitary basis cannot be rank-deficient"))
        sizes = Integers(mindim, maxdim)
        return new{T}(mindim, maxdim, sizes, unitary, deficiency)
    end
end

BasisVectors(mindim, maxdim; kwargs...) = BasisVectors{Float64}(mindim, maxdim; kwargs...)

function Base.:(==)(bv1::BasisVectors, bv2::BasisVectors)
    return bv1.sizes == bv2.sizes &&
        bv1.unitary == bv2.unitary &&
        bv1.deficiency == bv2.deficiency
end


function Base.show(io::IO, bv::BasisVectors)
    print(io, typename(bv), "(")
    print(io, "mindim=", string(bv.mindim), ", ")
    print(io, "maxdim=", string(bv.maxdim), "; ")
    print(io, "unitary=", string(bv.unitary), ", ")
    print(io, "deficiency=", string(bv.deficiency), ")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", bv::BasisVectors{T}) where T
    basis = bv.unitary ? "orthogonal basis vectors" : "basis vectors"
    field = T <: Complex ? "complex" : "real"
    sdesc = bv.deficiency == 0 ? "space" : "subspace"
    smin = bv.mindim - bv.deficiency
    smax = bv.maxdim - bv.deficiency

    desc = styled"$basis of a {code:$smin} to {code:$smax} dimensional $field $sdesc"
    if bv.deficiency > 0
        desc = styled"$desc in {code:$(bv.mindim)} to {code:$(bv.maxdim)} dimensions"
    end

    print(io, styled"""
    {code,underline:$BasisVectors}:

        Produce $(desc).
    """)
end

function produce!(tc::TestCase, bv::BasisVectors{T}) where {T}
    n = produce!(tc, bv.sizes)
    n > bv.deficiency || return zeros(T, n, n)

    V = randn(tc.rng, T, n, n)

    if bv.unitary
        Q, R = qr(V)
        V = Q * Diagonal(sign.(diag(R)))
    else
        V ./= transpose(map(norm, eachcol(V)))
    end

    rank(V) == size(V, 2) || reject(tc)

    if bv.deficiency > 0
        ndim = size(V, 2) - bv.deficiency
        for i in 1:ndim
            w = rand(ndim)
            w /= sum(w)
            v = sum(V[:, 1:ndim] .* w', dims=2)
            v /= norm(v)
            V[:, end-ndim+i] = v
        end
    end

    return V
end


struct SquareMatrices{
    T <: Union{AbstractFloat, Complex},
    S <: Union{RealNumbers, ComplexNumbers},
    B <: BasisVectors{T},
} <: Possibility{Matrix{T}}
    spectrum::S
    basis::B

    function SquareMatrices{T}(spectrum, basis) where {T}
        new{T, typeof(spectrum), typeof(basis)}(spectrum, basis)
    end
end

function SquareMatrices(spectrum::Union{RealNumbers, ComplexNumbers}, args...; kwargs...)
    T = postype(spectrum)
    basis = BasisVectors{T}(args...; kwargs...)
    return SquareMatrices{T}(spectrum, basis)
end

function SquareMatrices{T}(
    spectrum::Union{RealNumbers, ComplexNumbers},
    mindim::Integer,
    maxdim::Integer;
    kwargs...
) where {T}
    basis = BasisVectors{T}(mindim, maxdim; kwargs...)
    return SquareMatrices(spectrum, basis)
end

function Base.:(==)(sm1::SquareMatrices, sm2::SquareMatrices)
    return sm1.spectrum == sm2.spectrum && sm1.basis == sm2.basis
end

function Base.show(io::IO, sm::SquareMatrices)
    print(io, typename(sm), "(")
    print(io, ";")
    show(sm.spectrum)
    show(sm.basis)
    print(")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", sm::SquareMatrices)
    print(io, styled"""
    {code,underline:$SquareMatrices}:

        Produce square matrix with:
            * spectrum given by {code:$(sm.spectrum)}
            * basis given by {code:$(sm.basis)}
    """)
end

function produce!(tc::TestCase, sm::SquareMatrices{T,S}) where {T,S}
    V = produce!(tc, sm.basis)
    Vh = sm.basis.unitary ? V' : pinv(V)
    λ = Vector{eltype(T)}(undef, size(V, 2))
    for i in eachindex(λ)
        λ[i] = produce!(tc, sm.spectrum)
    end
    X = V * Diagonal(λ) * Vh
    if sm.basis.unitary
        X = Matrix(eltype(postype(sm)) <: Real ? Symmetric(X) : Hermitian(X))
    end
    # Ensure real matrix is generated when expected
    if T <: Real && eltype(X) <: Complex && all(imag.(X) ≈ 0)
        X = real.(X)
    end
    # Reject if spectrum is not real when expected to be real
    S <: RealNumbers && eltype(eigvals(X)) <: Complex && reject(tc)
    # Reject if basis rank is not correct
    basis_rank = size(X, 2) - sm.basis.deficiency
    rank(eigvecs(X)) == basis_rank || reject(tc)

    return X
end

end

[Seelengrab]

Thanks! At a glance, I'm suspecting the problem is that you're accessing the internals of TestCase (the RNG object) here - removing that and just using a regular randn call with the default RNG may help. I'll be able to take a closer look later today, could you also provide an example property that behaves wrongly?

As an aside, things like manually creating a Vector and filling it in a loop will shrink much worse than using the builtin Vectors{T} :)

[sztal]

Thanks a lot for the hints! I think I managed to reproduce the error in a rather simple test. Again, sorry for a massive code dump, but I wanted to include the entire reprex in one code snippet.

using Test
using LinearAlgebra
using Supposition, Supposition.Data


# POSSIBILITIES FOR GENERATING REAL AND COMPLEX NUMBERS ----------------------------------

module _Numbers
export RealNumbers, ComplexNumbers

import Supposition: example, TestCase, reject
import Supposition.Data: Possibility, Floats, produce!, postype
import StyledStrings: @styled_str

typename(obj) = String(typeof(obj).name.name)


struct RealNumbers{
    T <: AbstractFloat
} <: Possibility{T}
    floats::Floats{T}
    rtol::T

    function RealNumbers{T}(minimum=-T(Inf), maximum=+T(Inf); rtol=nothing) where {T}
        rtol = isnothing(rtol) ? sqrt(eps(T)) : rtol
        rtol > eps(T) ||
            throw(ArgumentError("'rtol' must be greater 'eps($T)'"))
        minimum < maximum ||
            throw(ArgumentError("'minimum' must be lower than 'maximum'"))
        floats = Floats{T}(; minimum=minimum, maximum=maximum, infs=false, nans=false)
        return new{T}(floats, rtol)
    end
end

RealNumbers(args...; kwargs...) = RealNumbers{Float64}(args...; kwargs...)

function Base.propertynames(r::RealNumbers, private::Bool=false)
    return (fieldnames(RealNumbers)..., :minimum, :maximum)
end

function Base.getproperty(r::RealNumbers, name::Symbol)
    if name == :minimum
        return r.floats.minimum
    elseif name == :maximum
        return r.floats.maximum
    else
        return getfield(r, name)
    end
end

Base.:(==)(r1::RealNumbers, r2::RealNumbers) = r1.floats == r2.floats

function Base.show(io::IO, r::RealNumbers)
    print(io, typename(r), "(")
    print(io, "; ")
    print(io, "minimum=", string(r.minimum), ", ")
    print(io, "maximum=", string(r.maximum), ")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", r::RealNumbers)
    left = isfinite(r.minimum) ? "[$(r.minimum)" : "(-∞"
    right = isfinite(r.maximum) ? "$(r.maximum)]" : "∞)"
    print(io, styled"""
    {code,underline:$RealNumbers}:

        Produce real number in {code:$(left), $(right)}.
    """)
end

function produce!(tc::TestCase, r::RealNumbers{T}) where {T}
    x = produce!(tc, r.floats)
    !issubnormal(x) || reject(tc)
    x = abs(x) < eps(T) ? zero(T) : x
    if x != 0 && !isnothing(r.rtol)
        abs(x) >= r.rtol || reject(tc)
    end
    return x
end


struct ComplexNumbers{T <: AbstractFloat} <: Possibility{Complex{T}}
    magnitudes::RealNumbers{T}

    function ComplexNumbers{T}(magnitudes::RealNumbers{T}) where {T}
        magnitudes.minimum >= 0 ||
            throw(ArgumentError("'magnitudes' must have non-negative minimum"))
        return new{T}(magnitudes)
    end
end

function ComplexNumbers{T}(; minabs=zero(T), maxabs=T(Inf), kwargs...) where {T<:Real}
    magnitudes = RealNumbers{T}(minabs, maxabs; kwargs...)
    return ComplexNumbers{T}(magnitudes)
end

function ComplexNumbers(; minabs=0.0, maxabs=Inf, kwargs...)
    minabs = convert(AbstractFloat, minabs)
    maxabs = convert(AbstractFloat, maxabs)
    T = promote_type(typeof(minabs), typeof(maxabs))
    return ComplexNumbers{T}(; minabs=minabs, maxabs=maxabs, kwargs...)
end

function Base.:(==)(cn1::ComplexNumbers, cn2::ComplexNumbers)
    return cn1.magnitudes == cn2.magnitudes
end

function Base.show(io::IO, cn::ComplexNumbers)
    print(io, typename(cn), "(")
    show(io, cn.magnitudes)
    print(io, ")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", cn::ComplexNumbers)
    min = cn.magnitudes.minimum
    max = cn.magnitudes.maximum
    print(io, styled"""
    {code,underline:$ComplexNumbers}:

        Produce complex number with magnitude from {code:$min} to {code:$max}.
    """)
end

function produce!(tc::TestCase, cn::ComplexNumbers{T}) where {T}
    z = randn(tc.rng, Complex{T})
    z /= abs(z)
    m = produce!(tc, cn.magnitudes)
    return z * m
end

end


# POSSIBILITIES FOR GENERATING SQUARE MATRICES -------------------------------------------

module _Matrices
export BasisVectors, SquareMatrices

import Supposition: example, TestCase
import Supposition.Data: Possibility, Integers, produce!, postype, reject
import StyledStrings: @styled_str

using LinearAlgebra:
    norm, rank, qr, diag, pinv, eigvecs, eigvals,
    Diagonal, Symmetric, Hermitian

using .._Numbers


typename(obj) = String(typeof(obj).name.name)

struct BasisVectors{
    T <: Union{AbstractFloat, Complex},
} <: Possibility{Matrix{T}}
    mindim::Int
    maxdim::Int
    sizes::Possibility{Int}
    unitary::Bool
    deficiency::Int

    function BasisVectors{T}(mindim, maxdim; unitary=false, deficiency=0) where {T}
        mindim >= 0 ||
            throw(ArgumentError("'mindim' cannot be negative"))
        maxdim >= mindim ||
            throw(ArgumentError("'maxdim' cannot be lower than 'mindim'"))
        deficiency > mindim &&
            throw(ArgumentError("rank 'deficiency' cannot be greater than 'mindim'"))
        unitary && deficiency > 0 &&
            throw(ArgumentError("unitary basis cannot be rank-deficient"))
        sizes = Integers(mindim, maxdim)
        return new{T}(mindim, maxdim, sizes, unitary, deficiency)
    end
end

BasisVectors(mindim, maxdim; kwargs...) = BasisVectors{Float64}(mindim, maxdim; kwargs...)

function Base.:(==)(bv1::BasisVectors, bv2::BasisVectors)
    return bv1.sizes == bv2.sizes &&
        bv1.unitary == bv2.unitary &&
        bv1.deficiency == bv2.deficiency
end


function Base.show(io::IO, bv::BasisVectors)
    print(io, typename(bv), "(")
    print(io, "mindim=", string(bv.mindim), ", ")
    print(io, "maxdim=", string(bv.maxdim), "; ")
    print(io, "unitary=", string(bv.unitary), ", ")
    print(io, "deficiency=", string(bv.deficiency), ")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", bv::BasisVectors{T}) where T
    basis = bv.unitary ? "orthogonal basis vectors" : "basis vectors"
    field = T <: Complex ? "complex" : "real"
    sdesc = bv.deficiency == 0 ? "space" : "subspace"
    smin = bv.mindim - bv.deficiency
    smax = bv.maxdim - bv.deficiency

    desc = styled"$basis of a {code:$smin} to {code:$smax} dimensional $field $sdesc"
    if bv.deficiency > 0
        desc = styled"$desc in {code:$(bv.mindim)} to {code:$(bv.maxdim)} dimensions"
    end

    print(io, styled"""
    {code,underline:$BasisVectors}:

        Produce $(desc).
    """)
end

function produce!(tc::TestCase, bv::BasisVectors{T}) where {T}
    n = produce!(tc, bv.sizes)
    n > bv.deficiency || return zeros(T, n, n)

    V = randn(tc.rng, T, n, n)

    if bv.unitary
        Q, R = qr(V)
        V = Q * Diagonal(sign.(diag(R)))
    else
        V ./= transpose(map(norm, eachcol(V)))
    end

    rank(V) == size(V, 2) || reject(tc)

    if bv.deficiency > 0
        ndim = size(V, 2) - bv.deficiency
        for i in 1:ndim
            w = rand(ndim)
            w /= sum(w)
            v = sum(V[:, 1:ndim] .* w', dims=2)
            v /= norm(v)
            V[:, end-ndim+i] = v
        end
    end

    return V
end


struct SquareMatrices{
    T <: Union{AbstractFloat, Complex},
    S <: Union{RealNumbers, ComplexNumbers},
    B <: BasisVectors{T},
} <: Possibility{Matrix{T}}
    spectrum::S
    basis::B

    function SquareMatrices{T}(spectrum, basis) where {T}
        new{T, typeof(spectrum), typeof(basis)}(spectrum, basis)
    end
end

function SquareMatrices(spectrum::Union{RealNumbers, ComplexNumbers}, args...; kwargs...)
    T = postype(spectrum)
    basis = BasisVectors{T}(args...; kwargs...)
    return SquareMatrices{T}(spectrum, basis)
end

function SquareMatrices{T}(
    spectrum::Union{RealNumbers, ComplexNumbers},
    mindim::Integer,
    maxdim::Integer;
    kwargs...
) where {T}
    basis = BasisVectors{T}(mindim, maxdim; kwargs...)
    return SquareMatrices(spectrum, basis)
end

function Base.:(==)(sm1::SquareMatrices, sm2::SquareMatrices)
    return sm1.spectrum == sm2.spectrum && sm1.basis == sm2.basis
end

function Base.show(io::IO, sm::SquareMatrices)
    print(io, typename(sm), "(")
    print(io, ";")
    show(sm.spectrum)
    show(sm.basis)
    print(")")
    nothing
end

function Base.show(io::IO, ::MIME"text/plain", sm::SquareMatrices)
    print(io, styled"""
    {code,underline:$SquareMatrices}:

        Produce square matrix with:
            * spectrum given by {code:$(sm.spectrum)}
            * basis given by {code:$(sm.basis)}
    """)
end

function produce!(tc::TestCase, sm::SquareMatrices{T,S}) where {T,S}
    V = produce!(tc, sm.basis)
    Vh = sm.basis.unitary ? V' : pinv(V)
    λ = Vector{eltype(T)}(undef, size(V, 2))
    for i in eachindex(λ)
        λ[i] = produce!(tc, sm.spectrum)
    end
    X = V * Diagonal(λ) * Vh
    if sm.basis.unitary
        X = Matrix(eltype(postype(sm)) <: Real ? Symmetric(X) : Hermitian(X))
    end
    # Ensure real matrix is generated when expected
    if T <: Real && eltype(X) <: Complex && all(imag.(X) ≈ 0)
        X = real.(X)
    end
    # Reject if spectrum is not real when expected to be real
    S <: RealNumbers && eltype(eigvals(X)) <: Complex && reject(tc)
    # Reject if basis rank is not correct
    basis_rank = size(X, 2) - sm.basis.deficiency
    rank(eigvecs(X)) == basis_rank || reject(tc)

    return X
end

end


# TEST DEFINITION ------------------------------------------------------------------------

using ._Numbers
using ._Matrices

reals = RealNumbers(-10, 10)
comps = ComplexNumbers(; maxabs=10)

shape = (4, 10)
dreal = SquareMatrices(reals, shape...)
dcomp = SquareMatrices(comps, shape...)
dsymm = map(Symmetric, SquareMatrices(reals, shape...; unitary=true))
dherm = map(Hermitian, SquareMatrices(comps, shape...; unitary=true))

matrices = dreal | dcomp | dsymm | dherm

@check eigen_exp_fail(X=matrices) = begin
    λ, V = eigen(X)
    expeigen = Eigen(exp.(λ), V)
    Matrix(expeigen) ≈ exp(X)
end

VERSION INFO

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 9 5900HX with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
  JULIA_PROJECT = @.
  JULIA_PKG_PRESERVE_TIERED_INSTALLED = true
  JULIA_REVISE = manual

PKG STATUS

  [7d9fca2a] Arpack v0.5.4
⌃ [86223c79] Graphs v1.11.2
⌃ [5078a376] LazyArrays v2.2.1
⌃ [fa8bd995] MetaGraphsNext v0.7.0
  [5a0628fe] Supposition v0.3.5
  [37e2e46d] LinearAlgebra v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [f489334b] StyledStrings v1.11.0
  [8dfed614] Test v1.11.0

[Seelengrab]

Alright, I figured it out! Using event!, we can check what the actual failing input was:

julia> function fail_prop(X)
           λ, V = eigen(copy(X))
           expeigen = Eigen(exp.(λ), V)
           Matrix(expeigen) ≈ exp(X)
       end
fail_prop (generic function with 1 method)

julia> sr = @check function eigen_exp_fail(X=matrices)
           event!("original input", copy(X))
           fail_prop(X)
       end
Found counterexample
  Context: eigen_exp_fail

  Arguments:
      X::Matrix{Float64} = [-14.434131830845924 -30.413084117462976 -30.499355176855126 -30.84403634581846; 10.715916233391175 20.545030288611116 22.551727646613045 23.321387572019976; 26.080774969555716 67.56928960639824 56.20952169565687 53.42343557789945; -27.958679440121188 -68.02462432664903 -59.95817346160535 -58.12739691863602]

  Events:
    original input
        [-4916.139124639685 -8754.808928780447 -15275.216656224593 6045.163021726459; 3634.371546617668 6469.695032425044 11290.184285848105 -4467.9617289796615; 60.948506843172815 108.04731747010656 188.9013980111351 -74.7282018981361; 1415.3261029795176 2515.575229237147 4392.982424800427 -1738.2642825617095]

# input that got recreated, which mysteriously passes?
julia> fail_prop([-14.434131830845924 -30.413084117462976 -30.499355176855126 -30.84403634581846; 10.715916233391175 20.545030288611116 22.551727646613045 23.321387572019976; 26.080774969555716 67.56928960639824 56.20952169565687 53.42343557789945; -27.958679440121188 -68.02462432664903 -59.95817346160535 -58.12739691863602])
true

# original failing input, which fails as expected!
julia> fail_prop([-4916.139124639685 -8754.808928780447 -15275.216656224593 6045.163021726459; 3634.371546617668 6469.695032425044 11290.184285848105 -4467.9617289796615; 60.948506843172815 108.04731747010656 188.9013980111351 -74.7282018981361; 1415.3261029795176 2515.575229237147 4392.982424800427 -1738.2642825617095])
false

When Supposition.jl shows that input in the final report, it recreates the object from scratch, exactly to prevent anything writing to the input data messing with the final output. The fact that this leads to a different object means that my hunch from above about your use of tc.rng being the culprit is very likely the root cause.

We can confirm this by removing all uses of tc.rng and just use the default task-local rng, i.e. removing calls like randn(tc.rng, ...) and just having randn(...) without passing the internal RNG of the TestCase. The default task-local RNG is already taken care of to be reproducible, so with that change we get:

julia> sr = @check db=false function eigen_exp_fail(X=matrices)
           event!("orig", copy(X))
           fail_prop(X)
       end
Found counterexample
  Context: eigen_exp_fail

  Arguments:
      X::Matrix{Float64} = [-86.26338147266239 8.497246492727074 631.4950303762605 -167.10315469984272; 465.82078981826413 -45.88498624734513 -3410.0624018477797 902.3541875421672; 341.16233357119506 -33.60568984079439 -2497.4944701194363 660.8748837717833; 1359.145933302111 -133.88065500905046 -9949.68998181909 2632.838746580092]

  Events:
    orig
        [-86.26338147266239 8.497246492727074 631.4950303762605 -167.10315469984272; 465.82078981826413 -45.88498624734513 -3410.0624018477797 902.3541875421672; 341.16233357119506 -33.60568984079439 -2497.4944701194363 660.8748837717833; 1359.145933302111 -133.88065500905046 -9949.68998181909 2632.838746580092]

julia> fail_prop([-86.26338147266239 8.497246492727074 631.4950303762605 -167.10315469984272; 465.82078981826413 -45.88498624734513 -3410.0624018477797 902.3541875421672; 341.16233357119506 -33.60568984079439 -2497.4944701194363 660.8748837717833; 1359.145933302111 -133.88065500905046 -9949.68998181909 2632.838746580092])
false

Which is exactly what I'd expect :) So in summary: Please don't access internals; they are internal for a reason! You can find the currently supported User-facing API here.

[sztal]

Thanks! Indeed, removing tc.rng usage the examples reported for failed tests are indeed failing.