Skip to content

Commit 22dc89e

Browse files
author
Amit N Subrahmanya
committed
Added weights, SIR
1 parent b12cd8a commit 22dc89e

File tree

11 files changed

+125
-347
lines changed

11 files changed

+125
-347
lines changed

‎.Rhistory

Whitespace-only changes.

‎src/+datools/+examples/l96_enkf.m

+5-4
Original file line numberDiff line numberDiff line change
@@ -43,8 +43,9 @@
4343

4444
ensembleGenerator = @(x) randn(natureODE.NumVars, x);
4545

46-
ensNs = 5:5:50;
46+
ensNs = 100:100:500;
4747
infs = 1.01:.01:1.05;
48+
infs = 0.1:0.1:0.7;
4849

4950
rmses = inf*ones(numel(ensNs), numel(infs));
5051

@@ -80,16 +81,16 @@
8081
inflation = inflationAll;
8182

8283
% No localization
83-
% r = 5;
84+
% r = 3;
8485
% d = @(t, y, i, j) modelODE.DistanceFunction(t, y, i, j);
8586
localization = [];
8687

8788
%localization= @(t, y, H) datools.tapering.gc(t, y, r, d, H);
8889

89-
%localization = @(t, y, Hi, k) datools.tapering.gcCTilde(t, y, r, d, Hi, k);
90+
% localization = @(t, y, Hi, k) datools.tapering.gcCTilde(t, y, Hi, r, d, k);
9091
%localization = @(t, y, Hi, k) datools.tapering.cutoffCTilde(t, y, r, d, Hi, k);
9192

92-
enkf = datools.statistical.ensemble.ETPF(model, ...
93+
enkf = datools.statistical.ensemble.SIR(model, ...
9394
'Observation', observation, ...
9495
'NumEnsemble', ensN, ...
9596
'ModelError', modelerror, ...

‎src/+datools/+examples/new.m

+60-136
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,18 @@
11
clear;
22
close all;
3-
figure;
4-
drawnow;
53

6-
Deltat = 0.05;
4+
rng(5);
75

6+
Deltat = 0.12;
87

9-
solvermodel = @(f, t, y) datools.utils.rk4(f, t, y, 1);
10-
11-
solvernature = @(f, t, y) datools.utils.rk4(f, t, y, 1);
12-
8+
solvermodel = @(f, t, y) datools.utils.rk4(f, t, y, 100);
9+
solvernature = @(f, t, y) datools.utils.rk4(f, t, y, 100);
1310

11+
% solvermodel = @(f, t, y) ode45(f, t, y);
12+
% solvernature = @(f, t, y) ode45(f, t, y);
1413

1514
natureODE = otp.lorenz63.presets.Canonical;
15+
natureODE.TimeSpan = [0, Deltat];
1616

1717
nvrs = natureODE.NumVars;
1818

@@ -22,173 +22,97 @@
2222
natureODE.Y0 = yy(end, :).';
2323

2424
modelODE = otp.lorenz63.presets.Canonical;
25+
modelODE.TimeSpan = [0, Deltat];
2526

2627

2728
model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
2829
nature = datools.Model('Solver', solvernature, 'ODEModel', natureODE);
2930

3031
naturetomodel = datools.observation.Linear(numel(nature0), 'H', speye(nvrs));
3132

32-
33-
observeindicies = 1:nvrs;
33+
observeindicies = 1;
3434

3535
nobsvars = numel(observeindicies);
3636

37-
R = (1/1)*speye(nobsvars);
37+
R = (8/1)*speye(nobsvars);
3838

3939
obserrormodel = datools.error.Gaussian('CovarianceSqrt', sqrtm(R));
4040
observation = datools.observation.Indexed(model.NumVars, ...
4141
'ErrorModel', obserrormodel, ...
4242
'Indicies', observeindicies);
4343

44-
%% Do the rest
45-
4644
% We make the assumption that there is no model error
4745
modelerror = datools.error.Error;
4846

49-
50-
5147
ensembleGenerator = @(x) randn(nvrs, x);
5248

53-
ensNs = 5:5:50;
54-
infs = 1.05:.05:1.4;
49+
ensN = 50;
50+
infl = 1.05;
51+
rej = 0.1;
5552

56-
rmses = inf*ones(numel(ensNs), numel(infs));
53+
% No localization
54+
localization = [];
5755

56+
meth = datools.statistical.ensemble.ETPF2(model, ...
57+
'ModelError', modelerror, ...
58+
'Observation', observation, ...
59+
'NumEnsemble', ensN, ...
60+
'EnsembleGenerator', ensembleGenerator, ...
61+
'Inflation', infl, ...
62+
'Rejuvenation', rej, ...
63+
'Localization', localization, ...
64+
'Parallel', false, ...
65+
'RIPIterations', 0);
5866

59-
maxallowerr = 2;
60-
61-
mm = min(rmses(:));
62-
if mm >= maxallowerr
63-
mm = 0;
64-
end
67+
spinup = 200;
68+
times = 11*spinup;
6569

66-
imagesc(ensNs, infs, rmses.'); caxis([mm, 1]); colorbar; set(gca,'YDir','normal');
67-
axis square; title('EnKF'); colormap('pink');
68-
xlabel('Ensemble Size'); ylabel('Inflation');
70+
do_enkf = true;
6971

70-
runsleft = find(rmses == inf);
72+
sse = 0;
7173

72-
for runn = runsleft.'
73-
[ensNi, infi] = ind2sub([numel(ensNs), numel(infs)], runn);
74+
for i = 1:times
75+
%fprintf('%d|', i);
76+
% forecast
7477

75-
fprintf('N: %d, inf: %.3f\n', ensNs(ensNi), infs(infi));
78+
nature.evolve();
7679

77-
ns = 1;
78-
sE = zeros(ns, 1);
80+
if do_enkf
81+
meth.forecast();
82+
end
7983

80-
inflationAll = infs(infi);
81-
ensN = ensNs(ensNi);
8284

83-
for sample = 1:ns
84-
% Set rng for standard experiments
85-
rng(17 + sample - 1);
86-
87-
inflation = inflationAll;
88-
89-
% No localization
90-
localization= [];
91-
92-
%fprintf('1\n');
93-
94-
enkf = datools.statistical.ensemble.ETKF(model, ...
95-
'Observation', observation, ...
96-
'NumEnsemble', ensN, ...
97-
'ModelError', modelerror, ...
98-
'EnsembleGenerator', ensembleGenerator, ...
99-
'Inflation', inflation, ...
100-
'Localization', localization, ...
101-
'Parallel', false, ...
102-
'RIPIterations', 0);
103-
104-
spinup = 200;
105-
times = 11*spinup;
106-
107-
mses = zeros(times - spinup, 1);
108-
109-
rmse = nan;
110-
111-
ps = '';
112-
113-
do_enkf = true;
114-
115-
for i = 1:times
116-
%fprintf('%d|', i);
117-
% forecast
118-
119-
nature.evolve();
120-
121-
if do_enkf
122-
enkf.forecast();
123-
end
124-
125-
126-
% observe
127-
xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
128-
y = enkf.Observation.observeWithError(model.TimeSpan(1), xt);
129-
130-
% analysis
131-
132-
try
133-
if do_enkf
134-
enkf.analysis(R, y);
135-
end
136-
catch
137-
do_enkf = false;
138-
end
139-
140-
xa = enkf.BestEstimate;
141-
142-
err = xt - xa;
143-
144-
if i > spinup
145-
146-
mses(i - spinup) = mean((xa - xt).^2);
147-
rmse = sqrt(mean(mses(1:(i - spinup))));
148-
149-
if rmse > maxallowerr || isnan(rmse) || mses(i - spinup) > 2*maxallowerr
150-
do_enkf = false;
151-
end
152-
end
153-
154-
155-
if ~do_enkf
156-
break;
157-
end
158-
159-
end
160-
161-
if isnan(rmse)
162-
rmse = 1000;
163-
end
164-
165-
if ~do_enkf
166-
sE(sample) = 1000;
167-
else
168-
sE(sample) = rmse;
169-
end
170-
171-
%fprintf('\n');
172-
end
85+
% observe
86+
xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
87+
y = meth.Observation.observeWithError(model.TimeSpan(1), xt);
17388

174-
resE = mean(sE);
89+
% analysis
17590

176-
if isnan(resE)
177-
resE = 1000;
91+
try
92+
if do_enkf
93+
meth.analysis(R, y);
94+
end
95+
catch
96+
do_enkf = false;
17897
end
17998

99+
xa = meth.BestEstimate;
180100

181-
rmses(ensNi, infi) = resE;
101+
err = xt - xa;
182102

183-
mm = min(rmses(:));
184-
if mm >= maxallowerr
185-
mm = 0;
103+
if i > spinup
104+
105+
sse = sse + sum((xa - xt).^2);
106+
rmse = sqrt(sse/(i - spinup)/nvrs);
107+
fprintf('step %d %.5f\n', i, rmse);
108+
109+
else
110+
111+
fprintf('step %d\n', i)
112+
186113
end
187114

188-
imagesc(ensNs, infs, rmses.'); caxis([mm, 1]); colorbar; set(gca,'YDir','normal');
189-
axis square; title('EnKF'); colormap('pink');
190-
xlabel('Ensemble Size'); ylabel('Inflation');
191-
drawnow;
192115
end
193116

194-
return;
117+
rmse
118+

0 commit comments

Comments
 (0)