Merge pull request #12 from DennisWayo/feature/reproducible-sine-noise-padding

Add reproducible sine noise and robust signal padding

Author: Gauthameshwar

src/signals/SignalConverters.jl

@@ -15,18 +15,20 @@ export signal_mps, signal_ztmps
 # Convert an input 1D array signal into ITensor MPS with the binary encoding of the indices
 function _array_to_tensor(x::AbstractVector; sites=undef)
     N = length(x)
-    n = round(Int, log2(N)) # no. of qubits that can encodes signal
+    N > 0 || throw(ArgumentError("_array_to_tensor: Input signal must be non-empty."))
+    n = ceil(Int, log2(N)) # no. of qubits that can encode signal
+    padded_length = 2^n
 
     ITensors.disable_warn_order()
     try
-        # If signal is not a power of 2, fill with 0s upto length 2^n with a warning
-        if N < 2^n
-            @warn "_array_to_tensor: Input signal length $N is not a power of 2. Filling with zeros upto length $(2^n). We recommend providing signals of length power-of-2 for best performance."
-            x_filled = zeros(2^n)
+        # If signal is not a power of 2, fill with 0s up to the next power of 2 with a warning.
+        if N != padded_length
+            @warn "_array_to_tensor: Input signal length $N is not a power of 2. Filling with zeros upto length $(padded_length). We recommend providing signals of length power-of-2 for best performance."
+            x_filled = zeros(eltype(x), padded_length)
             x_filled[1:N] .= x
             x = x_filled
         end
-        @assert length(x) == 2^n ||
+        @assert length(x) == padded_length ||
                 "_array_to_tensor: Length of signal vector must be a power of 2"
         if sites === undef
             sites = [Index(2; tags=@sprintf("site-%d", i)) for i in 1:n]

src/signals/Signals.jl

@@ -18,10 +18,12 @@ function _generate_signal(
     freq::Float64;
     phase::Float64=0.0,
     noise_level::Float64=0.0,
+    seed::Union{Nothing,Int}=nothing,
     kwargs...,
 )
     jvals = 0:(2^n-1)
-    return [sin(freq * dt * j + phase) + noise_level * randn() for j in jvals]
+    rng = isnothing(seed) ? Random.default_rng() : Xoshiro(seed)
+    return [sin(freq * dt * j + phase) + noise_level * randn(rng) for j in jvals]
 end
 
 function _generate_signal(
@@ -50,13 +52,15 @@ function _generate_signal(
     freq::Vector{Float64};
     phase::Vector{Float64}=zeros(length(freq)),
     noise_level::Float64=0.0,
+    seed::Union{Nothing,Int}=nothing,
     kwargs...,
 )
     jvals = 0:(2^n-1)
     length(freq) == length(phase) ||
         throw(ArgumentError("Frequency and phase vectors must be of the same length."))
+    rng = isnothing(seed) ? Random.default_rng() : Xoshiro(seed)
     return [
-        sum(sin(ω * dt * j + φ) for (ω, φ) in zip(freq, phase)) + noise_level * randn() for
+        sum(sin(ω * dt * j + φ) for (ω, φ) in zip(freq, phase)) + noise_level * randn(rng) for
         j in jvals
     ]
 end
@@ -162,14 +166,18 @@ Generate a length-`2^n` real signal of a specified type.
    Can be a scalar or a vector of frequencies. Defaults to `2π` if not provided.
 
 # kwargs (kind-dependent)
-- `phase::Float64`: The phase offset of the signal in radians. Applicable to `:sin` and `:cos` kinds. 
-   Defaults to `0.0`.
-- `decay_rate::Float64`: The rate of exponential decay. Only applicable to the `:decay` kind. 
-   A higher value results in faster decay. Defaults to `1.0`.
-- `noise_level::Float64`: The amplitude of random noise added to the signal. Applicable to all kinds. 
-   Defaults to `0.0` (no noise).
-- `seed::Int`: The random seed for reproducibility of noise generation. Applicable when `noise_level > 0`. 
-   Defaults to `nothing` (random seed).
+- `phase`: Phase offset in radians.
+  - For `:sin` with scalar frequency, use `phase::Float64` (default `0.0`).
+  - For `:sin` with vector frequency, use `phase::Vector{Float64}` (defaults to zeros).
+  - For `:sin_decay`, phase is optional (`Float64` or `Vector{Float64}` depending on frequency type).
+- `decay_rate`: Exponential damping rate. Required for `:sin_decay`
+  (`Float64` for scalar-frequency input, `Vector{Float64}` for vector-frequency input).
+- `noise_level::Float64`: The amplitude of Gaussian noise added to `:sin` signals. Defaults to `0.0`
+  (no noise).
+- `seed`: Random seed control:
+  - For `:sin` with `noise_level > 0`, pass `seed::Int` for reproducible noise; default `nothing`
+    (use current global RNG state).
+  - For `:random`, pass `seed::Int` (default `1234`).
 
 # Returns
 - `signal::Vector{Float64}`: A real-valued vector of length `2^n` representing the generated signal.

test/test_signal_converters.jl

@@ -19,6 +19,27 @@ using QILaplace.SignalConverters:
     ψ_padded, norm_padded = _array_to_tensor(x_non_pow2)
     @test length(inds(ψ_padded)) == 2  # rounds up to 4 = 2^2
     @test isapprox(norm_padded, sqrt(1.0 + 4.0 + 9.0); rtol=1e-12)
+
+    # Test with non-power-of-2 requiring ceil(log2(N)) (5 -> 8)
+    x_non_pow2_big = [1.0, 2.0, 3.0, 4.0, 5.0]
+    ψ_padded_big, norm_padded_big = _array_to_tensor(x_non_pow2_big)
+    @test length(inds(ψ_padded_big)) == 3  # rounds up to 8 = 2^3
+    @test isapprox(norm_padded_big, sqrt(sum(abs2, x_non_pow2_big)); rtol=1e-12)
+
+    # Non-power-of-2 complex signal should pad without type errors.
+    x_non_pow2_complex = ComplexF64[1.0 + 1.0im, 2.0 - 1.0im, 3.0 + 0.5im]
+    ψ_padded_complex, norm_padded_complex = _array_to_tensor(x_non_pow2_complex)
+    @test length(inds(ψ_padded_complex)) == 2
+    @test isapprox(norm_padded_complex, sqrt(sum(abs2, x_non_pow2_complex)); rtol=1e-12)
+
+    # Through public API, padded tail coefficients should be exactly zero.
+    ψ_non_pow2 = signal_mps(x_non_pow2_big; method=:svd, cutoff=0.0, maxdim=typemax(Int))
+    for i in 0:4
+        @test isapprox(coefficient(ψ_non_pow2, i), x_non_pow2_big[i+1]; atol=1e-12)
+    end
+    for i in 5:7
+        @test isapprox(coefficient(ψ_non_pow2, i), 0.0; atol=1e-12)
+    end
 end
 
 # ==================== Test _tensor_to_mps_svd ====================

test/test_signals.jl

@@ -78,6 +78,28 @@ end
     r2 = generate_signal(nrand, kind=:random, seed=42)
     @test r1 == r2
 
+    # seeded noise for :sin should be reproducible
+    s_noise_1 = generate_signal(
+        n, kind=:sin, dt=dt, freq=freq, phase=phase, noise_level=0.1, seed=2026
+    )
+    s_noise_2 = generate_signal(
+        n, kind=:sin, dt=dt, freq=freq, phase=phase, noise_level=0.1, seed=2026
+    )
+    s_noise_3 = generate_signal(
+        n, kind=:sin, dt=dt, freq=freq, phase=phase, noise_level=0.1, seed=2027
+    )
+    @test s_noise_1 == s_noise_2
+    @test s_noise_1 != s_noise_3
+
+    # seeded noise for vector-frequency :sin path should also be reproducible
+    s_noise_vec_1 = generate_signal(
+        n, kind=:sin, dt=dt, freq=[1.0, 2.0], phase=[0.0, 0.5], noise_level=0.1, seed=77
+    )
+    s_noise_vec_2 = generate_signal(
+        n, kind=:sin, dt=dt, freq=[1.0, 2.0], phase=[0.0, 0.5], noise_level=0.1, seed=77
+    )
+    @test s_noise_vec_1 == s_noise_vec_2
+
     # integer freq should work like float
     sint = generate_signal(n, kind=:sin, dt=dt, freq=1)
     sfloat = generate_signal(n, kind=:sin, dt=dt, freq=1.0)