Skip to content

Commit 2f6b044

Browse files
author
Abhinab Bhattacharjee
committed
Added 2d problems, standalone examples
1 parent 02a7fec commit 2f6b044

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+575
-223
lines changed

Diff for: ‎README.md

+1-1
Original file line numberDiff line numberDiff line change
@@ -2,4 +2,4 @@
22

33
A set of Data assimilation tools. Filter smoothers, ensemble and variational methods, localization, etc.
44

5-
An active work in progress.
5+
An active work in progress.

Diff for: ‎docs/CONTRIBUTING.md

+66-1
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,72 @@ To request a pull request, please follow the steps below
1111
* Create a draft pull request
1212
* Do not merge or modify the master branch
1313

14+
## Creating and running problems
15+
TODO
16+
1417
## Style Conventions
1518

16-
In order to prevent confusions and maintain a consistent coding style, the following styles should be strictly followed. The conventions match most of the MATLAB's own conventions. Failing to follow the conventions can delay in processing and approving new pull requests.
19+
In order to prevent confusions and maintain a consistent coding style, the following styles should be strictly followed.
20+
The conventions match most of the MATLAB's own conventions. Failing to follow the conventions can delay in processing
21+
and approving new pull requests.
22+
23+
### Line Formatting
24+
25+
Four spaces are used for indentation where required. A line should be kept less than 120 characters.
26+
27+
### Variables
28+
29+
Variable names should follow camelCasing convention.
30+
31+
```
32+
% Example
33+
N = 100;
34+
xValue = linspace(-5,5,N)
35+
zValue = @(x, sigma, n) exp(-x) + sig*randn(1,n)
36+
```
37+
38+
### Functions
39+
40+
Function names should be in camelCase. No special characters (for example '_') shall be used.
41+
42+
```
43+
% Example
44+
function val = myFun(args1, args2, ...)
45+
...
46+
end
47+
```
48+
49+
50+
### Structures
51+
52+
Structures should have camelCase property names
53+
54+
```
55+
% Example
56+
filter = struct('Type', 'Ensemble', 'Name', 'EnKF')
57+
```
58+
59+
60+
### Classes
61+
62+
Class names and properties should be in PascalCase. Acronyms should be capitalized. Methods will have camelCase as before.
63+
64+
```
65+
% Example
66+
classdef ExampleClass < handle
67+
properties
68+
PropertyOne
69+
PropertyTwo
70+
PropertyACRONYM
71+
end
72+
73+
methods
74+
function val = fun(...)
75+
...
76+
end
77+
end
78+
end
79+
80+
```
81+
1782

Diff for: ‎src/+datools/+examples/calc_B_l96.m renamed to ‎src/+datools/+examples/+oldscripts/calc_B_l96.m

File renamed without changes.

Diff for: ‎src/+datools/+examples/l63.m renamed to ‎src/+datools/+examples/+oldscripts/l63.m

File renamed without changes.

Diff for: ‎src/+datools/+examples/l63_RHF.m renamed to ‎src/+datools/+examples/+oldscripts/l63_RHF.m

+2-3
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
modelODE.TimeSpan = [0, Delta_t];
1919

2020
% Propogate
21-
[tt, yy] = ode45(natureODE.Rhs.F, [0, 10], nature0);
21+
[tt, yy] = ode45(natureODE.RHS.F, [0, 10], nature0);
2222
natureODE.Y0 = yy(end, :).';
2323

2424
model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
@@ -27,15 +27,14 @@
2727
naturetomodel = datools.observation.Linear(numel(nature0), 'H', ...
2828
speye(natureODE.NumVars));
2929

30-
N:50, inf:1.020
3130
%observeindicies = 1:natureODE.NumVars;
3231
observeindicies = 1:1:natureODE.NumVars;
3332

3433
nobsvars = numel(observeindicies);
3534

3635
R = (1 / 1) * speye(nobsvars);
3736

38-
obserrormodel = datools.error.Gaussian('CovarianceSqrt', sqrtm(R));
37+
obserrormodel = datools.uncertainty.Gaussian('CovarianceSqrt', sqrtm(R));
3938
%obserrormodel = datools.error.Tent;
4039
observation = datools.observation.Indexed(model.NumVars, ...
4140
'ErrorModel', obserrormodel, ...

Diff for: ‎src/+datools/+examples/l96_hybrid.m renamed to ‎src/+datools/+examples/+oldscripts/l96_hybrid.m

File renamed without changes.

Diff for: ‎src/+datools/+examples/qgexample.m renamed to ‎src/+datools/+examples/+oldscripts/qgexample.m

+24-18
Original file line numberDiff line numberDiff line change
@@ -2,19 +2,25 @@
22
close all;
33

44
% time steps
5-
Deltat = 0.0109;
5+
dt = 0.0109;
6+
filtername = 'EnKF';
7+
filterType = 'Ensemble';
68

79
% Time Stepping Methods (Use ode45 or write your own)
8-
solvermodel = @(f, t, y) ode45(f, t, y);
9-
solvernature = @(f, t, y) ode45(f, t, y);
10+
% solvermodel = @(f, t, y) ode45(f, t, y);
11+
% solvernature = @(f, t, y) ode45(f, t, y);
12+
% solvermodel = @(f, t, y) datools.utils.eDP54(f, t, y);
13+
% solvernature = @(f, t, y) datools.utils.eDP54(f, t, y);
14+
solvermodel = @(f, t, y) datools.utils.rk4ens(f, t, y, 1);
15+
solvernature = @(f, t, y) datools.utils.rk4ens(f, t, y, 1);
1016

1117
% Define ODE
12-
natureODE = otp.qg.presets.Canonical('Size', [63, 127]);
18+
natureODE = otp.quasigeostrophic.presets.Canonical('Size', [63, 127]);
1319
nature0 = natureODE.Y0;
14-
natureODE.TimeSpan = [0, Deltat];
20+
natureODE.TimeSpan = [0, dt];
1521

16-
modelODE = otp.qg.presets.Canonical('Size', [63, 127]);
17-
modelODE.TimeSpan = [0, Deltat];
22+
modelODE = otp.quasigeostrophic.presets.Canonical('Size', [63, 127]);
23+
modelODE.TimeSpan = [0, dt];
1824

1925
% initialize model
2026
model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
@@ -29,20 +35,20 @@
2935

3036
R = (1 / 1) * speye(nobsvars);
3137

32-
obserrormodel = datools.error.Gaussian('CovarianceSqrt', sqrtm(R));
38+
natureobserrormodel = datools.uncertainty.Gaussian('Covariance', R);
3339
%obserrormodel = datools.error.Tent;
3440
observation = datools.observation.Indexed(model.NumVars, ...
35-
'ErrorModel', obserrormodel, ...
41+
'Uncertainty', natureobserrormodel, ...
3642
'Indices', observeindicies);
3743

3844

3945
% We make the assumption that there is no model error
40-
modelerror = datools.error.Error;
46+
modelerror = datools.uncertainty.NoUncertainty;
4147

4248
%ensembleGenerator = @(N) randn(natureODE.NumVars, N);
4349

4450

45-
load('qgtrajectory.mat')
51+
load('+datools/+examples/+sandu/private/qgtrajectory.mat')
4652
y = y.';
4753
Nt = size(y, 2);
4854

@@ -123,9 +129,8 @@
123129
%localization = @(t, y, Hi, k) datools.tapering.cutoffCTilde(t, y, r, d, Hi, k);
124130

125131

126-
filter = datools.statistical.ensemble.EnKF(model, ...
127-
'Observation', observation, ...
128-
'NumEnsemble', ensN, ...
132+
filter = datools.filter.ensemble.(filtername)(model, ...
133+
'InitialEnsemble', ensembleGenerator(ensN)/10, ...
129134
'ModelError', modelerror, ...
130135
'EnsembleGenerator', ensembleGenerator, ...
131136
'Inflation', inflation, ...
@@ -136,8 +141,9 @@
136141

137142
%filterglobal = filter;
138143

139-
filter.setMean(natureODE.Y0);
140-
filter.scaleAnomalies(1/10);
144+
% filter.setMean(natureODE.Y0);
145+
% filter.scaleAnomalies(1/10);
146+
filter.MeanEstimate = natureODE.Y0;
141147

142148
% define steps and spinups
143149
%spinup = 5;
@@ -163,8 +169,8 @@
163169

164170

165171
% observe
166-
xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
167-
y = filter.Observation.observeWithError(model.TimeSpan(1), xt);
172+
xt = naturetomodel.observeWithoutError(nature.ODEModel.TimeSpan(1), nature.State);
173+
y = filter.Observation.observeWithError(model.ODEMOdel.TimeSpan(1), xt);
168174

169175
% Rank histogram (if needed)
170176
datools.utils.stat.RH(filter, xt);

Diff for: ‎src/+datools/+examples/+sandu/l96_experiments.m

+33-28
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
clear all; close all;
1+
clear all; close all; clc;
22

33
%% User Inputs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44
% Use this to run Lorenz96 experiments
@@ -58,26 +58,30 @@
5858
Dt = 0.05;
5959

6060
% Time Stepping Methods (Use ode45 or write your own)
61-
solvermodel = @(f, t, y) datools.utils.rk4(f, t, y, 1);
62-
solvernature = @(f, t, y) datools.utils.rk4(f, t, y, 1);
61+
% solvermodel = @(f, t, y) datools.utils.rk4(f, t, y, 1);
62+
% solvernature = @(f, t, y) datools.utils.rk4(f, t, y, 1);
63+
solvermodel = @(f, t, y) datools.utils.rk4ens(f, t, y, 1);
64+
solvernature = @(f, t, y) datools.utils.rk4ens(f, t, y, 1);
6365

64-
% Define ODE
66+
% Define ODE for the truth
6567
natureODE = otp.lorenz96.presets.Canonical;
6668
nature0 = randn(natureODE.NumVars, 1);
6769
natureODE.TimeSpan = [0, Dt];
6870

71+
% define ODE for the model
6972
modelODE = otp.lorenz96.presets.Canonical;
7073
modelODE.TimeSpan = [0, Dt];
7174

72-
% Propogate the model
75+
% Propogate the truth
7376
[tt, yy] = ode45(natureODE.RHS.F, [0, 10], nature0);
7477
natureODE.Y0 = yy(end, :).';
7578

76-
% initialize model
79+
% initialize model object
7780
model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
7881
nature = datools.Model('Solver', solvernature, 'ODEModel', natureODE);
7982

80-
% Observation Model
83+
% Observation object that maps from nature to model
84+
% define H. Need to have an option using linearized H (user must supply)
8185
naturetomodel = datools.observation.Linear(numel(nature0), 'H', ...
8286
speye(natureODE.NumVars));
8387

@@ -88,15 +92,14 @@
8892

8993
R = variance * speye(nobsvars);
9094

91-
% Observaton model (Gaussian here)
92-
obserrormodel = datools.error.Gaussian('CovarianceSqrt', sqrtm(R));
93-
95+
% Observaton model object (Gaussian here)
96+
obserrormodel = datools.uncertainty.Gaussian('Covariance', R);
9497
observation = datools.observation.Indexed(model.NumVars, ...
95-
'ErrorModel', obserrormodel, ...
98+
'Uncertainty', obserrormodel, ...
9699
'Indices', observeindicies);
97100

98101
% We make the assumption that there is no model error
99-
modelerror = datools.error.Error;
102+
modelerror = datools.uncertainty.NoUncertainty;
100103

101104
% This can be used to generate ensemble if needed
102105
ensembleGenerator = @(N) randn(natureODE.NumVars, N);
@@ -179,17 +182,15 @@
179182

180183
switch filtername
181184
case 'EnKF'
182-
filter = datools.statistical.ensemble.(filtername)(model, ...
183-
'Observation', observation, ...
184-
'NumEnsemble', ensN, ...
185-
'ModelError', modelerror, ...
185+
filter = datools.filter.ensemble.(filtername)(model, ...
186+
'InitialEnsemble', ensembleGenerator(ensN)/10, ...
186187
'EnsembleGenerator', ensembleGenerator, ...
187188
'Inflation', inflation, ...
188189
'Parallel', false, ...
189190
'RankHistogram', histvar, ...
190191
'Rejuvenation', rejuvenation);
191192
case 'ETKF'
192-
filter = datools.statistical.ensemble.(filtername)(model, ...
193+
filter = datools.filter.ensemble.(filtername)(model, ...
193194
'Observation', observation, ...
194195
'NumEnsemble', ensN, ...
195196
'ModelError', modelerror, ...
@@ -199,7 +200,7 @@
199200
'RankHistogram', histvar, ...
200201
'Rejuvenation', rejuvenation);
201202
case 'LETKF'
202-
filter = datools.statistical.ensemble.(filtername)(model, ...
203+
filter = datools.filter.ensemble.(filtername)(model, ...
203204
'Observation', observation, ...
204205
'NumEnsemble', ensN, ...
205206
'ModelError', modelerror, ...
@@ -209,7 +210,7 @@
209210
'RankHistogram', histvar, ...
210211
'Rejuvenation', rejuvenation);
211212
case 'SIR'
212-
filter = datools.statistical.ensemble.(filtername)(model, ...
213+
filter = datools.filter.ensemble.(filtername)(model, ...
213214
'Observation', observation, ...
214215
'NumEnsemble', ensN, ...
215216
'ModelError', modelerror, ...
@@ -219,7 +220,7 @@
219220
'RankHistogram', histvar, ...
220221
'Rejuvenation', rejuvenation);
221222
case 'ETPF'
222-
filter = datools.statistical.ensemble.(filtername)(model, ...
223+
filter = datools.filter.ensemble.(filtername)(model, ...
223224
'Observation', observation, ...
224225
'NumEnsemble', ensN, ...
225226
'ModelError', modelerror, ...
@@ -229,7 +230,7 @@
229230
'RankHistogram', histvar, ...
230231
'Rejuvenation', rejuvenation);
231232
case 'RHF'
232-
filter = datools.statistical.ensemble.(filtername)(model, ...
233+
filter = datools.filter.ensemble.(filtername)(model, ...
233234
'Observation', observation, ...
234235
'NumEnsemble', ensN, ...
235236
'ModelError', modelerror, ...
@@ -241,9 +242,10 @@
241242

242243
end
243244

244-
filter.setMean(natureODE.Y0);
245-
filter.scaleAnomalies(1/10);
246-
245+
% filter.setMean(natureODE.Y0);
246+
% filter.scaleAnomalies(1/10);
247+
248+
filter.MeanEstimate = natureODE.Y0;
247249
mses = zeros(steps-spinup, 1);
248250

249251
rmse = nan;
@@ -262,22 +264,25 @@
262264
end
263265

264266
% observe
265-
xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
266-
y = filter.Observation.observeWithError(model.TimeSpan(1), xt);
267+
xt = naturetomodel.observeWithoutError(nature.ODEModel.TimeSpan(1), nature.State);
268+
% y = filter.Observation.observeWithError(model.TimeSpan(1), xt);
269+
y = observation.observeWithError(xt);
270+
observation.Uncertainty.Mean = y;
271+
267272

268273
% Rank histogram (if needed)
269274
datools.utils.stat.RH(filter, xt);
270275

271276
% analysis
272277
% try
273278
if dofilter
274-
filter.analysis(R, y);
279+
filter.analysis(observation);
275280
end
276281
%catch
277282
% dofilter = false;
278283
%end
279284

280-
xa = filter.BestEstimate;
285+
xa = filter.MeanEstimate;
281286

282287
err = xt - xa;
283288

0 commit comments

Comments
 (0)