Skip to content

Commit 4fefc9d

Browse files
author
Abhinab Bhattacharjee
committed
Merge branch 'abhinab' of https://github.com/ComputationalScienceLaboratory/DATools into abhinab
2 parents 3ad3e3f + b94b5af commit 4fefc9d

File tree

2 files changed

+308
-1
lines changed

2 files changed

+308
-1
lines changed

‎src/+datools/+examples/+sandu/runnerscript.m

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717
% filtername = 'RHF';
1818

1919
%% User inputs
20-
user.modelname = 'Lorenz63';
20+
user.modelname = 'Lorenz96';
2121

2222
user.filtername = 'EnKF';
2323

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

Lines changed: 307 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,307 @@
1+
clear;
2+
close all;
3+
4+
% time steps
5+
Deltat = 0.0109;
6+
7+
% 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+
11+
% Define ODE
12+
natureODE = otp.qg.presets.Canonical('Size', [63, 127]);
13+
nature0 = natureODE.Y0;
14+
natureODE.TimeSpan = [0, Deltat];
15+
16+
modelODE = otp.qg.presets.Canonical('Size', [63, 127]);
17+
modelODE.TimeSpan = [0, Deltat];
18+
19+
% initialize model
20+
model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
21+
nature = datools.Model('Solver', solvernature, 'ODEModel', natureODE);
22+
23+
% Observation Model
24+
naturetomodel = datools.observation.Linear(numel(nature0), 'H', speye(natureODE.NumVars));
25+
26+
%observeindicies = 1:natureODE.NumVars;
27+
nobsvars = 150;
28+
observeindicies = round(linspace(1, natureODE.NumVars, nobsvars));
29+
30+
R = (1 / 1) * speye(nobsvars);
31+
32+
obserrormodel = datools.error.Gaussian('CovarianceSqrt', sqrtm(R));
33+
%obserrormodel = datools.error.Tent;
34+
observation = datools.observation.Indexed(model.NumVars, ...
35+
'ErrorModel', obserrormodel, ...
36+
'Indices', observeindicies);
37+
38+
39+
% We make the assumption that there is no model error
40+
modelerror = datools.error.Error;
41+
42+
%ensembleGenerator = @(N) randn(natureODE.NumVars, N);
43+
44+
load('qgtrajectory.mat')
45+
y = y.';
46+
Nt = size(y, 2);
47+
48+
ensembleGenerator = @(N) natureODE.Y0 + 0.1*y(:, randperm(Nt, N));
49+
50+
% ensNs = 100:100:500;
51+
% infs = 1.01:.01:1.05;
52+
53+
ensNs = [15, 25, 50, 100];
54+
infs = [1.01, 1.02, 1.05, 1.10];
55+
rejs = [1.01, 1.02, 1.03, 1.05];
56+
57+
% variables for which you need the rank histogram plot
58+
histvar = 1:1:nobsvars;
59+
rhplotval = inf * ones(numel(ensNs), numel(infs));
60+
61+
rmses = inf * ones(numel(ensNs), numel(infs));
62+
63+
maxallowerr = 10;
64+
65+
mm = min(rmses(:));
66+
67+
if mm >= maxallowerr
68+
mm = 0;
69+
end
70+
71+
% imagesc(ensNs, infs, rmses.'); caxis([mm, 1]); colorbar; set(gca,'YDir','normal');
72+
% axis square; title('ETKF'); colormap('hot');
73+
% xlabel('Ensemble Size'); ylabel('Inflation');
74+
75+
runsleft = find(rmses == inf);
76+
77+
f1 = figure;
78+
f2 = figure;
79+
f3 = figure;
80+
f4 = figure;
81+
82+
for runn = runsleft.'
83+
[ensNi, infi] = ind2sub([numel(ensNs), numel(infs)], runn);
84+
85+
fprintf('N: %d, inf: %.3f\n', ensNs(ensNi), infs(infi));
86+
87+
ns = 1;
88+
sE = zeros(ns, 1);
89+
90+
inflationAll = infs(infi);
91+
ensN = ensNs(ensNi);
92+
93+
for sample = 1:ns
94+
% Set rng for standard experiments
95+
rng(17+sample-1);
96+
97+
inflation = inflationAll;
98+
99+
% No localization
100+
101+
r = 4;
102+
d = @(t, y, i, j) modelODE.DistanceFunction(t, y, i, j);
103+
localization = [];
104+
105+
106+
%localization = @(t, y, H) datools.tapering.gc(t, y, r, d, H);
107+
%$localization = @(t, y, Hi, k) datools.tapering.gcCTilde(t, y, Hi, r, d, k);
108+
%localization = @(t, y, Hi, k) datools.tapering.cutoffCTilde(t, y, r, d, Hi, k);
109+
110+
111+
enkf = datools.statistical.ensemble.EnKF(model, ...
112+
'Observation', observation, ...
113+
'NumEnsemble', ensN, ...
114+
'ModelError', modelerror, ...
115+
'EnsembleGenerator', ensembleGenerator, ...
116+
'Inflation', inflation, ...
117+
'Localization', localization, ...
118+
'Parallel', false, ...
119+
'RankHistogram', histvar, ...
120+
'Rejuvenation', 0.1);
121+
122+
enkf.setMean(natureODE.Y0);
123+
enkf.scaleAnomalies(1/10);
124+
125+
% define steps and spinups
126+
spinup = 50;
127+
times = 350;
128+
129+
mses = zeros(times - spinup, 1);
130+
131+
rmse = nan;
132+
133+
ps = '';
134+
135+
do_enkf = true;
136+
137+
rmstempval = NaN * ones(1, times-spinup);
138+
139+
for i = 1:times
140+
% forecast
141+
nature.evolve();
142+
143+
if do_enkf
144+
enkf.forecast();
145+
end
146+
147+
148+
% observe
149+
xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
150+
y = enkf.Observation.observeWithError(model.TimeSpan(1), xt);
151+
152+
% Rank histogram (if needed)
153+
datools.utils.stat.RH(enkf, xt);
154+
155+
% analysis
156+
157+
% try
158+
if do_enkf
159+
enkf.analysis(R, y);
160+
end
161+
%catch
162+
% do_enkf = false;
163+
%end
164+
165+
xa = enkf.BestEstimate;
166+
167+
err = xt - xa;
168+
169+
if i > spinup
170+
171+
mses(i - spinup) = mean((xa - xt).^2);
172+
rmse = sqrt(mean(mses(1:(i - spinup))));
173+
174+
rmstempval(i - spinup) = rmse;
175+
176+
% if rmse > maxallowerr || isnan(rmse) || mses(i - spinup) > 2*maxallowerr
177+
% do_enkf = false;
178+
% end
179+
end
180+
181+
182+
if ~do_enkf
183+
break;
184+
end
185+
186+
end
187+
188+
if isnan(rmse)
189+
rmse = 1000;
190+
end
191+
192+
if ~do_enkf
193+
sE(sample) = 1000;
194+
else
195+
sE(sample) = rmse;
196+
end
197+
198+
end
199+
200+
resE = mean(sE);
201+
202+
if isnan(resE)
203+
resE = 1000;
204+
end
205+
206+
rmses(ensNi, infi) = resE;
207+
208+
[xs, pval, rhplotval(ensNi, infi)] = datools.utils.stat.KLDiv(enkf.RankValue(1, 1:end-1), ...
209+
(1 / ensN)*ones(1, ensN+1));
210+
211+
mm = min(rmses(:));
212+
mm = 0;
213+
214+
if mm >= maxallowerr
215+
mm = 0;
216+
end
217+
218+
rw = numel(infs) - 1 - floor((runn - 1)/numel(ensNs));
219+
cl = runn - floor((runn - 1)/numel(ensNs)) * numel(ensNs);
220+
figure(f1);
221+
set(gca, 'XTick', linspace(ensNs(1), ensNs(end), size(ensNs, 2)));
222+
set(gca, 'XTickLabel', ensNs);
223+
subplot(numel(infs), numel(ensNs), rw*numel(ensNs)+cl);
224+
hold all;
225+
z = enkf.RankValue(1, 1:end-1);
226+
z = z / sum(z);
227+
NN = numel(z);
228+
z = NN * z;
229+
bar(xs, z);
230+
plot(xs, pval, '-*r');
231+
set(gca, 'XTick', [xs(1), xs(end)]);
232+
set(gca, 'XTickLabel', [1, ensN + 1]);
233+
set(gca, 'YTick', []);
234+
set(gca, 'YTickLabel', []);
235+
han = axes(f1, 'visible', 'off');
236+
han.Title.Visible = 'on';
237+
han.XLabel.Visible = 'on';
238+
han.YLabel.Visible = 'on';
239+
%han.XTick.Visible = 'on';
240+
%han.XTickLabel.Visible = 'on';
241+
%set(han, 'XTick', linspace(ensNs(1), ensNs(end), size(ensNs,2)));
242+
%set(han, 'XTickLabel', ensNs);
243+
ylabel(han, 'Inflation');
244+
xlabel(han, 'Ensemble Size');
245+
title(han, 'Rank Histogram');
246+
drawnow;
247+
248+
figure(f2);
249+
imagesc(ensNs, infs, rmses.');
250+
caxis([0, 1]);
251+
colorbar;
252+
set(gca, 'YDir', 'normal');
253+
set(gca, 'XTick', linspace(ensNs(1), ensNs(end), size(ensNs, 2)));
254+
set(gca, 'XTickLabel', ensNs);
255+
set(gca, 'YTick', linspace(infs(1), infs(end), size(infs, 2)));
256+
set(gca, 'YTickLabel', infs);
257+
axis square;
258+
title('Rmse HeatMap');
259+
colormap('pink');
260+
xlabel('Ensemble Size');
261+
ylabel('Inflation');
262+
drawnow;
263+
264+
figure(f3);
265+
map = bone;
266+
map = map(1:2:end-1, :);
267+
pt = flipud(pink);
268+
map = [map; pt(2:2:end, :)];
269+
imagesc(ensNs, infs, rhplotval.');
270+
caxis([-0.1, 0.1]);
271+
colorbar;
272+
set(gca, 'YDir', 'normal');
273+
set(gca, 'XTick', linspace(ensNs(1), ensNs(end), size(ensNs, 2)));
274+
set(gca, 'XTickLabel', ensNs);
275+
set(gca, 'YTick', linspace(infs(1), infs(end), size(infs, 2)));
276+
set(gca, 'YTickLabel', infs);
277+
axis square;
278+
title('KLDiv');
279+
colormap(map);
280+
xlabel('Ensemble Size');
281+
ylabel('Inflation');
282+
drawnow;
283+
284+
figure(f4);
285+
subplot(numel(infs), numel(ensNs), rw*numel(ensNs)+cl);
286+
plot(spinup+1:1:times, rmstempval);
287+
xlim([spinup + 1, times]);
288+
ylim([0, 1]);
289+
set(gca, 'XTick', [spinup + 1, times])
290+
set(gca, 'XTickLabel', [spinup + 1, times])
291+
han = axes(f4, 'visible', 'off');
292+
han.Title.Visible = 'on';
293+
han.XLabel.Visible = 'on';
294+
han.YLabel.Visible = 'on';
295+
ylabel(han, 'Value');
296+
xlabel(han, 'Time Steps');
297+
title(han, 'RMSE');
298+
drawnow;
299+
300+
% imagesc(ensNs, infs, rmses.'); caxis([mm, 1]); colorbar; set(gca,'YDir','normal');
301+
% axis square; title('EnKF'); colormap('pink');
302+
% xlabel('Ensemble Size'); ylabel('Inflation');
303+
% drawnow;
304+
end
305+
% savdir = '/home/abhinab93/Documents/experiments/Lorenz96/EnKFr4/l96enkf.mat';
306+
% save(fullfile(savdir));
307+
return;

0 commit comments

Comments
 (0)