[Dart-dev] [7956] DART/trunk/DART_LAB/matlab: Changed the definition of spread and error in the oned_model.m to
nancy at ucar.edu
nancy at ucar.edu
Wed May 6 10:52:43 MDT 2015
Revision: 7956
Author: thoar
Date: 2015-05-06 10:52:42 -0600 (Wed, 06 May 2015)
Log Message:
-----------
Changed the definition of spread and error in the oned_model.m to
be consistent with RMSE and STD - as is used everywhere else.
Almost all the other changes are simply to satisfy M-lint
in preparation for Matlab versions beyond 2014a.
Made the 'green' line a little darker, made the lines (in general) thicker
for improved visibility, plotting, etc.,
Modified Paths:
--------------
DART/trunk/DART_LAB/matlab/gaussian_product.m
DART/trunk/DART_LAB/matlab/get_ens_rank.m
DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m
DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m
DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m
DART/trunk/DART_LAB/matlab/obs_increment_rhf.m
DART/trunk/DART_LAB/matlab/oned_ensemble.fig
DART/trunk/DART_LAB/matlab/oned_ensemble.m
DART/trunk/DART_LAB/matlab/oned_model.m
DART/trunk/DART_LAB/matlab/run_lorenz_63.m
DART/trunk/DART_LAB/matlab/run_lorenz_96.m
DART/trunk/DART_LAB/matlab/twod_ensemble.m
-------------- next part --------------
Modified: DART/trunk/DART_LAB/matlab/gaussian_product.m
===================================================================
--- DART/trunk/DART_LAB/matlab/gaussian_product.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/gaussian_product.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -1,17 +1,17 @@
function varargout = gaussian_product(varargin)
%% GAUSSIAN_PRODUCT demonstrates the product of two gaussian distributions.
%
-% This is fundamental to Kalman filters and to
-% ensemble data assimilation. Change the
-% parameters of the gaussian for the Prior (green) and the Observation
-% (red) and click on 'Plot Posterior'.
+% This is fundamental to Kalman filters and to ensemble
+% data assimilation. Change the parameters of the
+% gaussian for the Prior (green) and the Observation (red)
+% and click on 'Plot Posterior'.
%
-% The product of the two gaussians is a gaussian - in this case,
-% the 'Posterior'. If the parameters of the two gaussians are known,
-% the parameters of the resulting gaussian can be calculated.
+% The product (in this case, the 'Posterior') of two gaussians is a gaussian.
+% If the parameters of the two gaussians are known, the parameters of the
+% resulting gaussian can be calculated.
%
-% See also: oned_model, oned_ensemble, twod_ensemble, run_lorenz_63,
-% run_lorenz_96
+% See also: oned_model, oned_ensemble, twod_ensemble,
+% run_lorenz_63, run_lorenz_96
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -40,7 +40,7 @@
% --- Executes just before gaussian_product is made visible.
-function gaussian_product_OpeningFcn(hObject, eventdata, handles, varargin)
+function gaussian_product_OpeningFcn(hObject, ~, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
@@ -49,6 +49,9 @@
help gaussian_product
+% set random number seed to same value to generate known sequences
+rng('default')
+
% Choose default command line output for gaussian_product
handles.output = hObject;
@@ -64,7 +67,7 @@
% --- Outputs from this function are returned to the command line.
-function varargout = gaussian_product_OutputFcn(hObject, eventdata, handles)
+function varargout = gaussian_product_OutputFcn(~, ~, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
@@ -75,7 +78,7 @@
-function edit1_Callback(hObject, eventdata, handles)
+function edit1_Callback(~, ~, handles)
% hObject handle to edit1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -87,7 +90,7 @@
% --- Executes during object creation, after setting all properties.
-function edit1_CreateFcn(hObject, eventdata, handles)
+function edit1_CreateFcn(hObject, ~, ~)
% hObject handle to edit1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -98,7 +101,7 @@
set(hObject,'BackgroundColor','white');
end
-function edit2_Callback(hObject, eventdata, handles)
+function edit2_Callback(~, ~, handles)
% hObject handle to edit2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -110,7 +113,7 @@
% --- Executes during object creation, after setting all properties.
-function edit2_CreateFcn(hObject, eventdata, handles)
+function edit2_CreateFcn(hObject, ~, ~)
% hObject handle to edit2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -123,7 +126,7 @@
-function edit3_Callback(hObject, eventdata, handles)
+function edit3_Callback(~, ~, handles)
% hObject handle to edit3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -135,7 +138,7 @@
% --- Executes during object creation, after setting all properties.
-function edit3_CreateFcn(hObject, eventdata, handles)
+function edit3_CreateFcn(hObject, ~, ~)
% hObject handle to edit3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -148,7 +151,7 @@
-function edit4_Callback(hObject, eventdata, handles)
+function edit4_Callback(~, ~, handles)
% hObject handle to edit4 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -160,7 +163,7 @@
% --- Executes during object creation, after setting all properties.
-function edit4_CreateFcn(hObject, eventdata, handles)
+function edit4_CreateFcn(hObject, ~, ~)
% hObject handle to edit4 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -172,7 +175,7 @@
end
% --- Executes on button press in pushbutton1.
-function pushbutton1_Callback(hObject, eventdata, handles)
+function pushbutton1_Callback(~, ~, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -180,12 +183,13 @@
% Need to replot prior and obs, then compute posterior and plot
[prior_mean, prior_sd, obs_mean, obs_err_sd, is_err] = g_prod_plot(handles);
-% If there is an error, don't try to do posterior computation
-% But do zero out the posterior text values
-if(is_err) return;
- set(handles.text7, 'String', strcat('Posterior Mean = '));
- set(handles.text8, 'String', strcat('Posterior SD = '));
- set(handles.text9, 'String', strcat('Weight = '));
+% If there is an error, zero out the posterior text values
+% don't try to do posterior computation
+if(is_err)
+ set(handles.text7, 'String', strcat('Posterior Mean = '));
+ set(handles.text8, 'String', strcat('Posterior SD = '));
+ set(handles.text9, 'String', strcat('Weight = '));
+ return;
end
% Compute the posterior mean, sd and weight
Modified: DART/trunk/DART_LAB/matlab/get_ens_rank.m
===================================================================
--- DART/trunk/DART_LAB/matlab/get_ens_rank.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/get_ens_rank.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -8,7 +8,7 @@
% DART $Id$
s_ens = sort(ens);
-rank = max(find(s_ens < squeeze(x))) + 1;;
+rank = find(s_ens < squeeze(x),1,'last') + 1;
if(isempty(rank))
rank = 1;
end
Modified: DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -21,7 +21,7 @@
B = 8/3;
DELTAT = 0.01;
TIME_STEP_DAYS = 0;
-TIMES_STEP_SECONDS = 0;
+TIME_STEP_SECONDS = 0;
% Lorenz-63 fixed model parameters
Modified: DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -46,13 +46,17 @@
global FORCING
global MODEL_SIZE
+dt = zeros(1,MODEL_SIZE);
+
for j = 1:MODEL_SIZE
jp1 = j + 1;
- if(jp1 > MODEL_SIZE) jp1 = 1; end
+ if(jp1 > MODEL_SIZE), jp1 = 1; end
+
jm2 = j - 2;
- if(jm2 < 1) jm2 = MODEL_SIZE + jm2; end
+ if(jm2 < 1), jm2 = MODEL_SIZE + jm2; end
+
jm1 = j - 1;
- if(jm1 < 1) jm1 = MODEL_SIZE; end
+ if(jm1 < 1), jm1 = MODEL_SIZE; end
dt(j) = (x(jp1) - x(jm2)) * x(jm1) - x(j) + FORCING;
end
Modified: DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -17,7 +17,7 @@
FORCING = 8;
DELTA_T = 0.05;
TIME_STEP_DAYS = 0;
-TIMES_STEP_SECONDS = 0;
+TIME_STEP_SECONDS = 0;
% Lorenz-96 fixed model parameters
Modified: DART/trunk/DART_LAB/matlab/obs_increment_rhf.m
===================================================================
--- DART/trunk/DART_LAB/matlab/obs_increment_rhf.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/obs_increment_rhf.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -13,10 +13,15 @@
err = 0;
% Get the ensemble size
-ens_size = size(ensemble, 2);
-prior_sd = std(ensemble);
+ens_size = size(ensemble, 2);
+prior_sd = std(ensemble);
prior_var = prior_sd^2;
+% allocate space for some of the density calculation bits
+like_dense = zeros(1,ens_size); % likelihood density
+mass = zeros(1,ens_size);
+height = zeros(1,ens_size);
+
% Sort the ensemble members and keep the indices
[x, e_ind] = sort(ensemble);
@@ -34,15 +39,15 @@
% Variance of tails is just sample prior variance
% Mean is adjusted so that 1/(ens_size + 1) is outside
left_mean = x(1) + dist_for_unit_sd * prior_sd;
-left_var = prior_var;
-left_sd = prior_sd;
+left_var = prior_var;
+left_sd = prior_sd;
% Same for the right tail
right_mean = x(ens_size) - dist_for_unit_sd * prior_sd;
-right_var = prior_var;
-right_sd = prior_sd;
+right_var = prior_var;
+right_sd = prior_sd;
-% Eventually want to support options, for now
+% Eventually want to support options, for now
gaussian_likelihood_tails = false;
if(gaussian_likelihood_tails)
@@ -50,33 +55,34 @@
else
% Block to do flat tails for likelihood follows
% This removes assumptions about likelihood and cuts cost
- new_var_left = left_var;
- new_sd_left = left_sd;
- new_mean_left = left_mean;
+% new_var_left = left_var;
+ new_sd_left = left_sd;
+ new_mean_left = left_mean;
prod_weight_left = like(1);
- mass(1) = like(1) / (ens_size + 1);
+ mass(1) = like(1) / (ens_size + 1);
% Same for right tail
- new_var_right = right_var;
- new_sd_right = right_sd;
- new_mean_right = right_mean;
- prod_weight_right = like(ens_size);
+% new_var_right = right_var;
+ new_sd_right = right_sd;
+ new_mean_right = right_mean;
+ prod_weight_right = like(ens_size);
mass(ens_size + 1) = like(ens_size) / (ens_size + 1);
% End block for flat tail likelihood
end
% The mass in each interior box is the height times the width
% The height of the likelihood is like_dense
-% For the prior, mass is 1 / ((n+1) width), and mass = height x widt so
+% For the prior, mass is 1 / ((n+1) width), and mass = height x width so
% The height of the prior is 1 / ((n+1) width); Multiplying by width leaves 1/(n+1)
% In prior, have 1/(n+1) mass in each bin, multiply by mean likelihood
% density to get approximate mass in updated bin
+
for i = 2:ens_size
mass(i) = like_dense(i) / (ens_size + 1);
% Height of prior in this bin is mass/width; Only needed for trapezoidal
% If two ensemble members are the same, set height to -1 as flag
- if(x(i) == x(i-1))
+ if(x(i) == x(i-1))
height(i) = -1;
else
height(i) = 1 / ((ens_size + 1) * (x(i) - x(i-1)));
@@ -94,7 +100,7 @@
right_amp = prod_weight_right / mass_sum;
% Find cumulative mass at each box boundary and middle boundary
-cumul_mass(1) = 0;
+cumul_mass = zeros(1,ens_size+1);
for i = 1:ens_size + 1
cumul_mass(i+1) = cumul_mass(i) + nmass(i);
end
@@ -102,6 +108,8 @@
% Begin internal box search at bottom of lowest box
lowest_box = 1;
+new_ens = zeros(1,ens_size);
+obs_increments = zeros(1,ens_size);
% Find each new ensemble member's location
for i = 1:ens_size
% Each update ensemble member has 1/(n+1) mass before it
@@ -164,7 +172,7 @@
% Divide p by alpha to get the right place for weighted normal
np = p / alpha;
-% Find spot in standard normal
+% Find spot in standard normal
x = norm_inv(np);
% Add in the mean and normalize by sd
Modified: DART/trunk/DART_LAB/matlab/oned_ensemble.fig
===================================================================
(Binary files differ)
Modified: DART/trunk/DART_LAB/matlab/oned_ensemble.m
===================================================================
--- DART/trunk/DART_LAB/matlab/oned_ensemble.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/oned_ensemble.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -1,35 +1,35 @@
function varargout = oned_ensemble(varargin)
%% ONED_ENSEMBLE explore the details of ensemble data assimilation for a scalar.
%
-% Click on the 'Create New Ensemble' button to activate the interactive
+% Click on the 'Create New Ensemble' button to activate the interactive
% observation generation mechanism and lay down a set of 'observations'
-% representative of your ensemble. (Think: Some H() operator has
+% representative of your ensemble. (Think: Some H() operator has
% converted the model state to an expected observation.) This is done by
-% placing the cursor near the axis in the plot and clicking. When you
+% placing the cursor near the axis in the plot and clicking. When you
% have all the ensemble members you want, click in the grey area of
% the window outside of the white axis plot.
%
-% After you have an ensemble and an observation, click 'Update Ensemble'.
-% The algorithm is applied and the Posterior (blue) is plotted below the
+% After you have an ensemble and an observation, click 'Update Ensemble'.
+% The algorithm is applied and the Posterior (blue) is plotted below the
% Prior (green). The mean and standard deviation of the posterior are
% also printed on the plot.
-%
-% The type of ensemble Kalman filter update can be chosen using the
-% pulldown menu at the bottom.
%
+% The type of ensemble Kalman filter update can be chosen using the
+% pulldown menu at the bottom.
+%
% Checking the 'Show Inflation' box will also apply inflation to the
% prior before doing the update and will print the mean and standard
% deviation of the inflated prior and the resulting posterior. The
-% inflated prior and posterior are plotted on an axis below the
+% inflated prior and posterior are plotted on an axis below the
% axis for the uninflated ensemble.
%
% The 'EAKF' is a stochastic algorithm so repeated updates can be done
-% for the same prior and observation.
+% for the same prior and observation.
%
% change the Observation Error SD, lay down an ensemble pretty far away
% from the observation - have fun with it.
%
-% See also: gaussian_product, oned_model, twod_ensemble, run_lorenz_63,
+% See also: gaussian_product, oned_model, twod_ensemble, run_lorenz_63,
% run_lorenz_96
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
@@ -43,11 +43,11 @@
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
- 'gui_Singleton', gui_Singleton, ...
- 'gui_OpeningFcn', @oned_ensemble_OpeningFcn, ...
- 'gui_OutputFcn', @oned_ensemble_OutputFcn, ...
- 'gui_LayoutFcn', [] , ...
- 'gui_Callback', []);
+ 'gui_Singleton', gui_Singleton, ...
+ 'gui_OpeningFcn', @oned_ensemble_OpeningFcn, ...
+ 'gui_OutputFcn', @oned_ensemble_OutputFcn, ...
+ 'gui_LayoutFcn', [] , ...
+ 'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
@@ -61,7 +61,7 @@
% --- Executes just before oned_ensemble is made visible.
-function oned_ensemble_OpeningFcn(hObject, eventdata, handles, varargin)
+function oned_ensemble_OpeningFcn(hObject, ~, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
@@ -70,6 +70,9 @@
help oned_ensemble
+% set random number seed to same value to generate known sequences
+rng('default')
+
% Choose default command line output for oned_ensemble
handles.output = hObject;
@@ -103,19 +106,21 @@
set(handles.h_obs_plot, 'Color', 'r', 'Linestyle', '--', 'Linewidth', 2);
hold on
-% Plot an asterisk
-handles.h_obs_ast = plot(handles.observation, 0, 'r*', 'MarkerSize', 16);
+% Plot an asterisk
+handles.h_obs_ast = plot(handles.observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0);
% Set a basic plotting domain range that includes mean +/- 3 obs SDs
-lower = handles.observation - 3*handles.obs_error_sd;
-upper = handles.observation + 3*handles.obs_error_sd;
-axis([lower upper -0.4 1]);
+xlower = handles.observation - 3*handles.obs_error_sd;
+xupper = handles.observation + 3*handles.obs_error_sd;
+ylower = -0.4;
+yupper = 1.0;
+axis([xlower xupper ylower yupper]);
set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]);
set(gca, 'YTickLabel', [0 0.2 0.4 0.6 0.8]);
hold on
-plot([lower upper], [0 0], 'k', 'Linewidth', 2);
+plot([xlower xupper], [0 0], 'k', 'Linewidth', 2);
% Update handles structure
guidata(hObject, handles);
@@ -133,7 +138,7 @@
% --- Outputs from this function are returned to the command line.
-function varargout = oned_ensemble_OutputFcn(hObject, eventdata, handles)
+function varargout = oned_ensemble_OutputFcn(~, ~, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
@@ -147,7 +152,7 @@
% --- Executes on button press in pushbutton_create_new.
-function pushbutton_create_new_Callback(hObject, eventdata, handles)
+function pushbutton_create_new_Callback(hObject, ~, handles)
% hObject handle to pushbutton_create_new (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -176,69 +181,56 @@
hold on
% Set a basic plotting domain range that includes mean +/- 3 obs SDs
-lower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
-upper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
-axis([lower upper -0.4 1]);
+xlower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
+xupper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
+ylower = -0.4;
+yupper = 1.0;
+axis([xlower xupper ylower yupper]);
set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]);
set(gca, 'YTickLabel', [0 0.2 0.4 0.6 0.8]);
-% Messages should start 1/10 of the way across the screen
-x_message = lower + 0.1 * (upper - lower);
-h_click = text(x_message, 0.4, 'Click on x-axis to create member', 'FontSize', 16);
+% Messages are centered in the middle.
+xmid = (xupper + xlower) / 2.0;
+h_click = text(xmid, 0.8, {'Click inside graphics box to create member',...
+ '(only X value is used)'}, 'FontSize', 16, 'HorizontalAlignment', 'center');
-% Need to guarantee at least 2 ensemble members
-ens_size = 0;
+h_err_text = text(xmid, -0.15, 'An ensemble has to have at least 2 members.', ...
+ 'FontSize', 16, 'Visible', 'on', 'HorizontalAlignment', 'center','Color', 'r');
-h_err_text = text(x_message, 0.5, 'Click inside graphics box to select member', ...
- 'Color', 'r', 'FontSize', 16, 'Visible', 'off');
+h_finish = text(xmid, -0.15, 'Click outside of plot to finish', ...
+ 'Fontsize', 16, 'Visible', 'off', 'HorizontalAlignment', 'center');
-while ens_size < 2
- [xt, zt] = ginput(1);
+ens_size = 0;
+while ens_size < 100
+ [xt, yt] = ginput(1);
- % Clear out the error message if it's been made visible
- set(h_err_text, 'Visible', 'off');
+ if(xt >= xlower && xt <= xupper && yt >= ylower && yt <= yupper)
+ ens_size = ens_size + 1;
+ x(ens_size) = xt;
+ y(ens_size) = 0;
+ handles.h_ens_member(ens_size) = ...
+ plot(x(ens_size), y(ens_size), '*', 'MarkerSize', 16, 'Color', [0 0.73 0],'LineWidth',2.0);
- if(xt < lower | xt > upper)
- set(h_err_text, 'Visible', 'on');
- else
- ens_size = ens_size + 1;
- x(ens_size) = xt;
- y(ens_size) = 0;
- handles.h_ens_member(ens_size) = ...
- plot(x(ens_size), y(ens_size), '*', 'MarkerSize', 16, 'Color', [0 0.73 0]);
+ % Display the prior mean and sd
+ prior_mean = mean(x);
+ prior_sd = std(x);
+ set(handles.text2, 'String', ['Prior Mean = ', num2str(prior_mean)]);
+ set(handles.text3, 'String', ['Prior SD = ', num2str(prior_sd)]);
+ elseif (ens_size < 2)
+ set(h_err_text,'FontWeight','bold')
+ else
+ break;
+ end
- % Display the prior mean and sd
- prior_mean = mean(x);
- prior_sd = std(x);
- set(handles.text2, 'String', ['Prior Mean = ', num2str(prior_mean)]);
- set(handles.text3, 'String', ['Prior SD = ', num2str(prior_sd)]);
- end
+ % swap messages once you have a minimal ensemble.
+ if (ens_size == 2)
+ set(h_err_text, 'Visible', 'off');
+ set(h_finish, 'Visible', 'on')
+ end
end
-h_finish = text(x_message, 0.3, 'Click outside of plot to finish', 'Fontsize', 16);
-
-while ens_size < 100
- [xt, zt] = ginput(1);
- % Terminate by clicking outside of graph range
- if(xt > upper | xt < lower)
- break;
- else
- ens_size = ens_size + 1;
- x(ens_size) = xt;
- y(ens_size) = 0;
- handles.h_ens_member(ens_size) = ...
- plot(x(ens_size), y(ens_size), '*', 'MarkerSize', 16, 'Color', [0 0.73 0]);
-
- % Display the prior mean and sd
- prior_mean = mean(x);
- prior_sd = std(x);
- set(handles.text2, 'String', ['Prior Mean = ', num2str(prior_mean)]);
- set(handles.text3, 'String', ['Prior SD = ', num2str(prior_sd)]);
- end
-end
-
-% Ensemble created, comupte mean and sd, clean up and return
+% Ensemble created, compute mean and sd, clean up and return
% Set the global gui storage
handles.ens_size = ens_size;
handles.ens_members = x;
@@ -261,7 +253,7 @@
-function edit_observation_Callback(hObject, eventdata, handles)
+function edit_observation_Callback(hObject, ~, handles)
% hObject handle to edit_observation (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -289,21 +281,21 @@
% Only enable the update ensemble pushbutton if an ensemble has been created
if(handles.ens_size > 0)
- set(handles.pushbutton_update_ens, 'Enable', 'on');
+ set(handles.pushbutton_update_ens, 'Enable', 'on');
end
% Get the value of the observation
if(isfinite( str2double(get(hObject, 'String'))))
- observation = str2double(get(hObject, 'String'));
+ observation = str2double(get(hObject, 'String'));
else
- set(handles.edit_observation, 'String', '???');
+ set(handles.edit_observation, 'String', '???');
- % Disable other input to guarantee only one error at a time!
- set(handles.edit_obs_error_sd, 'Enable', 'off')
- set(handles.edit_inflation, 'Enable', 'off')
- set(handles.pushbutton_create_new, 'Enable', 'off')
- set(handles.pushbutton_update_ens, 'Enable', 'off')
- return
+ % Disable other input to guarantee only one error at a time!
+ set(handles.edit_obs_error_sd, 'Enable', 'off')
+ set(handles.edit_inflation, 'Enable', 'off')
+ set(handles.pushbutton_create_new, 'Enable', 'off')
+ set(handles.pushbutton_update_ens, 'Enable', 'off')
+ return
end
% Update the global storage
@@ -316,18 +308,20 @@
% Move the observation asterisk
set(handles.h_obs_ast, 'Visible', 'Off');
-handles.h_obs_ast = plot(handles.observation, 0, 'r*', 'MarkerSize', 16);
+handles.h_obs_ast = plot(handles.observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0);
% Set a basic plotting domain range that includes mean +/- 3 obs SDs
-lower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
-upper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
-axis([lower upper -0.4 1]);
+xlower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
+xupper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
+ylower = -0.4;
+yupper = 1.0;
+axis([xlower xupper ylower yupper]);
set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]);
set(gca, 'YTickLabel', [0 0.2 0.4 0.6 0.8]);
hold on
-plot([lower upper], [0 0], 'k', 'Linewidth', 2);
+plot([xlower xupper], [0 0], 'k', 'Linewidth', 2);
% Update handles structure
guidata(hObject, handles);
@@ -337,7 +331,7 @@
% --- Executes during object creation, after setting all properties.
-function edit_observation_CreateFcn(hObject, eventdata, handles)
+function edit_observation_CreateFcn(hObject, ~, ~)
% hObject handle to edit_observation (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -352,7 +346,7 @@
%----------------------------------------------------------------------
-function edit_obs_error_sd_Callback(hObject, eventdata, handles)
+function edit_obs_error_sd_Callback(hObject, ~, handles)
% hObject handle to edit_obs_error_sd (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -380,22 +374,22 @@
% Only enable the update ensemble pushbutton if an ensemble has been created
if(handles.ens_size > 0)
- set(handles.pushbutton_update_ens, 'Enable', 'on');
+ set(handles.pushbutton_update_ens, 'Enable', 'on');
end
% Get the value of the observation
if(isfinite( str2double(get(hObject, 'String'))) && ...
- str2double(get(hObject, 'String')) > 0)
- obs_error_sd = str2double(get(hObject, 'String'));
+ str2double(get(hObject, 'String')) > 0)
+ obs_error_sd = str2double(get(hObject, 'String'));
else
- set(handles.edit_obs_error_sd, 'String', '???');
+ set(handles.edit_obs_error_sd, 'String', '???');
- % Disable other input to guarantee only one error at a time!
- set(handles.edit_observation, 'Enable', 'off')
- set(handles.edit_inflation, 'Enable', 'off')
- set(handles.pushbutton_create_new, 'Enable', 'off')
- set(handles.pushbutton_update_ens, 'Enable', 'off')
- return
+ % Disable other input to guarantee only one error at a time!
+ set(handles.edit_observation, 'Enable', 'off')
+ set(handles.edit_inflation, 'Enable', 'off')
+ set(handles.pushbutton_create_new, 'Enable', 'off')
+ set(handles.pushbutton_update_ens, 'Enable', 'off')
+ return
end
% Update the value in global storage
@@ -407,9 +401,11 @@
set(handles.h_obs_plot, 'Color', 'r', 'Linestyle', '--', 'Linewidth', 2);
% Set a basic plotting domain range that includes mean +/- 3 obs SDs
-lower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
-upper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
-axis([lower upper -0.4 1]);
+xlower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
+xupper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
+ylower = -0.4;
+yupper = 1.0;
+axis([xlower xupper ylower yupper]);
set(handles.h_obs_plot, 'Color', 'r', 'Linestyle', '--', 'Linewidth', 2);
@@ -417,7 +413,7 @@
set(gca, 'YTickLabel', [0 0.2 0.4 0.6 0.8]);
hold on
-plot([lower upper], [0 0], 'k', 'Linewidth', 2);
+plot([xlower xupper], [0 0], 'k', 'Linewidth', 2);
% Update handles structure
guidata(hObject, handles);
@@ -427,7 +423,7 @@
% --- Executes during object creation, after setting all properties.
-function edit_obs_error_sd_CreateFcn(hObject, eventdata, handles)
+function edit_obs_error_sd_CreateFcn(hObject, ~, ~)
% hObject handle to edit_obs_error_sd (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -443,7 +439,7 @@
% --- Executes on selection change in popupmenu_filter_kind.
-function popupmenu_filter_kind_Callback(hObject, eventdata, handles)
+function popupmenu_filter_kind_Callback(~, ~, ~)
% hObject handle to popupmenu_filter_kind (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -456,7 +452,7 @@
% --- Executes during object creation, after setting all properties.
-function popupmenu_filter_kind_CreateFcn(hObject, eventdata, handles)
+function popupmenu_filter_kind_CreateFcn(hObject, ~, ~)
% hObject handle to popupmenu_filter_kind (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
@@ -472,7 +468,7 @@
% --- Executes on button press in pushbutton_update_ens.
-function pushbutton_update_ens_Callback(hObject, eventdata, handles)
+function pushbutton_update_ens_Callback(hObject, ~, handles)
% hObject handle to pushbutton_update_ens (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -498,15 +494,15 @@
filter_type = char(h_filter_kind.String(h_filter_kind.Value));
switch filter_type
- case 'EAKF'
- [obs_increments, err] = ...
- obs_increment_eakf(ensemble, handles.observation, handles.obs_error_sd^2);
- case 'EnKF'
- [obs_increments, err] = ...
- obs_increment_enkf(ensemble, handles.observation, handles.obs_error_sd^2);
- case 'RHF'
- [obs_increments, err] = ...
- obs_increment_rhf(ensemble, handles.observation, handles.obs_error_sd^2);
+ case 'EAKF'
+ [obs_increments, ~] = ...
+ obs_increment_eakf(ensemble, handles.observation, handles.obs_error_sd^2);
+ case 'EnKF'
+ [obs_increments, ~] = ...
+ obs_increment_enkf(ensemble, handles.observation, handles.obs_error_sd^2);
+ case 'RHF'
+ [obs_increments, ~] = ...
+ obs_increment_rhf(ensemble, handles.observation, handles.obs_error_sd^2);
end
% Add on increments to get new ensemble
@@ -517,9 +513,9 @@
% Plot lines connecting the prior and posterior ensemble members
for i = 1:size(ensemble, 2)
- x_line = [handles.ens_members(i), new_ensemble(i)];
- y_line = [0, -0.1];
- handles.h_update_lines(i) = plot(x_line, y_line, 'k');
+ x_line = [handles.ens_members(i), new_ensemble(i)];
+ y_line = [0, -0.1];
+ handles.h_update_lines(i) = plot(x_line, y_line, 'k');
end
% Add in a label of the updated mean and sd
@@ -534,18 +530,19 @@
% If the checkbox isn't set, return now
if(not(get(handles.checkbox_inflation, 'Value')))
- guidata(hObject, handles)
- return
+ guidata(hObject, handles)
+ return
end
% Plot the inflated prior ensemble
y = -0.2;
prior_mean = mean(handles.ens_members(1:handles.ens_size));
+inf_ens = zeros(1,handles.ens_size);
for i = 1: handles.ens_size
- inf_ens(i) = (handles.ens_members(i) - prior_mean) * sqrt(handles.inflation) + ...
- prior_mean;
- handles.h_inf_ens_member(i) = plot(inf_ens(i), y, '*', 'MarkerSize', 16, 'Color', [0 0.73 0]);
+ inf_ens(i) = (handles.ens_members(i) - prior_mean) * sqrt(handles.inflation) + ...
+ prior_mean;
+ handles.h_inf_ens_member(i) = plot(inf_ens(i), y, '*', 'MarkerSize', 16, 'Color', [0 0.73 0],'LineWidth',2.0);
end
% Update mean and sd of old posterior
@@ -558,15 +555,15 @@
% Get the update for the inflated ensemble
switch filter_type
- case 'EAKF'
- [obs_increments, err] = ...
- obs_increment_eakf(inf_ens, handles.observation, handles.obs_error_sd^2);
- case 'EnKF'
- [obs_increments, err] = ...
- obs_increment_enkf(inf_ens, handles.observation, handles.obs_error_sd^2);
- case 'RHF'
- [obs_increments, err] = ...
- obs_increment_rhf(inf_ens, handles.observation, handles.obs_error_sd^2);
+ case 'EAKF'
+ [obs_increments, ~] = ...
+ obs_increment_eakf(inf_ens, handles.observation, handles.obs_error_sd^2);
+ case 'EnKF'
+ [obs_increments, ~] = ...
+ obs_increment_enkf(inf_ens, handles.observation, handles.obs_error_sd^2);
+ case 'RHF'
+ [obs_increments, ~] = ...
+ obs_increment_rhf(inf_ens, handles.observation, handles.obs_error_sd^2);
end
% Add on increments to get new ensemble
@@ -577,20 +574,22 @@
% Plot lines connecting the prior and posterior ensemble members
for i = 1:size(ensemble, 2)
- x_line = [inf_ens(i), new_ensemble(i)];
- y_line = [-0.2, -0.3];
- handles.h_inf_lines(i) = plot(x_line, y_line, 'k');
+ x_line = [inf_ens(i), new_ensemble(i)];
+ y_line = [-0.2, -0.3];
+ handles.h_inf_lines(i) = plot(x_line, y_line, 'k');
end
% Set a basic plotting domain range that includes mean +/- 3 obs SDs
% Plus all inflated members
-lower = min(handles.observation - 3*handles.obs_error_sd, min(inf_ens));
-upper = max(handles.observation + 3*handles.obs_error_sd, max(inf_ens));
-axis([lower upper -0.4 1]);
+xlower = min(handles.observation - 3*handles.obs_error_sd, min(inf_ens));
+xupper = max(handles.observation + 3*handles.obs_error_sd, max(inf_ens));
+ylower = -0.4;
+yupper = 1.0;
+axis([xlower xupper ylower yupper]);
% Plot the axes for the two priors
-plot([lower upper], [0 0], 'k', 'Linewidth', 2);
-handles.h_inf_axis = plot([lower upper], [-0.2 -0.2], 'k', 'Linewidth', 2);
+plot([xlower xupper], [0 0], 'k', 'Linewidth', 2);
+handles.h_inf_axis = plot([xlower xupper], [-0.2 -0.2], 'k', 'Linewidth', 2);
% Update mean and sd of old posterior
update_inf_mean = mean(new_ensemble(1:handles.ens_size));
@@ -607,7 +606,7 @@
% --- Executes on button press in checkbox_inflation.
-function checkbox_inflation_Callback(hObject, eventdata, handles)
+function checkbox_inflation_Callback(~, ~, ~)
% hObject handle to checkbox_inflation (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -630,7 +629,7 @@
-function edit_inflation_Callback(hObject, eventdata, handles)
+function edit_inflation_Callback(hObject, ~, handles)
% hObject handle to edit_inflation (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
@@ -658,22 +657,22 @@
% Only enable the update ensemble pushbutton if an ensemble has been created
if(handles.ens_size > 0)
- set(handles.pushbutton_update_ens, 'Enable', 'on');
+ set(handles.pushbutton_update_ens, 'Enable', 'on');
end
% Get the value of the observation
if(isfinite( str2double(get(hObject, 'String'))) && ...
- str2double(get(hObject, 'String')) > 0)
- inflation = str2double(get(hObject, 'String'));
+ str2double(get(hObject, 'String')) > 0)
+ inflation = str2double(get(hObject, 'String'));
else
- set(handles.edit_inflation, 'String', '???');
+ set(handles.edit_inflation, 'String', '???');
- % Disable other input to guarantee only one error at a time!
- set(handles.edit_observation, 'Enable', 'off')
- set(handles.edit_obs_error_sd, 'Enable', 'off')
- set(handles.pushbutton_create_new, 'Enable', 'off')
- set(handles.pushbutton_update_ens, 'Enable', 'off')
- return
+ % Disable other input to guarantee only one error at a time!
+ set(handles.edit_observation, 'Enable', 'off')
+ set(handles.edit_obs_error_sd, 'Enable', 'off')
+ set(handles.pushbutton_create_new, 'Enable', 'off')
+ set(handles.pushbutton_update_ens, 'Enable', 'off')
+ return
end
% Update the value in global storage
@@ -685,9 +684,11 @@
set(handles.h_obs_plot, 'Color', 'r', 'Linestyle', '--', 'Linewidth', 2);
% Set a basic plotting domain range that includes mean +/- 3 obs SDs
-lower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
-upper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
-axis([lower upper -0.4 1]);
+xlower = min(handles.observation - 3*handles.obs_error_sd, min(handles.ens_members));
+xupper = max(handles.observation + 3*handles.obs_error_sd, max(handles.ens_members));
+ylower = -0.4;
+yupper = 1.0;
+axis([xlower xupper ylower yupper]);
set(handles.h_obs_plot, 'Color', 'r', 'Linestyle', '--', 'Linewidth', 2);
@@ -695,7 +696,7 @@
set(gca, 'YTickLabel', [0 0.2 0.4 0.6 0.8]);
hold on
-plot([lower upper], [0 0], 'k', 'Linewidth', 2);
+plot([xlower xupper], [0 0], 'k', 'Linewidth', 2);
% Update handles structure
guidata(hObject, handles);
@@ -707,7 +708,7 @@
% --- Executes during object creation, after setting all properties.
-function edit_inflation_CreateFcn(hObject, eventdata, handles)
+function edit_inflation_CreateFcn(hObject, ~, ~)
% hObject handle to edit_inflation (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
Modified: DART/trunk/DART_LAB/matlab/oned_model.m
===================================================================
--- DART/trunk/DART_LAB/matlab/oned_model.m 2015-05-06 15:09:05 UTC (rev 7955)
+++ DART/trunk/DART_LAB/matlab/oned_model.m 2015-05-06 16:52:42 UTC (rev 7956)
@@ -8,10 +8,10 @@
% ensemble sizes, model biases, etc. on-the-fly. The posterior
% of the state is indicated by blue asterisks, the states evolve along
% a trajectory indicated by the green lines to wind up at a prior state
-% for the assimilation - indicated by the green asterisks. After the
-% assimilation, the (posterior) state is indicated in blue and the
+% for the assimilation - indicated by the green asterisks. After the
+% assimilation, the (posterior) state is indicated in blue and the
% process is ready to repeat.
-%
+%
% ONED_MODEL opens two windows. A gui control window that also plots
% the most recent prior, posterior, and observation, and a figure
% window that plots time sequences of the assimilation, the RMS error,
@@ -23,19 +23,19 @@
%
% Since this is a 'perfect model' experiment, we know the true state,
% the amount of noise added to the observations, etc.; so it is possible to
-% calculate the error of the ensemble in addition to the spread. The
+% calculate the error of the ensemble in addition to the spread. The
% Truth is not (in general) the same as the observation!
%
-% This also introduces the concept of the 'rank histogram'. With N
+% This also introduces the concept of the 'rank histogram'. With N
% ensemble members, there will always be N+1 'bins' that encompass the
% state (to the left of the lowest ensemble member, to the right of the
% highest, and all the ones in-between). The rank histogram tracks the
% number of times the truth lies in each bin. With a good ensemble,
-% (i.e. a good assimilation system) the true state will be
-% indistinguishable from any of the ensemble members. This results
+% (i.e. a good assimilation system) the true state will be
+% indistinguishable from any of the ensemble members. This results
% in a flat rank histogram, given enough samples.
%
-% See also: gaussian_product, oned_ensemble, twod_ensemble, run_lorenz_63,
+% See also: gaussian_product, oned_ensemble, twod_ensemble, run_lorenz_63,
% run_lorenz_96
%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
@@ -49,11 +49,11 @@
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
- 'gui_Singleton', gui_Singleton, ...
- 'gui_OpeningFcn', @oned_model_OpeningFcn, ...
- 'gui_OutputFcn', @oned_model_OutputFcn, ...
- 'gui_LayoutFcn', [] , ...
- 'gui_Callback', []);
+ 'gui_Singleton', gui_Singleton, ...
+ 'gui_OpeningFcn', @oned_model_OpeningFcn, ...
+ 'gui_OutputFcn', @oned_model_OutputFcn, ...
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list