Revision

[patnr]

@cgrudz After merging the PR (#77), I also went through and re-did the example scripts. I hope you will find them improved in their simplicity. Please go through and see if you understand the new scripts, ensure that they make sense, and that they produce sensible results. How to checkout PRs
Update: relevant commit

Please also fill in the TODOs (use Ctrl-F here to find them).

[patnr]

Btw, AFAICT, there's no need to keep the deterministic and stochastic rk4 separate. I was thus able to simplify the rk4 code. Can you have a look?

472bcf8

[cgrudz]

Word, sorry I'm dragging my feet on this -- I have all my time going to submitting my new manuscript ASAP.  I'll try to get back to development in a week or two.  Also, once I get my manuscript up, we'll need to chat about integrating the SIEnKS into DAPPER. Cheers, Colin

…

On 8/31/21 8:14 AM, Patrick N. Raanes wrote: Btw, AFAICT, there's no need to keep the deterministic and stochastic rk4 separate. I was thus able to simplify the |rk4| code. Can you have a look? 472bcf8 <472bcf8> — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#79 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABXIPFXIPEFDEBFWASZ646DT7TWVPANCNFSM5CTCMSZQ>. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[patnr]

No worries. That's the way it goes. But I do need you to double check the validity of these changes. In the process, you will hopefully pick up some more of DAPPER.

[patnr]

Edit of the above: But I do need you to double check the validity of these changes when you find the time

[cgrudz]

Good news and bad news.

Good news, I submitted my manuscript and I had a chance to look over this today.

Bad news, the results that I reproduced with my earlier scripts are now broken in the examples provided, I think some of the changes are OK but some things appear to have gotten lost in this. I'm also not entirely sure based on the changes where this is occurring. Let me explain what the old script was reproducing, so that maybe you have a better sense of what this was intended to study.

In particular, we were looking at the following configurations:

High-precision (unrealistic for most twin experiments due to the extremely fine discretization)

• Ensemble 1: RK4 uses h=0.001
• Ensemble 2: EM uses h=0.001
• Truth twin: Taylor uses h=0.001
• observation error std and diffusion range {0.1, 0.25, 0.5, 0.75, 1.0}

The ensemble filters generated with the RK4 scheme and EM scheme perform roughly as good as each other, without the presence of filter divergence. This shows that it is possible to use the EM scheme with the extremely fine discretization step.

Low-precision (but statistically robust when comparing with the high-precision setting and using RK)

• Ensemble 1: RK4 uses h=0.01
• Ensemble 2: EM uses h=0.01
• Truth twin: Taylor uses h=0.005
• observation error std and diffusion range {0.1, 0.25, 0.5, 0.75, 1.0}

The ensemble filters generated by the RK4 scheme perform about as well as in the high-precision setting, but EM experiences filter divergence for low noise settings. On the other hand, in the high noise settings, RK and EM perform about as well as each other, if I recall, having an average RMSE around the 0.5 mark. This demonstrates that the EM scheme is not statistically robust as the RK scheme, but that when model uncertainty dominates the dynamics, we can still get away with a very low precision scheme like EM.

This second scenario is roughly visualized with the following plot online (https://gmd.copernicus.org/articles/13/1903/2020/gmd-13-1903-2020-f08-web.png) where it is plotting the difference between the low-precision EM scheme and a higher precision ensemble forecast.

[cgrudz]

In commit f9fbb0b, I produced the following with my script for the Euler-Maruyama scheme with h=0.01 and Taylor 0.005

    Diffusion  ObsNoise  |  rmse.a  1σ     rmv.a  1σ     rmse.f  1σ  
--  ---------  --------  -  -------------  ------------  ------------
 0       0.1       0.1   |   2     ±1      0.28  ±0.01     2    ±1   
 1       0.1       0.25  |   1.3   ±0.5    0.301 ±0.007    1.5  ±0.5 
 2       0.1       0.5   |   1.6   ±0.1    0.312 ±0.006    1.8  ±0.1 
 3       0.1       0.75  |   1.1   ±0.2    0.30  ±0.01     1.3  ±0.2 
 4       0.1       1     |   2.6   ±0.7    0.322 ±0.008    2.8  ±0.7 
 5       0.25      0.1   |   0.41  ±0.04   0.340 ±0.004    0.51 ±0.05
 6       0.25      0.25  |   0.55  ±0.05   0.349 ±0.002    0.67 ±0.06
 7       0.25      0.5   |   0.6   ±0.1    0.341 ±0.008    0.7  ±0.1 
 8       0.25      0.75  |   0.6   ±0.1    0.350 ±0.006    0.8  ±0.1 
 9       0.25      1     |   0.56  ±0.04   0.348 ±0.007    0.69 ±0.06
10       0.5       0.1   |   0.32  ±0.01   0.412 ±0.006    0.43 ±0.02
11       0.5       0.25  |   0.34  ±0.01   0.413 ±0.004    0.45 ±0.01
12       0.5       0.5   |   0.36  ±0.01   0.405 ±0.006    0.46 ±0.02
13       0.5       0.75  |   0.47  ±0.02   0.427 ±0.003    0.62 ±0.03
14       0.5       1     |   0.46  ±0.02   0.423 ±0.002    0.60 ±0.03
15       0.75      0.1   |   0.33  ±0.01   0.473 ±0.003    0.44 ±0.02
16       0.75      0.25  |   0.374 ±0.009  0.470 ±0.002    0.50 ±0.02
17       0.75      0.5   |   0.41  ±0.01   0.475 ±0.004    0.55 ±0.02
18       0.75      0.75  |   0.42  ±0.01   0.468 ±0.003    0.54 ±0.02
19       0.75      1     |   0.49  ±0.02   0.473 ±0.003    0.62 ±0.03
20       1         0.1   |   0.36  ±0.01   0.521 ±0.002    0.53 ±0.02
21       1         0.25  |   0.38  ±0.01   0.513 ±0.002    0.53 ±0.02
22       1         0.5   |   0.45  ±0.01   0.522 ±0.002    0.63 ±0.02
23       1         0.75  |   0.48  ±0.01   0.518 ±0.003    0.65 ±0.02
24       1         1     |   0.53  ±0.02   0.515 ±0.004    0.69 ±0.03

This is closer to what I remember seeing in my own own results, though it's possible there was still a bug here. However, the current results don't look much like this at all.

[cgrudz]

The RK scheme produced the following under the same conditions above

    Diffusion  ObsNoise  |  rmse.a  1σ     rmv.a  1σ     rmse.f  1σ   
--  ---------  --------  -  -------------  ------------  -------------
 0       0.1       0.1   |   0.140 ±0.006  0.289 ±0.006   0.175 ±0.008
 1       0.1       0.25  |   0.177 ±0.006  0.26  ±0.01    0.21  ±0.01 
 2       0.1       0.5   |   0.18  ±0.01   0.281 ±0.007   0.21  ±0.02 
 3       0.1       0.75  |   0.22  ±0.01   0.269 ±0.009   0.26  ±0.02 
 4       0.1       1     |   0.23  ±0.02   0.243 ±0.005   0.26  ±0.02 
 5       0.25      0.1   |   0.17  ±0.01   0.31  ±0.02    0.22  ±0.02 
 6       0.25      0.25  |   0.227 ±0.008  0.326 ±0.006   0.28  ±0.01 
 7       0.25      0.5   |   0.31  ±0.02   0.328 ±0.007   0.37  ±0.02 
 8       0.25      0.75  |   0.27  ±0.02   0.325 ±0.006   0.34  ±0.02 
 9       0.25      1     |   0.31  ±0.01   0.336 ±0.002   0.37  ±0.01 
10       0.5       0.1   |   0.25  ±0.01   0.401 ±0.005   0.33  ±0.02 
11       0.5       0.25  |   0.29  ±0.01   0.402 ±0.003   0.39  ±0.02 
12       0.5       0.5   |   0.33  ±0.01   0.406 ±0.003   0.41  ±0.02 
13       0.5       0.75  |   0.40  ±0.02   0.399 ±0.004   0.49  ±0.02 
14       0.5       1     |   0.39  ±0.02   0.399 ±0.004   0.49  ±0.03 
15       0.75      0.1   |   0.305 ±0.008  0.457 ±0.002   0.41  ±0.01 
16       0.75      0.25  |   0.33  ±0.01   0.459 ±0.003   0.45  ±0.01 
17       0.75      0.5   |   0.392 ±0.009  0.459 ±0.002   0.53  ±0.02 
18       0.75      0.75  |   0.40  ±0.01   0.453 ±0.003   0.52  ±0.02 
19       0.75      1     |   0.45  ±0.02   0.459 ±0.003   0.58  ±0.02 
20       1         0.1   |   0.344 ±0.008  0.506 ±0.003   0.49  ±0.01 
21       1         0.25  |   0.36  ±0.01   0.509 ±0.003   0.51  ±0.02 
22       1         0.5   |   0.42  ±0.01   0.506 ±0.002   0.57  ±0.02 
23       1         0.75  |   0.47  ±0.01   0.507 ±0.002   0.64  ±0.02 
24       1         1     |   0.49  ±0.01   0.515 ±0.002   0.68  ±0.02 

[cgrudz]

Under fine integration conditions described above, the RK scheme produced

    Diffusion  ObsNoise  |  rmse.a  1σ     rmv.a  1σ     rmse.f  1σ  
--  ---------  --------  -  -------------  ------------  ------------
 0       0.1       0.1   |   0.152 ±0.008  0.275 ±0.006    0.19 ±0.01
 1       0.1       0.25  |   0.17  ±0.01   0.23  ±0.01     0.19 ±0.02
 2       0.1       0.5   |   0.192 ±0.009  0.27  ±0.01     0.23 ±0.01
 3       0.1       0.75  |   0.22  ±0.01   0.276 ±0.008    0.26 ±0.01
 4       0.1       1     |   0.24  ±0.02   0.273 ±0.005    0.29 ±0.02
 5       0.25      0.1   |   0.21  ±0.01   0.340 ±0.003    0.26 ±0.01
 6       0.25      0.25  |   0.23  ±0.01   0.30  ±0.01     0.28 ±0.01
 7       0.25      0.5   |   0.26  ±0.01   0.328 ±0.005    0.32 ±0.02
 8       0.25      0.75  |   0.30  ±0.01   0.337 ±0.003    0.38 ±0.02
 9       0.25      1     |   0.31  ±0.01   0.326 ±0.006    0.37 ±0.02
10       0.5       0.1   |   0.26  ±0.01   0.396 ±0.003    0.34 ±0.02
11       0.5       0.25  |   0.287 ±0.007  0.402 ±0.004    0.37 ±0.01
12       0.5       0.5   |   0.33  ±0.01   0.411 ±0.003    0.42 ±0.02
13       0.5       0.75  |   0.38  ±0.01   0.405 ±0.003    0.48 ±0.02
14       0.5       1     |   0.37  ±0.02   0.404 ±0.003    0.46 ±0.02
15       0.75      0.1   |   0.32  ±0.01   0.454 ±0.002    0.45 ±0.01
16       0.75      0.25  |   0.313 ±0.006  0.461 ±0.003    0.42 ±0.02
17       0.75      0.5   |   0.40  ±0.01   0.461 ±0.003    0.53 ±0.02
18       0.75      0.75  |   0.42  ±0.01   0.459 ±0.003    0.55 ±0.02
19       0.75      1     |   0.44  ±0.01   0.457 ±0.005    0.56 ±0.02
20       1         0.1   |   0.36  ±0.01   0.506 ±0.003    0.52 ±0.02
21       1         0.25  |   0.370 ±0.008  0.500 ±0.003    0.50 ±0.01
22       1         0.5   |   0.43  ±0.01   0.508 ±0.003    0.57 ±0.03
23       1         0.75  |   0.48  ±0.01   0.503 ±0.003    0.62 ±0.02
24       1         1     |   0.50  ±0.01   0.507 ±0.002    0.65 ±0.01

[cgrudz]

Under the fine conditions described above, the EM scheme produced

    Diffusion  ObsNoise  |  rmse.a  1σ     rmv.a  1σ     rmse.f  1σ  
--  ---------  --------  -  -------------  ------------  ------------
 0       0.1       0.1   |   0.14  ±0.01   0.275 ±0.006    0.17 ±0.01
 1       0.1       0.25  |   0.16  ±0.01   0.275 ±0.009    0.20 ±0.02
 2       0.1       0.5   |   0.20  ±0.02   0.274 ±0.009    0.24 ±0.02
 3       0.1       0.75  |   0.25  ±0.01   0.29  ±0.01     0.31 ±0.02
 4       0.1       1     |   0.31  ±0.03   0.270 ±0.005    0.35 ±0.04
 5       0.25      0.1   |   0.20  ±0.02   0.324 ±0.004    0.24 ±0.02
 6       0.25      0.25  |   0.23  ±0.02   0.331 ±0.004    0.29 ±0.03
 7       0.25      0.5   |   0.253 ±0.008  0.332 ±0.007    0.32 ±0.01
 8       0.25      0.75  |   0.30  ±0.02   0.330 ±0.005    0.36 ±0.03
 9       0.25      1     |   0.30  ±0.01   0.325 ±0.006    0.35 ±0.02
10       0.5       0.1   |   0.244 ±0.007  0.400 ±0.005    0.32 ±0.01
11       0.5       0.25  |   0.29  ±0.02   0.403 ±0.003    0.39 ±0.02
12       0.5       0.5   |   0.34  ±0.01   0.395 ±0.004    0.41 ±0.02
13       0.5       0.75  |   0.36  ±0.01   0.397 ±0.003    0.44 ±0.01
14       0.5       1     |   0.40  ±0.01   0.401 ±0.004    0.49 ±0.02
15       0.75      0.1   |   0.313 ±0.007  0.459 ±0.003    0.43 ±0.01
16       0.75      0.25  |   0.33  ±0.01   0.461 ±0.002    0.43 ±0.01
17       0.75      0.5   |   0.37  ±0.01   0.465 ±0.002    0.49 ±0.02
18       0.75      0.75  |   0.42  ±0.01   0.462 ±0.002    0.54 ±0.02
19       0.75      1     |   0.45  ±0.02   0.460 ±0.006    0.58 ±0.04
20       1         0.1   |   0.338 ±0.008  0.510 ±0.003    0.49 ±0.01
21       1         0.25  |   0.368 ±0.009  0.509 ±0.002    0.52 ±0.01
22       1         0.5   |   0.45  ±0.02   0.509 ±0.003    0.62 ±0.02
23       1         0.75  |   0.45  ±0.02   0.510 ±0.003    0.59 ±0.02
24       1         1     |   0.54  ±0.02   0.512 ±0.002    0.71 ±0.03

[cgrudz]

Digging out my old simulation data from the manuscript, this is what I find for the coarse simulations (note the difference where I did the silly mistake of using the variance rather than std for the obs error):

EM
    Diffusion  ObsNoise  |  rmse.a 
--  ---------  --------  -  -------------  ------------  ------------
 0       0.1       0.1   |   0.4165107 
 1       0.1       0.25  |   0.93195499 
 2       0.1       0.5   |   1.91500919
 3       0.1       0.75  |  2.20296762
 4       0.1       1     |   2.30187736
 5       0.25      0.1   |  0.23268683
 6       0.25      0.25  |  0.33309646
 7       0.25      0.5   |  0.43717187
 8       0.25      0.75  |  0.51492269
 9       0.25      1     |   0.58213909
10       0.5       0.1   |  0.2174468 
11       0.5       0.25  | 0.30138099
12       0.5       0.5   | 0.38200366
13       0.5       0.75  | 0.44159024
14       0.5       1     |   0.48785175
15       0.75      0.1   |  0.23283048
16       0.75      0.25  | 0.31937611
17       0.75      0.5   |  0.40261691
18       0.75      0.75  | 0.4620334
19       0.75      1     |    0.50735597
20       1         0.1   |   0.24832587
21       1         0.25  | 0.34300861
22       1         0.5   |   0.43191628
23       1         0.75  | 0.49506043
24       1         1     |    0.5420419


RK
    Diffusion  ObsNoise  |  rmse.a 
--  ---------  --------  -  -------------  ------------  ------------
 0       0.1       0.1   |   0.11107565
 1       0.1       0.25  |   0.15544081 
 2       0.1       0.5   |   0.20184683, 
 3       0.1       0.75  |  0.23704249
 4       0.1       1     |   0.2701058
 5       0.25      0.1   |  0.1513685 
 6       0.25      0.25  |   0.2045175 
 7       0.25      0.5   |  0.2596008 
 8       0.25      0.75  |  0.29942023
 9       0.25      1     |   0.33385369
10       0.5       0.1   |  0.19469195 
11       0.5       0.25  | 0.26022757
12       0.5       0.5   | 0.32560759
13       0.5       0.75  | 0.37123433
14       0.5       1     |   0.40830518
15       0.75      0.1   |  0.22376344
16       0.75      0.25  | 0.30305462
17       0.75      0.5   |  0.3773218 
18       0.75      0.75  | 0.42866478
19       0.75      1     |     0.46934533
20       1         0.1   |   0.24352176
21       1         0.25  | 0.33453447
22       1         0.5   |   0.41865509
23       1         0.75  | 0.47717651
24       1         1     |    0.5243447

Again this is very similar to the previous commit, but things seem to be breaking now.

[patnr]

Hmm, ok, I will have a new look. I remember checking that the first script did not change the results. But for the second script I did not check.

[patnr]

Let me explain what the old script was reproducing, ...

IMO I had indeed understood this when making the changes. Is there somewhere in the current code that does not seem to conform to this purpose?

(note the difference where I did the silly mistake of using the variance rather than std for the obs error):

I don't see what you're referring to here

[patnr]

Hmmm.

Actually #79 contained one bugfix. The setting Obs['noise'] = xp.ObsNoise from here (merge commit of #77) only affects the truth simulation, because the original hmm (used by the DA methods) is unchanged.

And now I found another bug: we forgot that the rename of "Diffusion" to "diffusion" must also be done in the script! 🤕
Fixed in c6c1d60

I then computed results from #77, #79 and c6c1d60, and compared to your "old simulation data"
For easier comparison I increased KObs=1000 and only used 3 values of D and R.
Here's what I got (tables edited/merged by hand)

All tables are for **low-res (coarse)** simulations,
since that's what I have availble from your "old simulation data"

              #77      (g)old       #79    c6c1d60
EM
D    R    |  rmse.a  |  rmse.a |  rmse.a  | rmse.a
--------------------------------------------------
0.1  0.1  |     1.7  |  0.416  |   1.9    |  0.41
0.1  0.5  |     2.3  |  1.915  |   0.42   |  1.6
0.1  1    |     2.0  |  2.301  |   2.5    |  2.4
0.5  0.1  |     2.2  |  0.217  |   0.42   |  0.215
0.5  0.5  |     1.7  |  0.382  |   1.9    |  0.374
0.5  1    |     2.4  |  0.487  |   2.5    |  0.47
1    0.1  |     1.8  |  0.248  |   0.42   |  0.246
1    0.5  |     2.2  |  0.431  |   1.9    |  0.425
1    1    |     1.9  |  0.542  |   2.5    |  0.53

RK4
D    R    |  rmse.a  |  rmse.a |  rmse.a  |  rmse.a
---------------------------------------------------
0.1  0.1  |   0.145  |  0.111  |   0.106  |  0.105
0.1  0.5  |   0.198  |  0.201  |   0.193  |  0.188
0.1  1    |   0.256  |  0.270  |   0.26   |  0.25
0.5  0.1  |   0.139  |  0.194  |   0.106  |  0.192
0.5  0.5  |   0.203  |  0.325  |   0.193  |  0.317
0.5  1    |   0.270  |  0.408  |   0.26   |  0.394
1    0.1  |   0.141  |  0.243  |   0.106  |  0.242
1    0.5  |   0.209  |  0.418  |   0.193  |  0.411
1    1    |   0.259  |  0.524  |   0.26   |  0.510

The correspondence is now remarkably good! I think we can be pretty happy with this. A very close reproduction of your original results!

[cgrudz]

I don't see what you're referring to here

I had two axes in the figure that, to remain on the same scale, should have increased in both directions with respect to the standard deviation, rather than standard deviation / variance... Silly on my part, but remarkably silly that the referees and co-authors didn't catch this either... XD Not that this makes the results erroneous, just harder to compare based on the linear and quadratic scales side-by-side.

The numerics above do look really good -- as long as the fine simulation isn't so far different in either case (EM / RK) from the coarse RK simulation, then this validates that we can produce a statistically robust twin experiment in the suggested configuration in a way that doesn't bias the outcome away from the extremely fine ground truth. The EM scheme definitely biases this, as evidenced in a variety of places in the original manuscript, but the RK simulation makes for a very good out-of-the-box SDE ensemble integration scheme.

I'll make a pull of the current versions and look back over. Also, I'll start digging into this a bit to inspect how we might merge in the SIEnKS method. You are a recommended referee, in part because you're one of a handful of people who can review this without needing to read the long mathematical intro. I think if you aren't a referee, you can basically just jump to the punch line with the numerical benchmarks. In any case, I think you'll like the numerics, but if I rope you into reading the whole manuscript, you'll probably have some good critical feedback... hehehehe

[patnr]

I see. Well I didn't actually look it up. I have been more concerned with the actual reproduction than with the interpretation. I think I've spent several days total on these two PRs, so hopefully next time will be easier. I leave it to you to verify the fine / high-res case. Congrats on submitting! I don't know if I have time right now for a review. In any case, if I do, I will sign it, just FYI.

[cgrudz]

Hahahaha, indeed... This has been a fair amount of work on both of our sides I think, I hope that the next roll will come more quickly.  In any case, I think there's some fun additional pieces of research to put in this that will be good for all respective projects involved.

…

On 9/5/21 12:05 PM, Patrick N. Raanes wrote: I see. Well I didn't actually look it up. I have been more concerned with the actual reproduction than with the interpretation. I think I've spent several days total on these two PRs, so hopefully next time will be easier. I leave it to you to verify the fine / high-res case. Congrats on submitting! I don't know if I have time right now for a review. In any case, if I do, I will sign it, just FYI. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#79 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABXIPFQN5LZNDONUCJONVJDUAO5OBANCNFSM5CTCMSZQ>. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.