Add gaussian random noise

[tkoskela]

Adding functionality to generate gaussian random fields with the GaussianRandomFields.jl package.

Also initialises all particles to the true initial state.

Adding background noise to particles to address #6

TODO: Add noise to observations to address #17

[tkoskela]

@giordano could you investigate why the CI checks are failing on this one? I do not understand the first error message and google is not helping

[giordano]

I think the error was temporary, just rerunning the job made it go away. However there is another problem with the tests now: we're generating random numbers without fixing the seed (Random.seed!)

[tkoskela]

The random numbers are generated with a fixed seed. I create a random number generator and initialize it with params.random_seed

https://github.com/Team-RADDISH/TDAC.jl/blob/da742c350154f0a45c60c5f1759a6f88e5cf9be7/src/TDAC.jl#L218

It is then used to generate the random numbers for sampling the gaussian random fields

https://github.com/Team-RADDISH/TDAC.jl/blob/da742c350154f0a45c60c5f1759a6f88e5cf9be7/src/TDAC.jl#L159

Do you prefer doing it in a different way?

Ok, I admit the last line looks a bit obscure, maybe I could make it more readable.

[giordano]

Uh, oh... ok, I missed it 😬 I thought that was the reason why the tests are failing, they're failing also for me locally.

With which Julia version did generate the reference numbers for the test? Can you please rebase on master so that the test for Julia v1.3 will run even if the other fails?

[AlexandrosBeskos]

Hi,
let's be careful with the random seed.

1. We definitely do NOT want to use exactly the same random Gaussian field for each particle - these should be independent for each particle.
2. It is good to run the particle filter, or simulate data, using not the same seed all the time, so that we can see the performance of the algorithm for non-identical datasets.
   Thanks, Alex

[tkoskela]

Thanks for comments Alex!

1. The random fields should indeed be independent for each particle. I initially made a mistake where this was not the case, but I fixed it in da742c3. It is good you brought this up because it reminds me I need to address this when merging with MPI Parallelisation #11 to avoid each process having the same set of random numbers.

2. I understand the point, but it is still good to keep control of the random seed in the hands of the user so that our runs are reproducible. The way I did it in this PR is I added the random seed as an input parameter that the user can set for each run. Therefore if you want to run with a different random seed, you simply set it to a different value in the input file.

[AlexandrosBeskos]

Hi,
you did very well with the seed, nothing more is needed there.
Thanks, Alex

[tkoskela]

Hooray, tests pass!

[giordano]

I think this assertion is wrong. Consider for example

julia> using GaussianRandomFields, TDAC

julia> x = 1.:10.;

julia> y = 1.:10.;

julia> grf = TDAC.init_gaussian_random_field_generator(1.0,1.0,1.0,x,y,0);

julia> f = zeros(100);

julia> rnn = randn(100);

julia> TDAC.sample_gaussian_random_field!(f, grf, rnn);
ERROR: DimensionMismatch("length of random points vector must be equal to 1024")
Stacktrace:
 [1] #sample#10 at /home/mose/.julia/dev/GaussianRandomFields/src/gaussian_random_fields.jl:142 [inlined]
 [2] sample_gaussian_random_field!(::Array{Float64,1}, ::GaussianRandomField{CirculantEmbedding,CovarianceFunction{2,Matern{Float64}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},Array{Float64,2},Tuple{Array{Float64,2},FFTW.cFFTWPlan{Complex{Float64},-1,false,2}}}, ::Array{Float64,1}) at /home/mose/.julia/dev/TDAC/src/TDAC.jl:169
 [3] top-level scope at REPL[54]:1
 [4] eval(::Module, ::Any) at ./boot.jl:331
 [5] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [6] run_backend(::REPL.REPLBackend) at /home/mose/.julia/packages/Revise/C272c/src/Revise.jl:1075
 [7] top-level scope at none:0

Changing rnn to have 1024 elements:

julia> rnn = randn(1024);

julia> TDAC.sample_gaussian_random_field!(f, grf, rnn);
ERROR: AssertionError: length(random_numbers) == length(field)
Stacktrace:
 [1] sample_gaussian_random_field!(::Array{Float64,1}, ::GaussianRandomField{CirculantEmbedding,CovarianceFunction{2,Matern{Float64}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},Array{Float64,2},Tuple{Array{Float64,2},FFTW.cFFTWPlan{Complex{Float64},-1,false,2}}}, ::Array{Float64,1}) at /home/mose/.julia/dev/TDAC/src/TDAC.jl:168
 [2] top-level scope at REPL[56]:1
 [3] eval(::Module, ::Any) at ./boot.jl:331
 [4] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [5] run_backend(::REPL.REPLBackend) at /home/mose/.julia/packages/Revise/C272c/src/Revise.jl:1075
 [6] top-level scope at none:0

It's the output of GaussianRandomFields.sample that should have 100 elements, not random_numbers. Not sure how to quickly check the length of sample

[giordano]

Ok, I think the appropriate check would be

prod(length.(grf.pts)) == length(field)

but I feel like I had to dig too much deep into GaussianRandomFields internals to find this. I'd suggest to just remove the check and let broadcast fail naturally if there is a dimension mismatch.

[tkoskela]

Thanks for digging into this! I guess I tried to take a short-cut here and it was wrong.

The appropriate check would be to use what the other instance of sample_gaussian_random_field() uses to get the number of random numbers from the generator, randdim(grf)

The error message from just letting it fail is more helpful though :)

[tkoskela]

By the way, do you think it would be appropriate to use the GaussianRandomFields namespace explicitly when calling functions from it?

[giordano]

do you think it would be appropriate to use the GaussianRandomFields namespace explicitly when calling functions from it?

I personally don't mind too much. The sample name is exported from that module, so we could omit it, but I think it's fine also to qualify the call, to make it easier to understand where it's coming from.