Upon review of the unwrap code this afternoon, I determined that this is due to use of rng() inside of a parallel loop, which leads to race-like out-of-order reliability value generation.

Possible, but I can't reproduce with

julia> measure(s) = ((Random.seed!(10); unwrap(s; dims=1:2)) - (Random.seed!(10); unwrap(s; dims=1:2))).^2 |> mean |> sqrt
measure (generic function with 1 method)

julia> ps = angle.(rand(ComplexF64, 1000, 1000));

julia> measure(ps)
0.0

Be sure to start julia with additional threads - julia -t 8 or somesuch amount. When I start from the command line it usually only starts with one in the pool.

Yes, that was run on a threaded session.

It looks like this is down to the semantics of seed! and the implementation of whatever the default RNG is - I am going to guess the RNG is cloned from its starting state for each soft-thread and then generates the same first value for each pixel, whereas when the rng keyword is used, the RNG is shared and does not. Since unwrap accepts an rng kwarg it should still behave in a stable fashion given some RNG instance, rather than only if used with the task-local RNG.

I cannot reproduce with the default TaskLocalRNG, MersenneTwister, or StableRNG if I do it like this

julia> measure(s, rng, seed) = (unwrap(s; dims=1:2, rng=rng(seed)) - unwrap(s; dims=1:2, rng=rng(seed))).^2 |> mean |> sqrt
measure (generic function with 1 method)

julia> ps = angle.(rand(ComplexF64, 1000, 1000));

julia> measure(ps, StableRNG, 10) == measure(ps, MersenneTwister, 10) == measure(ps, Random.default_rng, 10) == 0
true

However, there is indeed a difference with

julia> ps = rand(1000, 1000) * 1234;

julia> @test measure(ps, StableRNG, 10) == measure(ps, MersenneTwister, 10) == measure(ps, Random.default_rng, 10) == 0
Test Failed at REPL[11]:1
  Expression: measure(ps, StableRNG, 10) == measure(ps, MersenneTwister, 10) == measure(ps, Random.default_rng, 10) == 0
   Evaluated: 294.78516174681107 == 160.76991914304145 == 177.74442713272344 == 0

ERROR: There was an error during testing

Inspecting the output of different runs of unwrap here, I noticed that they differ by multiples of 2pi. The documentation states that unwrap requires the values of the input array to be in a certain range, which is probably why I get these results. Maybe your input is similarly invalid?

Although yes, this is orthogonal to reproducibility, because even with an invalid input we should maybe get the same invalid output...

The theory about invalid input is interesting - I was initially going to suggest something about the nature of the random data before I got distracted with the potentially wrong conclusion about task-local RNG.

Per the DSP.jl docs,

A common usage for unwrapping across a singleton dimension is for a phase measurement over time, such as when comparing successive frames of a short-time-fourier-transform, as each frame is wrapped to stay within (-pi, pi].

A common usage for unwrapping across multiple dimensions is for a phase measurement of a scene, such as when retrieving the phase information of of an image, as each pixel is wrapped to stay within (-pi, pi].

Indeed, the range of my FT phase angle data is within (-pi, pi]

julia> any(-pi .>= ps .> pi)  
false

Anyways, correctness of inputs aside - yes - the behavior should be consistent between applications and means of seeding.

That inequality looks wrong. Indeed,

julia> f(x) = -pi >= x > pi
f (generic function with 1 method)

julia> @code_llvm debuginfo=:none f(100)
; Function Signature: f(Int64)
; Function Attrs: uwtable
define i8 @julia_f_9255(i64 signext %"x::Int64") #0 {
top:
  ret i8 0
}

Could you try instead

all(-pi .< ps .<= pi) # returns false if there is invalid input

Ah, you are correct. I've been staring at the computer for too long.

So, the minimum of the live phase angle data is -pi in accordance with the range of angle (-pi <= angle() <= pi), which does put it out of domain by a small margin.