Skip to content

Commit be5b3e9

Browse files
author
Abhinab Bhattacharjee
committed
resolved first set of conflicts
2 parents faaf1f7 + 19717f0 commit be5b3e9

File tree

7 files changed

+284
-7
lines changed

7 files changed

+284
-7
lines changed

‎.gitignore

+4
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,10 @@
22
*.asv
33
*.png
44
*.fig
5+
<<<<<<< HEAD
56
*.pdf
67
*.mat
78
*.html
9+
=======
10+
.DS_Store
11+
>>>>>>> master

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

+112
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
clear;
2+
close all;
3+
4+
rng(5);
5+
6+
Deltat = 0.12;
7+
8+
solvermodel = @(f, t, y) datools.utils.rk4(f, t, y, 100);
9+
solvernature = @(f, t, y) datools.utils.rk4(f, t, y, 100);
10+
11+
% solvermodel = @(f, t, y) ode45(f, t, y);
12+
% solvernature = @(f, t, y) ode45(f, t, y);
13+
14+
natureODE = otp.lorenz63.presets.Canonical;
15+
natureODE.TimeSpan = [0, Deltat];
16+
17+
nvrs = natureODE.NumVars;
18+
19+
nature0 = randn(nvrs, 1);
20+
21+
[tt, yy] = ode45(natureODE.RHS.F, [0, 10], nature0);
22+
natureODE.Y0 = yy(end, :).';
23+
24+
modelODE = otp.lorenz63.presets.Canonical;
25+
modelODE.TimeSpan = [0, Deltat];
26+
27+
28+
model = datools.Model('Solver', solvermodel, 'ODEModel', modelODE);
29+
nature = datools.Model('Solver', solvernature, 'ODEModel', natureODE);
30+
31+
% naturetomodel = datools.observation.Linear(numel(nature0), 'H', speye(nvrs));
32+
naturetomodel = datools.observation.Indexed(numel(nature0), 'Indices', 1:nvrs);
33+
34+
observeindicies = 1:3;
35+
36+
nobsvars = numel(observeindicies);
37+
38+
R = (8 / 1) * speye(nobsvars);
39+
dR = decomposition(R, 'chol');
40+
41+
obserrormodel = datools.error.Gaussian('CovarianceSqrt', sqrtm(R));
42+
observation = datools.observation.Indexed(model.NumVars, ...
43+
'ErrorModel', obserrormodel, ...
44+
'Indices', observeindicies);
45+
46+
% We make the assumption that there is no model error
47+
modelerror = datools.error.Error;
48+
49+
% No localization
50+
localization = [];
51+
52+
B = [0.86, 0.86, -0.02; ...
53+
0.86, 1.1149, -0.01; ...
54+
-0.02, -0.01, 1.02];
55+
56+
B = 8*B;
57+
58+
meth = datools.variational.ThreeDVar(model, ...
59+
'ModelError', modelerror, ...
60+
'Observation', observation, ...
61+
'InitialState', nature.State, ...
62+
'BackgroundCovariance', B, ...
63+
'OptimizationType', 'lbfgs');
64+
65+
spinup = 200;
66+
times = 2200;
67+
68+
do_filter = true;
69+
70+
sse = 0;
71+
rmse = 0;
72+
73+
for i = 1:times
74+
75+
nature.evolve();
76+
77+
if do_filter
78+
meth.forecast();
79+
end
80+
81+
82+
% observe
83+
xt = naturetomodel.observeWithoutError(nature.TimeSpan(1), nature.State);
84+
y = meth.Observation.observeWithError(model.TimeSpan(1), xt);
85+
86+
% analysis
87+
%try
88+
if do_filter
89+
meth.analysis(dR, y);
90+
end
91+
%catch
92+
% do_enkf = false;
93+
%end
94+
95+
xa = meth.BestEstimate;
96+
97+
err = xt - xa;
98+
99+
if i > spinup
100+
101+
sse = sse + sum((xa - xt).^2);
102+
rmse = sqrt(sse/(i - spinup)/nvrs);
103+
104+
end
105+
106+
if rem(i, 10) == 0
107+
fprintf('step = %d rmse = %.5f\n', i, rmse);
108+
end
109+
110+
end
111+
112+
rmse

‎src/+datools/+statistical/+ensemble/ETKF.m

+4-4
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ function analysis(obj)
2020

2121
xfm = mean(xf, 2);
2222
Af = xf - repmat(xfm, 1, ensN);
23-
Af = inflation * Af/sqrt(ensN - 1);
23+
Af = inflation*Af/sqrt(ensN - 1);
2424
xf = repmat(xfm, 1, ensN) + Af*sqrt(ensN - 1);
2525

2626
Hxf = obj.ObservationOperator.observeWithoutError(tc, xf);
@@ -29,18 +29,18 @@ function analysis(obj)
2929
HAf = HAf/sqrt(ensN - 1);
3030

3131
dS = decomposition((HAf * HAf.')+ R, 'chol');
32+
3233
dR = decomposition(R, 'chol');
3334

34-
temp = dS \ HAf;
35-
T = sqrtm(eye(ensN)-(HAf.' * temp));
35+
T = sqrtm(eye(ensN) - (HAf.'*(dS\HAf)));
3636

3737
Aa = Af * T;
3838
xam = xfm + ((Aa * (HAf * T).') * (dR \ (obj.Observation.Y - Hxfm)));
3939
xa = sqrt(ensN-1) .* Aa + repmat(xam, 1, ensN);
4040

4141
obj.Ensemble = xa;
4242

43-
obj.Weights = ones(ensN, 1) / ensN;
43+
obj.Weights = ones(ensN, 1)/ensN;
4444

4545
end
4646

‎src/+datools/+statistical/+ensemble/EnF.m

+2-2
Original file line numberDiff line numberDiff line change
@@ -48,11 +48,11 @@
4848
%addParameter(p, 'ModelError', datools.error.Error);
4949
addParameter(p, 'NumEnsemble', 1);
5050
addParameter(p, 'Inflation', 1);
51-
addParameter(p, 'Rejuvenation', []);
51+
addParameter(p, 'Rejuvenation', 0);
5252
addParameter(p, 'Localization', []);
5353
addParameter(p, 'Parallel', false);
5454
addParameter(p, 'RIPIterations', 0);
55-
addParameter(p, 'RankHistogram', [])
55+
addParameter(p, 'RankHistogram', []);
5656
addParameter(p, 'ResamplingThreshold', 0.5);
5757
addParameter(p, 'EnsembleGenerator', @(x, y) randn(x, y));
5858

‎src/+datools/+statistical/+ensemble/EnKF.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ function analysis(obj)
3939
else
4040
H = obj.ObservationOperator.linearization(tc, xfm);
4141
rhoHt = obj.Localization(tc, xfm, H);
42-
HrhoHt = H * rhoHt;
42+
HrhoHt = eye(size(H)) * rhoHt;
4343
end
4444

4545
PfHt = rhoHt .* ((1 / (ensN - 1)) * (Af * (HAf.')));

‎src/+datools/+tapering/gaussRloc.m

+14
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
function Ctilde = gaussRloc(t, y, H, r, d, k)
2+
3+
ki = d(t, y, k, H) / r;
4+
CtildeD = gauss(ki).';
5+
6+
Ctilde = spdiags(CtildeD, 0, numel(H), numel(H));
7+
8+
end
9+
10+
function g = gauss(k)
11+
12+
g = exp(-(1/2)*(k.^2));
13+
14+
end

‎src/+datools/+variational/ThreeDVar.m

+147
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
classdef ThreeDVar < handle
2+
%
3+
4+
properties
5+
Model % type of ODE solver (ode45/Runge Kutta) and the model (eg: Lorenz63)
6+
ModelError % type err
7+
Observation % type of obervation
8+
State % current state values for all the states
9+
B % background covariance
10+
OptAlg % lbfgs or newton
11+
end
12+
13+
properties (Dependent)
14+
BestEstimate % Current estimate of the particles/ensembles
15+
end
16+
17+
properties (Access=protected)
18+
BDecomposition
19+
end
20+
21+
methods
22+
function obj = ThreeDVar(varargin)
23+
%ENF The constructor initializes the properties/attributes
24+
%
25+
% OBJ = ENF(VARARGIN) accepts variable length argument list
26+
% VARARGIN and updates the properties/attributes of the object
27+
% (OBJ) of this class or a derived class
28+
29+
optAlgValFcn = @(x) (strcmp(x, 'lbfgs') || strcmp(x, 'newton'));
30+
31+
p = inputParser;
32+
p.KeepUnmatched = true;
33+
addRequired(p, 'Model', @(x) isa(x, 'datools.Model'));
34+
addParameter(p, 'ModelError', datools.error.Error);
35+
addParameter(p, 'InitialState', []);
36+
addParameter(p, 'BackgroundCovariance', []);
37+
addParameter(p, 'OptimizationType', 'lbfgs', optAlgValFcn);
38+
parse(p, varargin{:});
39+
40+
s = p.Results;
41+
42+
obj.Model = s.Model;
43+
obj.ModelError = s.ModelError;
44+
obj.State = s.InitialState;
45+
obj.B = s.BackgroundCovariance;
46+
obj.BDecomposition = decomposition(obj.B, 'chol');
47+
obj.OptAlg = s.OptimizationType;
48+
49+
kept = p.Unmatched;
50+
51+
p = inputParser;
52+
addParameter(p, 'Observation', datools.observation.Observation(s.Model.NumVars));
53+
parse(p, kept);
54+
55+
s = p.Results;
56+
57+
obj.Observation = s.Observation;
58+
end
59+
60+
function forecast(obj)
61+
%FORECAST Method to propagate the model forward in time
62+
%
63+
% FORECAST(OBJ) propoagates the model one step in time
64+
% using a user defined time integration method
65+
66+
[~ , yend] = obj.Model.solve([], obj.State);
67+
68+
obj.State = obj.ModelError.adderr(obj.Model.TimeSpan(end), yend);
69+
70+
71+
end
72+
73+
function analysis(obj, R, y)
74+
75+
% A constrained problem like double pendulum will need a
76+
% constraint coming in from the object.
77+
78+
dB = obj.BDecomposition;
79+
80+
if strcmp(class(R), "decomposition")
81+
dR = R;
82+
else
83+
dR = decomposition(R, 'chol');
84+
end
85+
86+
xb = obj.State;
87+
obs = obj.Observation;
88+
89+
t = obj.Model.TimeSpan(end);
90+
91+
H = @(x) obs.observeWithoutError(t, x);
92+
Hadjoint = @(x) obs.linearization(t, x)';
93+
94+
J = @(x) cost(x, xb, dB, y, dR, H, Hadjoint);
95+
96+
if strcmp(obj.OptAlg, 'lbfgs')
97+
opts = optimoptions('fmincon','Display','none', ...
98+
'SpecifyObjectiveGradient',true, ...
99+
'HessianApproximation', {'lbfgs', 50});
100+
else
101+
102+
hessv = @(Hax, v) dB\v + Hax*(dR\(Hax'*v));
103+
opts = optimoptions('fmincon','Display','none', ...
104+
'Algorithm','trust-region-reflective', ...
105+
'SpecifyObjectiveGradient',true, ...
106+
'HessianMultiplyFcn', hessv, ...
107+
'SubproblemAlgorithm', 'cg');
108+
109+
end
110+
111+
xa = fmincon(J, xb, [], [], [], [], [], [], [], opts);
112+
113+
obj.State = xa;
114+
115+
function [c, g, hessinfo] = cost(x, xb, dB, y, dR, H, Hadjoint)
116+
117+
dx = x - xb;
118+
dy = H(x) - y;
119+
Hax = Hadjoint(x);
120+
121+
Bix = dB\dx;
122+
Riy = dR\dy;
123+
124+
c = 0.5*((dx'*Bix) + (dy'*Riy));
125+
g = Bix + Hax*Riy;
126+
hessinfo = Hax;
127+
128+
end
129+
130+
131+
end
132+
133+
function x = get.BestEstimate(obj)
134+
%GET.BESTESTIMATE Method to estimate of ensemble values
135+
%
136+
% X = GET.BESTESTIMATE(OBJ) uses an in-built getter method,
137+
% derived from handle class, to return the best estimate of
138+
% the information from the current ensembles of states and
139+
% its corresponding weights
140+
141+
x = obj.State;
142+
143+
end
144+
145+
end
146+
147+
end

0 commit comments

Comments
 (0)