Add Gaussian Noise Process submodule

.github/workflows/test.yml

@@ -13,6 +13,9 @@ on:
   # Allows you to run this workflow manually from the Actions tab
   workflow_dispatch:
 
+env:
+  PYTHON : python3
+
 jobs:
   test:
     runs-on: ${{ matrix.os }}
@@ -27,6 +30,10 @@ jobs:
 
     steps:
       - uses: actions/checkout@v2
+      - name: Prepare python
+        run: |
+          python3 -m pip install numpy
+          python3 -m pip install scipy
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.julia-version }}

Project.toml

@@ -4,14 +4,21 @@ authors = ["Joseph Broz <[email protected]>"]
 version = "0.4.2"
 
 [deps]
+DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 
 [compat]
@@ -27,7 +34,6 @@ julia = "1.3"
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
-StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

src/IonSim.jl

@@ -1,3 +1,4 @@
+__precompile__() # this module is safe to precompile
 module IonSim
 
 using QuantumOptics
@@ -33,6 +34,35 @@ An abstract type for specialized bases, which are unique to IonSim.
 abstract type IonSimBasis <: Basis end
 export IonSimBasis
 
+module noise
+using FFTW
+using LinearAlgebra: dot
+using IterativeSolvers: lsqr, lsqr!
+using Random
+using StaticArrays
+
+using PyCall
+const so = PyNULL()
+const np = PyNULL()
+function __init__()
+    #TODO: figure out if this fails gracefully when you don't have scipy/numpy installed
+    #probably should let you use all the code except the noisevector stuff
+    copy!(so, pyimport_conda("scipy.optimize", "scipy"))
+    copy!(np, pyimport_conda("numpy", "numpy"))
+    return
+end
+
+include("noise/noisevector.jl")
+include("noise/ARMA_NoiseProcess.jl")
+include("noise/ARMA_Structure.jl")
+include("noise/CARMA_NoiseProcess.jl")
+include("noise/CARMA_Structure.jl")
+include("noise/Gaussian_types.jl")
+include("noise/Gaussian_interface.jl")
+
+export GaussianProcess, CARMAProcess, ARMAProcess, ARMA, NoiseVector, calculate_noise!
+end
+
 include("constants.jl")
 include("ions.jl")
 include("vibrational_modes.jl")

src/lasers.jl

@@ -1,25 +1,26 @@
 using .PhysicalConstants: c
-
+using .noise: NoiseVector
 export Laser
 
 """
     Laser(;λ=nothing, E=0, Δ=0, ϵ=(x̂+ŷ)/√2, k=ẑ, ϕ=0, pointing::Array{Tuple{Int,Real}})
-        
+
 The physical parameters defining laser light.
 **args**
 * `λ::Union{Real,Nothing}`: the wavelength of the laser in meters
-* `E::Union{Function,Real}`: magnitude of the E-field in V/m
+* `E::Union{Function,Real, NoiseVector}`: magnitude of the E-field in V/m
 * `Δ`: static detuning from f = c/λ in [Hz]
 * `ϵ::NamedTuple`: (ϵ.x, ϵ.y, ϵ.z), polarization direction, requires norm of 1
 * `k::NamedTuple`: (k.x, k.y, k.z), propagation direction, requires norm of 1
-* `ϕ::Union{Function,Real}`: time-dependent phase. of course, this can also be used to model a 
+* `ϕ::Union{Function,Real}`: time-dependent phase. of course, this can also be used to model a
     time-dependent detuning. Units are in radians. Note: if this is set to a function of time,
     then when constructing a Hamiltonian with the `hamiltonian` function, the units of time
     will be as specified by the `timescale` keyword argument.
 * `pointing`: an array of `Tuple{Int,Real}` for describing ion-laser pointing configuration.
     (first element of the tuple is the index for an ion and the second element is the scaling
     factor for the laser's Efield which must be between 0 and 1).
 """
+
 mutable struct Laser
     λ::Union{Real, Nothing}
     E::Function
@@ -36,7 +37,10 @@ mutable struct Laser
         k = ẑ,
         ϕ::Tϕ = 0,
         pointing = Array{Tuple{Int, <:Real}}(undef, 0)
-    ) where {TE, Tϕ}
+    ) where {
+        TE <: Union{NoiseVector, Real, Function},
+        Tϕ <: Union{NoiseVector, Real, Function}
+    }
         rtol = 1e-6
         @assert isapprox(norm(ϵ), 1, rtol = rtol) "!(|ϵ| = 1)"
         @assert isapprox(norm(k), 1, rtol = rtol) "!(|k| = 1)"
@@ -50,8 +54,22 @@ mutable struct Laser
         for s in scaling
             @assert 0 <= s <= 1 "must have s ∈ [0,1]"
         end
-        TE <: Number ? Et(t) = E : Et = E
-        Tϕ <: Number ? ϕt(t) = ϕ : ϕt = ϕ
+        if TE <: Number
+            Et = t -> E
+        elseif TE <: NTuple{2, Vector}
+            Et = t -> NoiseVector(E[1], E[2])
+        else
+            Et = E
+        end
+
+        if Tϕ <: Number
+            ϕt = t -> ϕ
+        elseif TE <: NTuple{2, Vector}
+            ϕt = t -> NoiseVector(ϕ[1], ϕ[2])
+        else
+            ϕt = ϕ
+        end
+
         return new(λ, Et, Δ, ϵ, k, ϕt, pointing)
     end
     # for copying
@@ -82,12 +100,12 @@ function Base.setproperty!(L::Laser, s::Symbol, v::Tv) where {Tv}
         @assert isapprox(norm(v), 1, rtol = rtol) "!(|ϵ| = 1)"
         # if ! isapprox(ndot(L.k, v), 0, rtol=rtol)
         #     @warn "!(ϵ ⟂ k)"
-        # end 
+        # end
     elseif s == :k
         @assert isapprox(norm(v), 1, rtol = rtol) "!(|k| = 1)"
         # if ! isapprox(ndot(v, L.ϵ), 0, rtol=rtol)
         #     @warn "!(ϵ ⟂ k)"
-        # end 
+        # end
     elseif s == :pointing
         b = Tv <: Vector{Tuple{Int64, Float64}} || Tv <: Vector{Tuple{Int64, Int64}}
         @assert b "type != Vector{Tuple{Int,Real}}"
@@ -99,7 +117,13 @@ function Base.setproperty!(L::Laser, s::Symbol, v::Tv) where {Tv}
             @assert 0 <= s <= 1 "must have s ∈ [0,1]"
         end
     elseif s == :E || s == :ϕ
-        Tv <: Number ? vt(t) = v : vt = v
+        if Tv <: Number
+            vt = t -> v
+        elseif Tv <: NTuple{2, Vector}
+            vt = t -> NoiseVector(v[1], v[2])
+        else
+            vt = v
+        end
         Core.setproperty!(L, s, vt)
         return
     end

src/noise/ARMA_NoiseProcess.jl

@@ -0,0 +1,186 @@
+using DiffEqNoiseProcess
+using LinearAlgebra: dot
+using Random
+using PyCall
+
+@inline wiener_randn(rng::AbstractRNG, ::Type{T}) where {T} = randn(rng, T)
+function arma_step!(x, mean, sigma, rng, T)
+    return x[1] = mean + sigma * wiener_randn(rng, T)
+end
+
+# == Generators of Discrete AR/MA time series == #
+mutable struct AR_Process{T1, T2, T3, T4, T5}
+    phi::T1
+    p::T2
+    sigma::T3
+    mean::T3
+    c::T3
+    dt::T3
+    t::T3
+    filter::T4
+    past::T5
+end
+
+mutable struct MA_Process{T1, T2, T3, T4, T5}
+    psi::T1
+    q::T2
+    sigma::T3
+    mean::T3
+    dt::T3
+    t::T3
+    et::T4
+    filter::T5
+end
+
+@inline function (X::AR_Process)(dW, W, dt, u, p, t, rng)
+    T = (typeof(dW) <: AbstractArray) ? dW : typeof(dW)
+
+    # == derivative of step function == #
+    ϵ = t isa Union{Rational, Integer} ? 0 : 100eps(typeof(t))
+    t += dt
+    if t < X.t + ϵ
+        return 0 #dW = 0 between steps
+    end
+
+    # == AR(p) Step == #
+    while t >= X.t + ϵ
+        mean = X.c + dot(X.phi, X.past)
+        pushfirst!(X.past, 0.0)
+        pop!(X.past)
+
+        arma_step!(X.past, mean, X.sigma, rng, T)
+        while X.filter(X.past[1] - X.mean) #Filter Function
+            arma_step!(X.past, mean, X.sigma, rng, T)
+        end
+
+        X.t += X.dt
+    end
+
+    return X.past[1] * max(dt, X.dt)
+end
+
+@inline function (X::MA_Process)(dW, W, dt, u, p, t, rng)
+    T = (typeof(dW) <: AbstractArray) ? dW : typeof(dW)
+
+    # == derivative of step function == #
+    ϵ = t isa Union{Rational, Integer} ? 0 : 100eps(typeof(t))
+    t += dt
+    if t < X.t + ϵ
+        return 0 #dW = 0 between steps
+    end
+
+    # == MA(∞) Step == #
+    Xt = 0
+    while t >= X.t + ϵ
+        pushfirst!(X.et, 0.0) #add new first value
+        pop!(X.et) #remove last value
+        arma_step!(X.et, 0.0, X.sigma, rng, T)
+        Xt = X.mean + dot(X.psi, X.et)
+
+        while X.filter(Xt - X.mean) #Filter Function
+            arma_step!(X.et, 0, X.sigma, rng, T)
+            Xt = X.mean + dot(X.psi, X.et)
+        end
+
+        X.t += X.dt
+    end
+    return Xt * max(dt, X.dt)
+end
+
+# == Bridge Function == #
+@inline function AR_STEP_BRIDGE(X, dW, W, W0, Wh, q, h, u, p, t, rng)
+    # setup_next_step! correction
+    if isapprox(W0, W.curW) # line 97
+        return AR_STEP_BRIDGE(X, dW, W, 1, Wh, q, h, u, p, t, rng)
+    end
+
+    if t isa Union{Rational, Integer}
+        i = searchsortedlast(W.t, h * floor(t / h))
+    else
+        i = searchsortedlast(W.t, h * floor(t / h) + eps(typeof(t)))
+    end
+    return W[i] - (1 - q) * W0
+end
+
+# reject_step! correction
+# https://github.com/SciML/DiffEqNoiseProcess.jl/blob/master/src/noise_interfaces/noise_process_interface.jl
+@inline function AR_STEP_BRIDGE(X, dW, W, W0::Int, Wh, q, h, u, p, t, rng)
+
+    # line 203 correction
+    # q = dtnew / W.dt -> h = W.dt (not dtnew)
+    if isempty(W.S₂) && isapprox(h, q * W.dt) # if h=dtnew
+        h = W.dt
+    end
+
+    tnew = t + q * h
+    tfloor = (
+        X.dt * floor(round(t / X.dt, digits = 8)),
+        X.dt * floor(round(tnew / X.dt, digits = 8))
+    )
+    ϵ = t isa Union{Rational, Integer} ? 0 : sqrt(eps(typeof(t)))
+
+    while tnew < (X.t - X.dt) + ϵ
+        X.t -= X.dt
+        popfirst!(X.past)
+        push!(X.past, 0.0)
+        # X.past .= [X.past[2:end]; 0.0]
+    end
+
+    # both in same step
+    if isapprox(tfloor[1], tfloor[2])
+        dWnew = 0.0
+        # both in different steps
+    else
+        dWnew = X.past[1] * min(q * h, X.dt)
+    end
+
+    return dWnew
+end
+
+@inline function MA_STEP_BRIDGE(X, dW, W, W0, Wh, q, h, u, p, t, rng)
+    # setup_next_step! correction
+    if isapprox(W0, W.curW) # line 97
+        return MA_STEP_BRIDGE(X, dW, W, 1, Wh, q, h, u, p, t, rng)
+    end
+
+    if t isa Union{Rational, Integer}
+        i = searchsortedlast(W.t, h * floor(t / h))
+    else
+        i = searchsortedlast(W.t, h * floor(t / h) + eps(typeof(t)))
+    end
+    return W[i] - (1 - q) * W0
+end
+
+# reject_step! correction
+# https://github.com/SciML/DiffEqNoiseProcess.jl/blob/master/src/noise_interfaces/noise_process_interface.jl
+@inline function MA_STEP_BRIDGE(X, dW, W, W0::Int, Wh, q, h, u, p, t, rng)
+
+    # line 203 correction
+    # q = dtnew / W.dt -> h = W.dt (not dtnew)
+    if isempty(W.S₂) && isapprox(h, q * W.dt) # if h=dtnew
+        h = W.dt
+    end
+
+    tnew = t + q * h
+    tfloor = (
+        X.dt * floor(round(t / X.dt, digits = 8)),
+        X.dt * floor(round(tnew / X.dt, digits = 8))
+    )
+    ϵ = t isa Union{Rational, Integer} ? 0 : sqrt(eps(typeof(t)))
+
+    while tnew < (X.t - X.dt) + ϵ
+        X.t -= X.dt
+        popfirst!(X.et) # remove last innovation
+        push!(X.et, 0.0) # add zero to end
+    end
+
+    # both in same step
+    if isapprox(tfloor[1], tfloor[2])
+        dWnew = 0.0
+        # both in different steps
+    else
+        dWnew = (X.mean + dot(X.psi, X.et)) * min(q * h, X.dt)
+    end
+
+    return dWnew
+end