Add gaussian random noise

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+GaussianRandomFields = "e4b2fa32-6e09-5554-b718-106ed5adafe9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -15,8 +16,8 @@ YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 julia = "1.3"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "HDF5"]

extra/Plot_tdac_output.ipynb

@@ -27,7 +27,8 @@
     "nx = 200\n",
     "ny = 200\n",
     "N = nx * ny\n",
-    "dims = (nx,ny)"
+    "dims = (nx,ny)\n",
+    "plt.rcParams[\"figure.figsize\"] = (9,9)"
    ]
   },
   {

src/TDAC.jl

@@ -1,6 +1,6 @@
 module TDAC
 
-using Random, Distributions, Statistics, Distributed, Base.Threads, YAML
+using Random, Distributions, Statistics, Distributed, Base.Threads, YAML, GaussianRandomFields
 
 export tdac, main
 
@@ -125,10 +125,72 @@ function resample!(state_resampled::AbstractMatrix{T}, state::AbstractMatrix{T},
 
 end
 
-# Add (0,1) normal distributed white noise to state
-function add_noise!(state, amplitude)
+function get_axes(nx::Int, ny::Int, dx::Real, dy::Real)
+
+    x = range(0, length=nx, step=dx)
+    y = range(0, length=ny, step=dy)
+
+    return x,y
+end
+
+# Initialize a gaussian random field generating function using the Matern covariance kernel
+# and circulant embedding generation method
+# TODO: Could generalise this
+function init_gaussian_random_field_generator(lambda::T,
+                                              nu::T,
+                                              sigma::T,
+                                              x::AbstractVector{T},
+                                              y::AbstractVector{T},
+                                              pad::Int) where T
+
+    # Let's limit ourselves to two-dimensional fields
+    dim = 2
+
+    cov = CovarianceFunction(dim, Matern(lambda, nu, σ = sigma))
+    grf = GaussianRandomField(cov, CirculantEmbedding(), x, y, minpadding=pad)
+
+end
+
+# Get a random sample from gaussian random field grf using random number generator rng
+function sample_gaussian_random_field!(field::AbstractVector{T},
+                                       grf::GaussianRandomFields.GaussianRandomField,
+                                       rng::Random.AbstractRNG) where T
+
+    field .= @view(GaussianRandomFields.sample(grf, xi=randn(rng, randdim(grf)))[:])
+
+end
+
+# Get a random sample from gaussian random field grf using random_numbers
+function sample_gaussian_random_field!(field::AbstractVector{T},
+                                       grf::GaussianRandomFields.GaussianRandomField,
+                                       random_numbers::AbstractVector{T}) where T
+
+    field .= @view(GaussianRandomFields.sample(grf, xi=random_numbers)[:])
+
+end
+
+# Add a gaussian random field to each variable in the state vector of one particle
+function add_random_field!(state::AbstractVector{T},
+                           grf::GaussianRandomFields.GaussianRandomField,
+                           rng::Random.AbstractRNG,
+                           nvar::Int,
+                           dim_grid::Int) where T
+
+    random_field = Vector{Float64}(undef, dim_grid)
 
-    @. state += randn() * amplitude
+    for ivar in 1:nvar
+
+        sample_gaussian_random_field!(random_field, grf, rng)
+        @view(state[(nvar-1)*dim_grid+1:nvar*dim_grid]) .+= random_field
+
+    end
+
+end
+
+# Add a (0,1) normal distributed random number, scaled by amplitude, to each element of vec
+function add_noise!(vec::AbstractVector{T}, rng::Random.AbstractRNG, amplitude::T) where T
+
+    @. vec += amplitude * randn((rng,), T)
 
 end
 
@@ -164,6 +226,12 @@ function tdac(params)
     state, state_true, state_avg, state_resampled, weights, obs_real, obs_model, ist, jst = init_tdac(params.dim_state,
                                                                                                       params.nobs,
                                                                                                       params.nprt)
+
+    x, y = get_axes(params.nx, params.ny, params.dx, params.dy)
+
+    background_grf = init_gaussian_random_field_generator(params.lambda, params.nu, params.sigma, x, y, params.padding)
+
+    rng = Random.MersenneTwister(params.random_seed)
 
     # Set up tsunami model
     gg, hh, hm, hn, fm, fn, fe = LLW2d.setup(params.nx, params.ny, params.bathymetry_setup)
@@ -172,6 +240,9 @@ function tdac(params)
     eta = reshape(@view(state_true[1:params.dim_grid]), params.nx, params.ny)
     LLW2d.initheight!(eta, hh, params.dx, params.dy, params.source_size)
 
+    # Initialize all particles to the true initial state
+    state .= state_true
+
     # set station positions
     LLW2d.set_stations!(ist,
                         jst,
@@ -207,7 +278,9 @@ function tdac(params)
         Threads.@threads for ip in 1:params.nprt
 
             tsunami_update!(@view(state[:,ip]), params.nx, params.ny, params.dx, params.dy, params.dt, hm, hn, fn, fm, fe, gg)
+            add_random_field!(@view(state[:,ip]), background_grf, rng, params.n_state_var, params.dim_grid)
             get_obs!(@view(obs_model[:,ip]), @view(state[:,ip]), params.nx, ist, jst)
+            add_noise!(@view(obs_model[:,ip]), rng, params.obs_noise_amplitude)
 
         end
 

src/params.jl

@@ -30,15 +30,22 @@ Parameters for TDAC run. Arguments:
 * `inv_rr` : Inverse of length scale, stored for performance
 * `source_size` : Initial condition parameter
 * `bathymetry_setup` : Bathymetry set-up.
+* `lambda` : Length scale for Matérn covariance kernel (could be same as rr)
+* `nu` : Smoothess parameter for Matérn covariance kernel
+* `sigma` : Marginal standard deviation for Matérn covariance kernel
+* `padding` : Min padding for circulant embedding gaussian random field generator
+* `obs_noise_amplitude`: Multiplier for noise added to observations of the true state
+* `random_seed` : Seed number for the pseudorandom number generator
 
 """
 Base.@kwdef struct tdac_params{T<:AbstractFloat}
 
     nx::Int = 200
     ny::Int = 200
-
+
+    n_state_var::Int = 3
     dim_grid::Int = nx * ny
-    dim_state::Int = 3 * dim_grid
+    dim_state::Int = n_state_var * dim_grid
     dx::T = 2.0e3
     dy::T = 2.0e3
 
@@ -63,6 +70,14 @@ Base.@kwdef struct tdac_params{T<:AbstractFloat}
 
     source_size::T = 3.0e4
     bathymetry_setup::T = 3.0e4
+
+    lambda::T = 1.0e4
+    nu::T = 2.5
+    sigma::T = 1.0
+    padding::Int = 100
+    obs_noise_amplitude::T = 1.0
+
+    random_seed::Int = 12345
 end
 
 end

test/runtests.jl

@@ -1,5 +1,5 @@
 using TDAC
-using LinearAlgebra, Test, HDF5
+using LinearAlgebra, Test, HDF5, Random
 
 @testset "LLW2d" begin
     using TDAC.LLW2d
@@ -126,6 +126,15 @@ end
     TDAC.tsunami_update!(x, nx, ny, dx, dy, dt, hm, hn, fm, fn, fe, gg)
     @test sum(eta, dims=1) ≈ [0.9140901416339269 1.7010577375770561 0.9140901416339269 0.06356127284539884 0.0 0.0 0.0 0.0 0.0 0.0]
     @test sum(eta, dims=2) ≈ [0.9068784611641829; 1.6999564781646717; 0.9204175965604575; 0.06554675780099671; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
+
+    # Test gaussian random field sampling
+    x = 1.:2.
+    y = 1.:2.
+    grf = TDAC.init_gaussian_random_field_generator(1.0,1.0,1.0,x,y,0)
+    f = zeros(4)
+    rnn = [9.,9.,9.,9.]
+    TDAC.sample_gaussian_random_field!(f,grf,rnn)
+    @test f ≈ [17.07868874658702, 3.812477565296764, 3.8124775652967635, 1.8023374615813794]
 end
 
 @testset "TDAC integration tests" begin