nansencenter / DAPPER

Data Assimilation with Python: a Package for Experimental Research
https://nansencenter.github.io/DAPPER
MIT License
341 stars 119 forks source link

Test operator #88

Closed dafeda closed 2 years ago

dafeda commented 2 years ago

When I read code, I tend to write tests in order to get a feel for how things work. Here are a few simple tests I wrote in case you are interested. Also, I thinks some of the lambdas could make do without arguments.

dafeda commented 2 years ago

I believe the args to the lambdas need to stay, because they will be called with those args, for example in ExtKF.

Is the idea that we should be able to run ExtKF with model set to None, something like:

import numpy as np

import dapper.mods as modelling
from dapper.mods.LA import Fmat, homogeneous_1D_cov
from dapper.mods.Lorenz96 import LPs
import dapper.da_methods as da

tseq = modelling.Chronology(dt=1, dko=5, T=300, BurnIn=-1, Tplot=100)
Nx = 10
Fm = Fmat(Nx, -1, 1, tseq.dt)

def step(x, t, dt):
    assert dt == tseq.dt
    return x @ Fm.T

Dyn = {
    "M": Nx,
    #    'model': step,
    "model": None,
    # 'linear': lambda x, t, dt: Fm,
    "noise": 0,
}
Dyn_op = modelling.Operator(**Dyn)

X0 = modelling.GaussRV(mu=np.zeros(Nx), C=homogeneous_1D_cov(Nx, Nx / 8, kind="Gauss"))

Ny = 4
jj = modelling.linspace_int(Nx, Ny)
Obs = modelling.partial_Id_Obs(Nx, jj)
Obs["noise"] = 0.01

HMM = modelling.HiddenMarkovModel(Dyn_op, Obs, tseq, X0, LP=LPs(jj))

xp = da.ExtKF()

xx, yy = HMM.simulate()

xp.assimilate(HMM, xx, yy)

When I try running this I get:

ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)

The tests are nice. Apologies that the code around operator initialization is a bit messy. FYI there's nothing really interesting going on there, just bookkeeping and convenience. But I understand that you want to get a handle on it.

I've have never written non-trivial code that's not messy 😄

patnr commented 2 years ago

The idea was just to default to the identity if you don't specify the observation operator. I.e. it will assume you want to observe (directly) the full state variable. This is a fairly common case, although it might be too much against the "explicit is better than implicit" principle.

But it could also be used for the dynamical operator. In fact, that is what I've don in the latest commit which sets up a test case with effectively no dynamics. This is interesting, and something I've wanted to do for some time, so thank you for the impetus.

Note that there was a bug with one of the lambdas' args. This has now been eliminated with the removal of Id_mat as a whole.

dafeda commented 2 years ago

Thanks for the explanation. I prefer explicit when I don't need much flexibility, but it's not always easy to know how flexible the code should be.

The latest commit is nice! I'll rebase next week and finish these simple tests.

dafeda commented 2 years ago

I've rebased and cleaned up a bit. Great if you would have a look @patricknraanes .

dafeda commented 2 years ago

I don't think test_ExtKF is necessary, since the ExtKF is already being tested in tests/test_example2.py, while the HMMs get tested in test_demos.py and test_HMMs.py

Otherwise it looks all good!

I wrote test_ExtKF in order to catch the bug fixed here: https://github.com/nansencenter/DAPPER/commit/6ee1e3e7d4a6668e560ccb01ee231574546708ff

Neither of tests/test_example2.py, test_demos.py or test_HMMs.py will fail if the fix linked to above is reverted. The newly added test test_operator_defaults() will fail but for other reasons. I am happy to remove test_ExtKF as I don't much like writing tests without explicit asserts, but now at least you are aware why I wrote it in the first place.

patnr commented 2 years ago

Thank you for the PR!