Initial SDE Code Pull Request

[cgrudz]

Hey Y'all,

Re #70, apologies for delays on this, I am finally getting back around to merging some changes that have been sitting for a while. In this request, I have added the general Runge-Kutta SDE scheme with a fairly simple change to the existing version in the integration schemes in mods. This includes adding an "s" parameter for the instantaneous standard deviation of the noise (the diffusion coefficent) assuming additive noise, and re-writes "order" into a "stages" parameter. This has to due with the discrepancy between the number of stages and the order of the method for the SDE models. This was tested versus the test cases and I found that there was only an issue with a special case in bocquet19, and otherwise the re-naming seems to have been propagated successfully.

I am still working on the full set up of the L96s model on which these changes depend to a certain extent, but I should have this model integrated in the coming weeks as the next stage. This will help me get more into DAPPER development now that I'm preparing my course materials consistently.

Cheers,
Colin

[patnr]

Looks pretty good.

I have marked/suggested/discussed some changes.

One thing that's missing is a script demonstrating these features, ideally reproducing something. Is this coming with the fuller setup?

[patnr]

Suggested change

	The rule has strong / weak convergence order 1.0 for generic SDEs and order 4.0
	convergence for ODEs when stages=4. For stages=1, this becomes the Euler-Maruyama
	schemefor SDEs (s > 0.0) with strong / weak convergence order 1.0 for SDEs with
	additive noise as defined in the below. See `bib.grudzien2020numerical`.
	This scheme has order-4 convergence for ODEs when `stages=4`.
	For SDEs with additive noise as defined below, it has strong / weak convergence
	of order 1.0 (for any value of `stages`). See `bib.grudzien2020numerical`.

Also, when is it strong, and when is it weak ?

Please also add a note saying that the stochastic integrator should not be used if the DA method has its own way of dealing with the noise, such as for gaussian mixtures, particle filters, and the classic/extended KF.

[cgrudz]

Actually regarding this -- how do we set options so that line 244 of DAPPER/dapper/mods/__init__.py

xx[k] = Dyn(xx[k-1], t-dt, dt) + np.sqrt(dt)*Dyn.noise.sample(1)

doesn't add additional noise? I think this may be the only compatibility issue, is the noise a function that we have to supply in the Dyn dict? Should I just set a lambda function to return zero in general?

Otherwise this model configuration can be used with multiplicative or additive inflation / Gaussian mixtures / etc... to handle issues like sampling error or to increase the variance of the estimator. The example that I am writing will be treated as a perfect-random model configuration in which the truth twin and each ensemble member will be propagated by the same random model under different outcomes. Note, the observations are still discrete-in-time so that this is not equivalent to the Kalman-Bucy filter, but rather represents the case in which the truth twin and ensemble members are actually evolved with respect to statistically equivalent, random diffeomorphisms. In this case, many of the usual issues of sampling error due to nonlinearity etc. still apply, though the SDE evolution tends to regularize some stochastic estimators like particle filters which are more degenerate in totally deterministic systems. This also tends to regularize some of the issues with the EnKF's need to use random rotations, etc. to numerically regularize the collapse of variances.

[patnr]

Using HMM = modelling.HiddenMarkovModel(Dyn={..., 'noise': 0}, ...) should do it. But that means that your diffusion/noise must be specified elsewhere. Presumably, for perfect-yet-stochastic-model-case, it should just be inherent to the step function. I believe, for now at least, that is the easiest solution.

[cgrudz]

That makes sense to me -- I'll get a draft of this running shortly.

[patnr]

Actually, maybe you do not need to create dapper/mods/L96s but rather just add a file to dapper/mods/Lorenz96 that is called stoch_integrator.py (for example) containing l96s_tay2_step and nothing more?

Or am I missing something? Please discuss hereunder.

[cgrudz]

The main reason for using this scheme is for the high-accuracy of the truth twin simulation, where this scheme is only designed for the Lorenz-96 model with additive noise that has a scalar covariance at all times. The key is that strong convergence needs to be the criterion for whether a simulation of the truth twin will be consistent with the equations of motion in generating an observation process in a perfect-random model. Because this is a specific sub-class of stochastic Lorenz-96 models (though one that is probably most commonly used in the SDE configuration) I decided to name it the L96-s model where the s refers to the stochasticity, and particularly the diffusion coefficient which defines an "instantaneous" standard deviation of the noise process.

[patnr]

Got you, I think. You still might be able to avoid duplicating Lorenz96.py/extras.py and large parts of __init__.py by importing them from Lorenz96 into L96s, I belive.

[patnr]

Please put the simplest case (s==0) up top. No need for .0 AFAIK.

[cgrudz]

In this case s should generally be a float though, not int. I don't know if it makes a difference in Python because of the type inference / coercion that occurs but this is just a habit from more static typing.

[patnr]

I don't think there is good conformity in the code. But at least for python 3 it should not matter and I tend to lean to using just 0.

[cgrudz]

Definitely want to add the example, the trick is that this requires having two distinct time-steppers for the ensemble and the the truth twin respectively.  Is there any example where this is performed elsewhere in DAPPER?  I'll incorporate the additional suggestions below later today, along with a rough draft of the example and edits I've been working on, add this all to the pull request. Otherwise, strong / weak convergence refers to whether we are referring to convergence in path or convergence in distribution.

…

On 7/15/21 2:02 AM, Patrick N. Raanes wrote: ***@***.**** requested changes on this pull request. Looks good, although I have marked/suggested/discussed some changes. One thing that's missing is a script demonstrating these features, ideally reproducing something. Is this something you are thinking to add? ------------------------------------------------------------------------ In .gitignore <#76 (comment)>: > +############################## +### Vim ### +############################## +*.swp + ############################## IMHO this belongs in |$HOME/.gitignore| ------------------------------------------------------------------------ In dapper/mods/integration.py <#76 (comment)>: > + The rule has strong / weak convergence order 1.0 for generic SDEs and order 4.0 + convergence for ODEs when stages=4. For stages=1, this becomes the Euler-Maruyama + schemefor SDEs (s > 0.0) with strong / weak convergence order 1.0 for SDEs with + additive noise as defined in the below. See `bib.grudzien2020numerical`. ⬇️ Suggested change - The rule has strong / weak convergence order 1.0 for generic SDEs and order 4.0 - convergence for ODEs when stages=4. For stages=1, this becomes the Euler-Maruyama - schemefor SDEs (s > 0.0) with strong / weak convergence order 1.0 for SDEs with - additive noise as defined in the below. See `bib.grudzien2020numerical`. + This scheme has order-4 convergence for ODEs when `stages=4`. + For SDEs with additive noise as defined below, it has strong / weak convergence + of order 1.0 (for any value of `stages`). See `bib.grudzien2020numerical`. Also, when is it strong, and when is it weak ? Please also add a note saying that the stochastic integrator should not be used if the DA method has its own way of dealing with the noise, such as for gaussian mixtures, particle filters, and the classic/extended KF. ------------------------------------------------------------------------ In dapper/mods/L96s/__init__.py <#76 (comment)>: > +######################################################################################## +######################################################################################## Sorry, I want more style conformity. Please use a box-like approach |####################### # This is a heading # ####################### | For smaller sub-headings use |# -------------------------| ------------------------------------------------------------------------ In dapper/mods/L96s/__init__.py <#76 (comment)>: > +def dxdt_autonomous(x): + return (shift(x, 1)-shift(x, -2))*shift(x, -1) - x + + +def dxdt(x): + return dxdt_autonomous(x) + Force + + +def step(x0, t, dt): + return rk4(lambda t, x: dxdt(x), x0, np.nan, dt) + + +######################################################################################## +# 2nd order strong taylor SDE step + +def l96s_tay2_step(x, t, dt): Actually, maybe you do not need to create |dapper/mods/L96s| but rather just add a file to |dapper/mods/Lorenz96| that is called |stoch_integrator.py| (for example) containing |l96s_tay2_step| and nothing more? Or am I missing something? Please discuss hereunder. ------------------------------------------------------------------------ In dapper/mods/integration.py <#76 (comment)>: > + The diffusion coeffient for models with additive noise. Default: 0 for + deterministic integration. ⬇️ Suggested change - The diffusion coeffient for models with additive noise. Default: 0 for - deterministic integration. + The diffusion coeffient for models with additive noise. + Default: 0, yielding deterministic integration. ------------------------------------------------------------------------ In dapper/mods/integration.py <#76 (comment)>: > + The number of stages of the RK method. Default: 4. When stages=1, this becomes + Euler / Euler-Maruyama. ⬇️ Suggested change - The number of stages of the RK method. Default: 4. When stages=1, this becomes - Euler / Euler-Maruyama. + The number of stages of the RK method. + When stages=1, this becomes the Euler (-Maruyama) scheme. + Default: 4. ------------------------------------------------------------------------ In dapper/mods/integration.py <#76 (comment)>: > Returns ------- ndarray State vector at the new time, `t+dt` """ - if order >=1: k1 = dt * f(t , x) # noqa - if order >=2: k2 = dt * f(t+dt/2, x+k1/2) # noqa - if order ==3: k3 = dt * f(t+dt , x+k2*2-k1) # noqa - if order ==4: # noqa - k3 = dt * f(t+dt/2, x+k2/2) # noqa - k4 = dt * f(t+dt , x+k3) # noqa - if order ==1: return x + k1 # noqa - elif order ==2: return x + k2 # noqa - elif order ==3: return x + (k1 + 4*k2 + k3)/6 # noqa - elif order ==4: return x + (k1 + 2*(k2 + k3) + k4)/6 # noqa - else: raise NotImplementedError # noqa + + if s > 0.0: Please put the simplest case (|s==0|) up top. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#76 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABXIPFTFWWCLWFWJEU4CRDTTX2P4BANCNFSM5AMNRCYQ>.

[patnr]

the trick is that this requires
having two distinct time-steppers for the ensemble and the the truth
twin respectively. Is there any example where this is performed
elsewhere in DAPPER?

I think this will do it (in the case of examples/basic_1.py):

@@ -15,6 +15,8 @@ import dapper.da_methods as da
 # #### Load experiment setup: the hidden Markov model (HMM)
 
 from dapper.mods.Lorenz63.sakov2012 import HMM
+HMM2 = HMM.copy()
+HMM2.step = my_custom_step_function
 
 # #### Generate the same random numbers each time this script is run
 
@@ -35,7 +37,7 @@ xp = da.EnKF('Sqrt', N=10, infl=1.02, rot=True)
 
 # #### Assimilate yy, knowing the HMM; xx is used to assess the performance
 
-xp.assimilate(HMM, xx, yy, liveplots=not nb)
+xp.assimilate(HMM2, xx, yy, liveplots=not nb)
 
 # #### Average the time series of various statistics

[patnr]

BTW, I suggest reviewing the changes on GitHub rather than by email. You get pretty coloring, and my up-to-date comments (often edited for tone 🤣 )

[cgrudz]

Ahahahaha sounds good, I'm just working from home at the moment while I take care of some things, once I get to the office later I'll get back into serious dev mode.

…

On 7/15/21 8:40 AM, Patrick N. Raanes wrote: BTW, I suggest reviewing the changes on GitHub rather than by email. You get pretty coloring, and my up-to-date comments (often edited for tone 🤣 ) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#76 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABXIPFVKTBW2FL4YTITU3QLTX36O7ANCNFSM5AMNRCYQ>.

[cgrudz]

Good point, I'll have a look at this to see what redundancy can be eliminated.

…

On 7/15/21 12:22 PM, Patrick N. Raanes wrote: ***@***.**** commented on this pull request. ------------------------------------------------------------------------ In dapper/mods/L96s/__init__.py <#76 (comment)>: > +def dxdt_autonomous(x): + return (shift(x, 1)-shift(x, -2))*shift(x, -1) - x + + +def dxdt(x): + return dxdt_autonomous(x) + Force + + +def step(x0, t, dt): + return rk4(lambda t, x: dxdt(x), x0, np.nan, dt) + + +######################################################################################## +# 2nd order strong taylor SDE step + +def l96s_tay2_step(x, t, dt): Got you, I think. You still might be able to avoid duplicating |Lorenz96.py/extras.py| and large parts of |__init__.py| by importing them from |Lorenz96| into |L96s|, I belive. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#76 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABXIPFUEVVSXVQYJYL7IDSDTX4YQNANCNFSM5AMNRCYQ>.

[cgrudz]

the trick is that this requires
having two distinct time-steppers for the ensemble and the the truth
twin respectively. Is there any example where this is performed
elsewhere in DAPPER?

I think this will do it (in the case of examples/basic_1.py):

@@ -15,6 +15,8 @@ import dapper.da_methods as da
 # #### Load experiment setup: the hidden Markov model (HMM)
 
 from dapper.mods.Lorenz63.sakov2012 import HMM
+HMM2 = HMM.copy()
+HMM2.step = my_custom_step_function
 
 # #### Generate the same random numbers each time this script is run
 
@@ -35,7 +37,7 @@ xp = da.EnKF('Sqrt', N=10, infl=1.02, rot=True)
 
 # #### Assimilate yy, knowing the HMM; xx is used to assess the performance
 
-xp.assimilate(HMM, xx, yy, liveplots=not nb)
+xp.assimilate(HMM2, xx, yy, liveplots=not nb)
 
 # #### Average the time series of various statistics

Ah ha, so it seems to be working so far, but I think I notice now a possible issue with the next step. Part of the configuration that I'm trying to build requires that the model twin is actually using a more coarse discretization in time, e.g.,

• truth twin h = 0.005
• ensemble twin h=0.01
• \Delta_obs = 0.1

In this case, these align on the same absolute times when observations are given, but these have different partitions of the time-interval for the dynamic propagation.

I'll look into this to see if there's a good way to implement this with minimal invasion to the code, this can probably be made to work within the example scrip itself if I just end up needing to subset the finer discretization.

[patnr]

Closed in favour of #77