Multivariate_Localization_Functions/optimal_univariate_nu_save_output.m at main · zcstanley/Multivariate_Localization_Functions · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
%% Saves output from EnKF with different values of nu for Askey and Wendland
%
% Author: Zofia Stanley

%% Set fixed parameters
params = struct('K',36,'J',10,'F',10,'a',10,'b',10,'h',2);  %L96 params
Nt = 3000;              % Number of assimilation cycles
sigma2Y = 0.005;        % Y obs error variance
sigma2X = 0.28;         % X obs error variance
Ne = 20;                % Ensemble size
x_position = 'middle';  % Where is X_k located? (middle or first)
rInf = sqrt(1.21);      % Inflation factor
Adapt_Inf = true;       % Use adaptive inflation or constant inflation?
True_Fcast_Err = false; % Save true forecast error correlations?
Ntrial = 3;             % Number of trials
askey_nu = [1, 1.5, 2, 2.5];    % nu for univariate Askey
wendland_nu = askey_nu+1;       % nu for univariate Wendland
start = 1001;           % first assimilation cycle considered in RMSE computations

%% Run Experiments

% 1. All Y, No X
fprintf('All Y, No X.\n')
dtObs = 0.005;         % Time between assimilation cycles
Frac_Obs_Y = 1;        % Fraction of Y variables that are observed
Frac_Obs_X = 0;        % Fraction of X variables that are observed
R = 15;                % Univariate localization radius

% Initialize arrays
RMSE_Y_YnoX_UVA = zeros(1, length(askey_nu));
RMSE_X_YnoX_UVA = zeros(1, length(askey_nu));
RMSE_Y_YnoX_UVW = zeros(1, length(wendland_nu));
RMSE_X_YnoX_UVW = zeros(1, length(wendland_nu));

% 1A. Askey
fprintf('\nAskey\n')
loc_fun_name = 'askey' ;
loc_params = struct('rYY', R, 'rXX', R, 'rXY', R, 'gammaYY', 0, 'gammaXX', 0, 'gammaXY', 0, 'beta', 1);
for ii = 1:length(askey_nu)
    nu = askey_nu(ii);
    loc_params.nu = nu;
    [RMSE_Y, RMSE_X] = L96_EnKF_Multitrial(params, Nt, dtObs, sigma2Y, sigma2X,...
                                Frac_Obs_Y, Frac_Obs_X, Ne, x_position, rInf, Adapt_Inf, ...
                                loc_fun_name, loc_params, True_Fcast_Err, Ntrial);
    RMSE_Y_YnoX_UVA(ii) = mean(reshape(RMSE_Y(:, start:end), 1, []));
    RMSE_X_YnoX_UVA(ii) = mean(reshape(RMSE_X(:, start:end), 1, []));
end
save('optimal_univariate_nu.mat', 'askey_nu', 'RMSE_Y_YnoX_UVA', 'RMSE_X_YnoX_UVA')

% 1B. Wendland
fprintf('\nWendland\n')
loc_fun_name = 'wendland' ;
loc_params = struct('rYY', R, 'rXX', R, 'rXY', R, 'gammaYY', 0, 'gammaXX', 0, 'gammaXY', 0, 'k', 1, 'beta', 1);
for ii = 1:length(wendland_nu)
    nu = wendland_nu(ii);
    loc_params.nu = nu;
    [RMSE_Y, RMSE_X] = L96_EnKF_Multitrial(params, Nt, dtObs, sigma2Y, sigma2X,...
                                Frac_Obs_Y, Frac_Obs_X, Ne, x_position, rInf, Adapt_Inf, ...
                                loc_fun_name, loc_params, True_Fcast_Err, Ntrial);
    RMSE_Y_YnoX_UVW(ii) = mean(reshape(RMSE_Y(:, start:end), 1, []));
    RMSE_X_YnoX_UVW(ii) = mean(reshape(RMSE_X(:, start:end), 1, []));
end
save('optimal_univariate_nu.mat', 'wendland_nu', 'RMSE_Y_YnoX_UVW', 'RMSE_X_YnoX_UVW', '-append')

% 2. All X, No Y
fprintf('\nAll X, No Y.\n')
dtObs = 0.05;          % Time between assimilation cycles
Frac_Obs_Y = 0;        % Fraction of Y variables that are observed
Frac_Obs_X = 1;        % Fraction of X variables that are observed
R = 20;                % Univariate localization radius

% Initialize arrays
RMSE_Y_XnoY_UVA = zeros(1, length(askey_nu));
RMSE_X_XnoY_UVA = zeros(1, length(askey_nu));
RMSE_Y_XnoY_UVW = zeros(1, length(wendland_nu));
RMSE_X_XnoY_UVW = zeros(1, length(wendland_nu));


% 2A. Askey
fprintf('\nAskey\n')
loc_fun_name = 'askey' ;
loc_params = struct('rYY', R, 'rXX', R, 'rXY', R, 'gammaYY', 0, 'gammaXX', 0, 'gammaXY', 0, 'beta', 1);
for ii = 1:length(askey_nu)
    nu = askey_nu(ii);
    loc_params.nu = nu;
    [RMSE_Y, RMSE_X] = L96_EnKF_Multitrial(params, Nt, dtObs, sigma2Y, sigma2X,...
                                Frac_Obs_Y, Frac_Obs_X, Ne, x_position, rInf, Adapt_Inf, ...
                                loc_fun_name, loc_params, True_Fcast_Err, Ntrial);
    RMSE_Y_XnoY_UVA(ii) = mean(reshape(RMSE_Y(:, start:end), 1, []));
    RMSE_X_XnoY_UVA(ii) = mean(reshape(RMSE_X(:, start:end), 1, []));
end
save('optimal_univariate_nu.mat', 'RMSE_Y_XnoY_UVA', 'RMSE_X_XnoY_UVA', '-append')

% 2B. Wendland
fprintf('\nWendland\n')
loc_fun_name = 'wendland' ;
loc_params = struct('rYY', R, 'rXX', R, 'rXY', R, 'gammaYY', 0, 'gammaXX', 0, 'gammaXY', 0, 'k', 1, 'beta', 1);
for ii = 1:length(wendland_nu)
    nu = wendland_nu(ii);
    loc_params.nu = nu;
    [RMSE_Y, RMSE_X] = L96_EnKF_Multitrial(params, Nt, dtObs, sigma2Y, sigma2X,...
                                Frac_Obs_Y, Frac_Obs_X, Ne, x_position, rInf, Adapt_Inf, ...
                                loc_fun_name, loc_params, True_Fcast_Err, Ntrial);
    RMSE_Y_XnoY_UVW(ii) = mean(reshape(RMSE_Y(:, start:end), 1, []));
    RMSE_X_XnoY_UVW(ii) = mean(reshape(RMSE_X(:, start:end), 1, []));
end
save('optimal_univariate_nu.mat', 'RMSE_Y_XnoY_UVW', 'RMSE_X_XnoY_UVW', '-append')


% 3. Both X and Y
fprintf('\nBoth X and Y.\n')
dtObs = 0.005;          % Time between assimilation cycles
Frac_Obs_Y = 0.75;      % Fraction of Y variables that are observed
Frac_Obs_X = 0.75;      % Fraction of X variables that are observed
sigma2Y = 0.01;         % Y obs error variance
sigma2X = 0.57;         % X obs error variance
R = 15;                 % Univariate localization radius

% Initialize arrays
RMSE_Y_BothXY_UVA = zeros(1, length(askey_nu));
RMSE_X_BothXY_UVA = zeros(1, length(askey_nu));
RMSE_Y_BothXY_UVW = zeros(1, length(wendland_nu));
RMSE_X_BothXY_UVW = zeros(1, length(wendland_nu));

% 3A. Askey
fprintf('\nAskey\n')
loc_fun_name = 'askey' ;
loc_params = struct('rYY', R, 'rXX', R, 'rXY', R, 'gammaYY', 0, 'gammaXX', 0, 'gammaXY', 0, 'beta', 1);
for ii = 1:length(askey_nu)
    nu = askey_nu(ii);
    loc_params.nu = nu;
    [RMSE_Y, RMSE_X] = L96_EnKF_Multitrial(params, Nt, dtObs, sigma2Y, sigma2X,...
                                Frac_Obs_Y, Frac_Obs_X, Ne, x_position, rInf, Adapt_Inf, ...
                                loc_fun_name, loc_params, True_Fcast_Err, Ntrial);
    RMSE_Y_BothXY_UVA(ii) = mean(reshape(RMSE_Y(:, start:end), 1, []));
    RMSE_X_BothXY_UVA(ii) = mean(reshape(RMSE_X(:, start:end), 1, []));
end
save('optimal_univariate_nu.mat', 'RMSE_Y_BothXY_UVA', 'RMSE_X_BothXY_UVA', '-append')

% 3B. Wendland
fprintf('\nWendland\n')
loc_fun_name = 'wendland' ;
loc_params = struct('rYY', R, 'rXX', R, 'rXY', R, 'gammaYY', 0, 'gammaXX', 0, 'gammaXY', 0, 'k', 1, 'beta', 1);
for ii = 1:length(wendland_nu)
    nu = wendland_nu(ii);
    loc_params.nu = nu;
    [RMSE_Y, RMSE_X] = L96_EnKF_Multitrial(params, Nt, dtObs, sigma2Y, sigma2X,...
                                Frac_Obs_Y, Frac_Obs_X, Ne, x_position, rInf, Adapt_Inf, ...
                                loc_fun_name, loc_params, True_Fcast_Err, Ntrial);
    RMSE_Y_BothXY_UVW(ii) = mean(reshape(RMSE_Y(:, start:end), 1, []));
    RMSE_X_BothXY_UVW(ii) = mean(reshape(RMSE_X(:, start:end), 1, []));
end
save('optimal_univariate_nu.mat', 'RMSE_Y_BothXY_UVW', 'RMSE_X_BothXY_UVW', '-append')