I am writing unit tests for Matlab and Julia functions (Julia translated from Matlab), and both use the randn() function. Is there a way to get the same sequence of random numbers both in Matlab and Julia?
Without that it is difficult to write unit tests that give the same results for both languages.
Example of a function I want to test:
"""
getWindDirT(::Direction_Constant_wErrorCov, WindDir, iT)
Return wind direction in SOWFA-deg for the requested turbine(s).
# Arguments
- `WindDir`: Struct with fields
- `Data::Float64`: wind direction value
- `CholSig::AbstractMatrix`: Cholesky factor of covariance matrix (nT x nT)
- `iT`: Vector of turbine indices (can be any indexable collection)
# Returns
- `phi`: Vector of wind directions for the selected turbines, including random perturbation
"""
function getWindDirT(::Direction_Constant_wErrorCov, WindDir, iT)
n = length(iT)
phi = fill(WindDir.Data, n)
# randn(n) gives a vector of n normal random numbers
# WindDir.CholSig should
@assert issquare(WindDir.CholSig)
phi .= phi .+ WindDir.CholSig * randn(n)
return phi
end
Unit test in Julia:
using Random
Random.seed!(1234)
dir_mode = Direction_Interpolation_wErrorCov()
# Example wind direction data (time, phi)
wind_data = [
0.0 10.0;
5.0 20.0;
10.0 30.0
]
# Example Cholesky factor (for 2 turbines)
chol_sig = cholesky([1.0 0.5; 0.5 1.0]).L
# Create WindDir instance
WindDir = WindDirMatrix(wind_data, chol_sig)
# Example turbine indices (for 2 turbines)
iT = [1, 2]
# Example time
t = 7.0
# Call the function
phi = getWindDirT(dir_mode, WindDir, iT, t)
@test size(phi) == (2,1)
@test phi[1] ā 24.92903352636283
@test phi[2] ā 24.363944731838128
But if I translate this unit test to Matlab it will not work, because the sequence of random numbers is different.
const DNG = Iterators.Stateful(parse(Float64, x) for x in eachline("rng-output.csv"))
randn() = popfirst!(DNG)
randn(n) = collect(Iterators.take(DNG, n))
Edit: oh, but thatāll only work if your randn calls are entirely in the scope of your unit test, which it looks like theyāre called from your source. In such a case, I think youād want to explicitly pass a RNG object instead of using the global oneā¦ and then youād need to do a bit more work to satisfy the RNG interface. Thatās more complicated, but I think there may be some existing RNG types that could do the ticket?
I would wager that the exact method of going from the MersenneTwisterās ānativeā output (32 or 52 or 64 raw bits) to a normal distribution is going to be sufficiently complicated to not match. Saving the outputs to a one-column CSV would likely be simplest. For the purposes of unit-testing, you could probably get away with something like this:
struct DeterministicNumberGenerator{T}
iter::T
end
DeterministicNumberGenerator(filename::String) = DeterministicNumberGenerator(Iterators.Stateful(parse(Float64, x) for x in eachline(filename)))
Random.randn(dng::DeterministicNumberGenerator) = popfirst!(dng.iter)
Random.randn(dng::DeterministicNumberGenerator, n) = collect(Iterators.take(dng.iter, n))
And then making your signatures accept an rng argument, calling randn(rng, n).
Well, it might be possible, but I want to have the unit tests of our Julia package FLORIDyn.jl running on Github without access to Matlab. So this is not a good option. The other way round it would be easier, but still more complicated than just sharing a file with values.
IMO if the test depends on randomness, the āfix seed and use approxā combo is unreliable and might miss edge cases (what if all other seeds break everything?).
You could follow a statistical approach:
In Matlab, run the testās code with different seeds many times (1000x, the more the better). Modify the code to compute averages of quantities of interest (phi in your case) or write out all 1000 phis to a file.
In Julia, run the code the same number of times, again computing the mean. Suppose phi_mean holds the mean of the vector phi. Then, @test abs(phi_mean[1] - phi_mean_matlab[1]) < 0.01. Same for the other component. The 0.01 here is an arbitrary small number that somehow depends on the number of samples (more samples => could use 0.005 as the threshold).
Ideally, put the entire sample computed by Matlab in the file with your Julia tests and test the hypothesis E(phi_Julia[1] - phi_Matlab[1]) == 0 (same for phi[2]) with a two-sample test using HypothesisTests.jl.
Since thereās a Normal distribution involved, perhaps itās possible to compute the expected value exactly. Then you donāt even need Matlab: simply test whether phi_mean is sufficiently close to that true expectation.
Octave is easier to access though, and thereās some code here to patch over initialisation differences between it and MATLAB, and have Octave produce a MATLAB compatible sequence. So thatās another option worth keeping in mind.
Not much to add to the discussion other than that I asked this exact question over a decade ago now, but in the context of Julia and R. I got a response from Dirk Eddelbuettel here.
I am now using this Matlab script to create a sequence of random numbers:
clear all
rng(1234);
N=1000;
vec=randn(1,N);
filename = 'test/randn.mat';
save(filename, 'vec')
And this code to create a file based pseudo-random number generator in Julia:
# Copyright (c) 2025 Uwe Fechner
# SPDX-License-Identifier: BSD-3-Clause
using FLORIDyn
using Test
using LinearAlgebra
using MAT
using Random
if basename(pwd()) == "test"
cd("..")
end
filename="test/randn.mat"
randn_vec=vec(matread(filename)["vec"])
mutable struct FileRNG <: AbstractRNG
const data::Vector{Float64}
idx::Int
end
FileRNG(data::Vector{Float64}) = FileRNG(data, 1)
function Base.rand(rng::FileRNG)
x = rng.data[rng.idx]
rng.idx += 1
return x
end
function Base.rand(rng::FileRNG, ::Type{Float64})
return rand(rng)
end
function Base.rand(rng::FileRNG, ::Type{T}) where {T}
error("FileRNG only supports Float64")
end
Random.randn(rng::FileRNG) = rand(rng)
rng = FileRNG(randn_vec)
FLORIDyn.set_rng(rng)
include("test_windfield.jl")
In my package I defined:
# global variables
RNG::AbstractRNG = Random.default_rng()
function set_rng(rng)
global RNG
RNG = rng
end
This means, normally the default RNG is used, but for the unit tests the FileRNG is used.