diff --git a/CHANGELOG.rst b/CHANGELOG.rst index ede2aad257..64b31d3ec7 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,12 @@ individual files. The changes are now listed with the most recent at the top. +**January 23 2025 :: DART_LAB QCEFF. Tag v11.10.0** + +- Updated DART_LAB to include QCEFF +- NSF NCAR RDA links updated +- /glade RDA locations updated + **January 23 2025 :: Pangu-DART. Tag v11.9.0** - Pangu-Weather ML model DART interface diff --git a/conf.py b/conf.py index 335f9c7925..f9724da98a 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.9.0' +release = '11.10.0' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/guide/DART_LAB/DART_LAB.rst b/guide/DART_LAB/DART_LAB.rst index 33eeb3be51..b59d261d36 100644 --- a/guide/DART_LAB/DART_LAB.rst +++ b/guide/DART_LAB/DART_LAB.rst @@ -15,11 +15,11 @@ DART tutorial presentations Here are the PDF files for the presentation part of the tutorial: -- :download:`Section 1: ` The basics in 1D. -- :download:`Section 2: ` How should observations of a state variable impact an unobserved - state variable? Multivariate assimilation. -- :download:`Section 3: ` Sampling error and localization. -- :download:`Section 4: ` The Ensemble Kalman Filter (Perturbed Observations). +- :download:`Section 1: ` Ensemble Data Assimilation Concepts in 1D. +- :download:`Section 2: ` How Should Observations Impact an Unobserved + State Variable? Multivariate Assimilation. +- :download:`Section 3: ` Inflation and Localization to Improve Performance. +- :download:`Section 4: ` Nonlinear and Non-Gaussian Extensions. - :download:`Section 5: ` Adaptive Inflation. Matlab hands-on exercises @@ -31,11 +31,16 @@ options. A valid `Matlab `__ license The exercises use the following functions: +- bounded_oned_ensemble - gaussian_product -- oned_model +- oned_cycle - oned_ensemble +- oned_model +- oned_model_inf - run_lorenz_63 - run_lorenz_96 +- run_lorenz_96_inf - twod_ensemble +- twod_ppi_ensemble To run these, cd into the DART_LAB/matlab directory, start matlab, and type the names at the prompt. diff --git a/guide/DART_LAB/matlab/bounded_oned_ensemble.m b/guide/DART_LAB/matlab/bounded_oned_ensemble.m new file mode 100644 index 0000000000..fc27058300 --- /dev/null +++ b/guide/DART_LAB/matlab/bounded_oned_ensemble.m @@ -0,0 +1,1183 @@ +function bounded_oned_ensemble + +%% BOUNDED_ONED_ENSEMBLE explores the details of ensemble data assimilation for a +% nonnegative scalar variable. +% +% Push on the 'Create New Ensemble' button to activate the interactive +% mechanism to generate a prior ensemble estimate of the observation. +% This is done by 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. As you add ensemble members, +% the ensemble prior mean and prior standard deviation are updated in the upper +% left corner of the plot. Members must be nonnegative, and a warning will +% be displayed if an attempt is made to create a negative member. +% +% Once you have created a prior ensemble (green asterisks), click the +% 'Update Ensemble' button. With the default settings, this will apply the EAKF +% algorithm to produce a posterior ensemble (blue asterisks) just below the axis. +% The mean and standard deviation of the posterior are also printed on the plot. +% The EAKF algorithm can produce negative posterior ensemble members. +% +% Two other QCEFF ensemble filter variants can be selected with the pushbuttons +% at the lower right. The gamma filter fits a gamma distribution to the prior +% ensemble and uses a gamma likelihood corresponding to the mean and standard +% deviation of the likelihood. The BNRHF fits a Bounded Normal Rank Histogram +% Prior to the prior ensemble and uses a normal likelihood. +% Selecting one of these and pressing 'Update Ensemble' will produce +% the posterior ensemble and posterior statistics using the selected +% filter algorithm. For all filter types, the prior (green) and posterior (blue) +% continuous distributions are plotted as well as the likelihood function (red). +% +% Checking the 'Apply Inflation' box will also apply inflation to the +% prior ensemble 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 +% axis for the uninflated ensemble. The value of the covariance inflation +% applied can be adjusted using the slider or by typing in a value in the +% 'Inflation Amount' box. The inflation is applied in a probit probability +% integral transformed space. +% +% The mean and standard deviation of the observation likelihood +% can be changed in the red box. +% +% See also: gaussian_product.m oned_cycle.m oned_ensemble.m oned_model.m +% oned_model_inf.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +help bounded_oned_ensemble + +atts = stylesheet; % get the default fonts and colors + +% Set random number seed to same value to generate known sequences +rng('default') + +% Initialize the basic structure we will be using to hold everything. +handles.ens_size = 0; +handles.ens_members = 0; +handles.h_obs_plot = []; +handles.h_update_ens = []; +handles.h_prior_pdf = []; +handles.h_post_pdf = []; +handles.h_ens_member = []; +handles.h_obs_ast = []; +handles.h_update_lines = []; +handles.plot_inflation = false; +handles.h_inf_ens_member = []; +handles.h_inf_up_ens = []; +handles.h_inf_lines = []; +handles.h_inf_axis = []; + +%% ----------------------------------------------------------------------------- + +% Specify the figure size in pixels. After that, all positions are +% specified as fractions (units=Normalized). That way, the objects +% scale proportionally as the figure gets resized. +figXmin = 450; % The horizontal position of the entire figure, in pixels +figYmin = 250; % The vertical position of the entire figure, in pixels +figWidth = 670; % The width of the entire figure, in pixels +figHeight = 500; % The height of the entire figure, in pixels + +handles.figure1 = figure('Position', [figXmin figYmin figWidth figHeight], ... + 'Units', 'Pixels', ... + 'Name', 'bounded_oned_ensemble', ... + 'Color', atts.background); + +%% ----------------------------------------------------------------------------- +%Position has units of normalized and therefore the components must be +%Fractions of figure. Also, the actual text has FontUnits which are +%normalized, so the Font Size is a fraction of the text box/ edit box that +%contains it. Making the font size and box size normalized allows the text +%size and box size to change proportionately if the user resizes the +%figure window. This is true for all text/edit boxes and buttons. + +%% ----------------------------------------------------------------------------- +% Set up a parent container so we can move the one container around instead of +% trying to manipulate the positions of all the components. + +handles.observation_panel = uipanel('BackgroundColor',atts.red, ... + 'BorderType','none', ... + 'Units', 'Normalized', ... + 'Position',[0.66 0.744 0.325 0.2]); + +handles.ui_text_observation = uicontrol(handles.observation_panel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.035 0.600 0.600 0.275], ... + 'String', 'Likelihood mean' , ... + 'BackgroundColor', atts.red, ... + 'ForegroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'FontWeight', 'Bold', ... + 'HorizontalAlignment', 'center'); + +handles.ui_edit_observation = uicontrol(handles.observation_panel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', ... + 'Position', [0.675 0.562 0.300 0.373], ... + 'String', '1', ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'Callback', @edit_observation_Callback); +handles.observation = str2double(get(handles.ui_edit_observation,'String')); + +handles.ui_text_obs_error_sd = uicontrol(handles.observation_panel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.035 0.100 0.600 0.275], ... + 'String', 'Likelihood SD', ... + 'BackgroundColor', atts.red,... + 'ForegroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'center'); + +handles.ui_edit_obs_error_sd = uicontrol(handles.observation_panel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', ... + 'Position', [0.675 0.062 0.300 0.373], ... + 'String', '1', ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'Callback', @edit_obs_error_sd_Callback); +handles.obs_error_sd = str2double(get(handles.ui_edit_obs_error_sd,'String')); + +%% ----------------------------------------------------------------------------- +% Try to center the boxes with what is on top. + +box = get(handles.observation_panel,'Position'); +center = box(1) + box(3)/2.0; +mywid = 0.250; +myleft = center - mywid/2.0; + +handles.ui_button_create_new_ens = uicontrol('Style', 'pushbutton', ... + 'Units', 'Normalized', ... + 'Position', [myleft 0.645 mywid 0.075], ... + 'String', 'Create New Ensemble', ... + 'BackgroundColor', 'White', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'Callback', @button_create_new_ens_Callback); + +mywid = 0.200; +myleft = center - mywid/2.0; +handles.ui_button_update_ens = uicontrol('style', 'pushbutton', ... + 'Units', 'Normalized', ... + 'Enable','Off', ... + 'Position', [myleft 0.549 mywid 0.075], ... + 'String', 'Update Ensemble' , ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.4, ... + 'Callback', @button_update_ens_Callback); + +%% ----------------------------------------------------------------------------- + +handles.inflation_panel = uipanel('BackgroundColor',atts.lightblue, ... + 'BorderType','none', ... + 'Position',[0.661 0.282 0.32 0.248]); + +handles.ui_checkbox_inflation = uicontrol(handles.inflation_panel, ... + 'Style', 'checkbox', ... + 'Units', 'Normalized', ... + 'Position', [0.037 0.624 0.900 0.335], ... + 'String', 'Apply Inflation', ... + 'BackgroundColor', atts.lightblue,... + 'ForegroundColor', 'k', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'Callback', @inflation_toggle_Callback); + +handles.ui_slider_inflation = uicontrol(handles.inflation_panel, ... + 'Style', 'slider', ... + 'Units', 'Normalized', ... + 'Position', [0.058 0.411 0.893 0.18], ... + 'Value', 1,... + 'Max', 5,... + 'Min', 1,... + 'Sliderstep',[0.05 0.2], ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.5, ... + 'Enable', 'Off', ... + 'Callback', @slider_Callback); + +handles.ui_text_inflation = uicontrol(handles.inflation_panel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.018 0.006 0.671 0.29], ... + 'String', 'Inflation Amount', ... + 'BackgroundColor', atts.lightblue,... + 'ForegroundColor', 'k', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'Enable', 'Off', ... + 'FontWeight','Bold'); + +handles.ui_edit_inflation_label = uicontrol(handles.inflation_panel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', ... + 'Position', [0.718 0.052 0.263 0.258], ... + 'String', get(handles.ui_slider_inflation,'Value'), ... + 'BackgroundColor', 'White', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.5, ... + 'Enable', 'Off', ... + 'Callback', @edit_inflation_Callback); +handles.inflation = str2double(get(handles.ui_edit_inflation_label,'String')); + +%% ----------------------------------------------------------------------------- +% doesn't seem to do a great job of centering ... but the distribute is nice. + +hlist = [handles.observation_panel, ... + handles.ui_button_create_new_ens, ... + handles.ui_button_update_ens, ... + handles.inflation_panel]; + +align(hlist,'Center','Distribute') + +%% ----------------------------------------------------------------------------- + +handles.ui_radio_button_group = uibuttongroup ('BackgroundColor', atts.background, ... + 'BorderType', 'none', ... + 'Position',[0.776 0.013 0.200 0.256]); + +handles.ui_radio_button_eakf = uicontrol(handles.ui_radio_button_group, ... + 'Style', 'radio button', ... + 'Units', 'Normalized', ... + 'Position', [0.007 0.667 670 0.253], ... + 'String', 'EAKF', ... + 'BackgroundColor', atts.background, ... + 'Foreground', 'Black', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontWeight','Bold', ... + 'FontSize', 0.5, ... + 'HandleVisibility', 'On', ... + 'Callback', @all_buttons_Callback); + +handles.ui_radio_button_gamma = uicontrol(handles.ui_radio_button_group, ... + 'Style', 'radio button', ... + 'Units', 'Normalized', ... + 'Position', [0.007 0.347 670 0.253], ... + 'String', 'Gamma', ... + 'BackgroundColor', atts.background, ... + 'Foreground', 'Black', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontWeight','Bold', ... + 'FontSize', 0.5, ... + 'HandleVisibility', 'Off', ... + 'Callback', @all_buttons_Callback); + +handles.ui_radio_button_rhf = uicontrol(handles.ui_radio_button_group, ... + 'Style', 'radio button', ... + 'Units', 'Normalized', ... + 'Position', [0.007 0.06 670 0.253], ... + 'String', 'BNRHF', ... + 'BackgroundColor', atts.background, ... + 'Foreground', 'Black', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontWeight','Bold', ... + 'FontSize', 0.5, ... + 'HandleVisibility', 'Off', ... + 'Callback', @all_buttons_Callback); + +%% ----------------------------------------------------------------------------- + +handles.axes = axes ('Units', 'Normalized', ... + 'Position', [30/figWidth 30/figWidth 0.6000 0.9000], ... + 'FontName', atts.fontname, ... + 'FontSize', atts.fontsize, ... + 'Color','White'); + +% This section specifies the annotation for the values on the main graph. +% By specifying them this way, we can turn them On/Off at will and specify +% that they scale when the object is resized. +% FIXME ... these should be part of some uipanel that is on the graphic. + +handles.ui_text_prior_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.852 0.233 0.065], ... + 'String', 'Prior Mean = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'right'); + +handles.ui_text_inflated_prior_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.806 0.233 0.065], ... + 'String', 'Inflated = ', ... + 'Visible', 'Off', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_prior_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.762 0.233 0.065], ... + 'String', 'Prior SD = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'right'); + +handles.ui_text_inflated_prior_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.713 0.233 0.065], ... + 'String', 'Inflated = ', ... + 'Visible', 'Off', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor',atts.green,... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_post_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.852 0.270 0.065], ... + 'String', 'Posterior Mean = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_inflated_post_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.806 0.270 0.065], ... + 'String', 'Inflated = ', ... + 'Visible', 'Off', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_post_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.762 0.270 0.065], ... + 'String', 'Posterior SD = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_inflated_post_sd = uicontrol('Style', 'text', ... + 'Units' , 'Normalized', ... + 'Position' , [0.340 0.713 0.270 0.065], ... + 'String' , 'Inflated = ', ... + 'Visible' , 'Off', ... + 'BackgroundColor' , 'White', ... + 'ForegroundColor' , atts.blue, ... + 'FontUnits' , 'normalized', ... + 'FontName' , atts.fontname, ... + 'FontSize' , 0.4, ... + 'FontWeight' , 'Bold', ... + 'HorizontalAlignment', 'right'); + +% justify all the strings: + +set(handles.ui_text_prior_mean , 'String', 'Prior Mean = '); +set(handles.ui_text_prior_sd , 'String', 'Prior SD = '); +set(handles.ui_text_inflated_prior_mean , 'String', 'Inflated = '); +set(handles.ui_text_inflated_prior_sd , 'String', 'Inflated = '); + +set(handles.ui_text_post_mean , 'String', 'Posterior Mean = '); +set(handles.ui_text_post_sd , 'String', 'Posterior SD = '); +set(handles.ui_text_inflated_post_mean , 'String', 'Inflated = '); +set(handles.ui_text_inflated_post_sd , 'String', 'Inflated = '); + +% FIXME These are error messages that may or may not be present. +% We create them and then turn them off. When we need to, IF we need to, +% they can be turned back on. + +handles.ui_text_inf_err_print = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [40/figWidth 260/figHeight 400/figWidth 40/figHeight], ... + 'String', 'ERROR: Inflation value must be between 1 and 5.', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.red, ... + 'FontSize', atts.fontsize, ... + 'FontWeight', 'Bold', ... + 'FontName', atts.fontname, ... + 'Visible', 'Off'); + +handles.ui_text_obs_err_print = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [40/figWidth 260/figHeight 400/figWidth 40/figHeight], ... + 'String', 'ERROR: Observation value must be numeric.', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.red, ... + 'FontSize', atts.fontsize, ... + 'FontWeight', 'Bold', ... + 'FontName', atts.fontname, ... + 'Visible', 'Off'); + +handles.ui_text_obs_sd_err_print = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [40/figWidth 260/figHeight 400/figWidth 40/figHeight], ... + 'String', 'ERROR: Obs. Error SD value must be numeric.', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.red, ... + 'FontSize', atts.fontsize, ... + 'FontWeight', 'Bold', ... + 'FontName', atts.fontname, ... + 'Visible','Off'); + +% Go ahead and plot the initial observational error distribution + +handles.h_obs_plot = plot_gaussian(handles.observation, handles.obs_error_sd, 1); +set(handles.h_obs_plot, 'Color', atts.red, 'Linestyle', '--', 'Linewidth', 2.0); +hold on + +% 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 +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, 'FontSize', atts.fontsize) +set(gca, 'XGrid', 'on') +set(gca, 'color', [0.8, 0.8, 0.8]); + +% Plot an axis line +plot([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); +% Plot a vertical boundary at 0 since this is bounded +plot([0 0], [0 0.7], 'k--', 'Linewidth', 3.0) +set(gca, 'FontUnits', 'Normalized'); +title('bounded_oned_ensemble','Interpreter','none') + +% Format the information and error messages for ensemble creation +h_click = text(0, 0.6, {'Click inside grey graphics box to create member', ... + '(only X value is used)'}, 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center', ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_err_text = text(0, -0.15, 'An ensemble has to have at least 2 members.', ... + 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center','Color', atts.red, ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_negative_text = text(0, -0.25, 'Ensemble members must be nonnegative.', ... + 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center','Color', atts.red, ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_finish = text(0, -0.15, 'Click outside of grey to finish', ... + 'Fontsize', atts.fontsize, 'HorizontalAlignment', 'center', ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); + + +%% ----------------------------------------------------------------------------- + + function button_create_new_ens_Callback(~,~) + + % Disable the update ensemble button and all other active buttons + set(handles.ui_button_update_ens, 'Enable', 'Off'); + set(handles.ui_edit_observation, 'Enable', 'Off'); + set(handles.ui_edit_obs_error_sd, 'Enable', 'Off'); + set(handles.ui_edit_inflation_label, 'Enable', 'Off'); + + % Clear out any old ensemble members if they exist + set(handles.h_ens_member, 'Visible', 'Off'); + set(handles.h_inf_ens_member, 'Visible', 'Off'); + + set(handles.h_update_lines, 'Visible', 'Off'); + set(handles.h_inf_lines, 'Visible', 'Off'); + set(handles.h_inf_axis, 'Visible', 'Off'); + + % Turn Off any old update points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_inf_up_ens, 'Visible', 'Off'); + set(handles.h_inf_ens_member, 'Visible', 'Off'); + + clear_ui_labels; + + hold on + + % Set a basic plotting domain range that includes mean +/- 3 obs SDs + 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]); + + % Plot the observation distribution + update_obs_dist_plot(); + + % Messages are centered in the middle. + xmid = (xupper + xlower) / 2.0; + set(h_click, 'Position', [xmid, 0.6, 0], 'Visible', 'On'); + set(h_err_text, 'Position', [xmid, -0.15, 0], 'Visible', 'On'); + set(h_negative_text, 'Position', [xmid, -0.25, 0], 'Visible', 'Off'); + set(h_finish, 'Position', [xmid, -0.15, 0], 'Visible', 'Off'); + + ens_size = 0; + + while ens_size < 100 + [xt, yt] = ginput(1); + + % Bounded below at 0, ensemble members must be nonnegative + if(xt >= 0.0 && xt <= xupper && yt >= ylower && yt <= yupper) + ens_size = ens_size + 1; + x(ens_size) = xt; %#ok + y(ens_size) = 0; %#ok + handles.h_ens_member(ens_size) = ... + plot(x(ens_size), y(ens_size), '*', 'MarkerSize', 16, 'Color', atts.green,'LineWidth',2.0); + + % Display the prior mean and sd + prior_mean = mean(x); + prior_sd = std(x); + str1 = sprintf('Prior Mean = %.4f', prior_mean); + set(handles.ui_text_prior_mean, 'String', str1); + str1 = sprintf('Prior SD = %.4f', prior_sd); + set(handles.ui_text_prior_sd, 'String', str1); + % Turn off the warning about no negative members + set(h_negative_text, 'Visible', 'Off'); + + elseif(xt >= xlower && xt < 0.0 && yt >= ylower && yt <= yupper) + % This is a negative x value so not allowed + set(h_negative_text, 'Visible', 'On'); + + elseif (ens_size < 2) + set(h_err_text,'FontWeight','bold') + + else + % This is the case to exit with 2 or more members; turn off the bounded message + set(h_negative_text, 'Visible', 'Off'); + break; + + 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 + + % Ensemble created, compute mean and sd, clean up and return + % Set the global gui storage + handles.ens_size = ens_size; + handles.ens_members = x; + + % Turn Off the data entry messages + set(h_click, 'Visible', 'Off'); + set(h_finish, 'Visible', 'Off'); + + % Enable the update ensemble button + set(handles.ui_button_update_ens, 'Enable', 'On'); + set(handles.ui_edit_observation, 'Enable', 'On'); + set(handles.ui_edit_obs_error_sd, 'Enable', 'On'); + set(handles.ui_edit_inflation_label, 'Enable', 'On'); + + end + + +%% ----------------------------------------------------------------------------- + + function inflation_toggle_Callback (~, ~) + + enabled = get(handles.ui_checkbox_inflation, 'Value'); + if (enabled) + set(handles.ui_slider_inflation, 'Enable', 'On'); + set(handles.ui_text_inflation, 'Enable', 'On'); + set(handles.ui_edit_inflation_label, 'Enable', 'On'); + set(handles.ui_text_inflated_prior_mean, 'Visible', 'On'); + set(handles.ui_text_inflated_post_mean, 'Visible', 'On'); + set(handles.ui_text_inflated_prior_sd, 'Visible', 'On'); + set(handles.ui_text_inflated_post_sd, 'Visible', 'On'); + + else + set(handles.ui_slider_inflation, 'Enable', 'Off'); + set(handles.ui_text_inflation, 'Enable', 'Off'); + set(handles.ui_edit_inflation_label, 'Enable', 'Off'); + set(handles.ui_text_inflated_prior_mean, 'Visible', 'Off'); + set(handles.ui_text_inflated_post_mean, 'Visible', 'Off'); + set(handles.ui_text_inflated_prior_sd, 'Visible', 'Off'); + set(handles.ui_text_inflated_post_sd, 'Visible', 'Off'); + end + + + % Plot the updated distribution + update_obs_dist_plot(); + + end + +%% ----------------------------------------------------------------------------- + + function all_buttons_Callback(~, ~) + + % Plot the updated distribution + update_obs_dist_plot(); + + end + +%% ----------------------------------------------------------------------------- + + function slider_Callback (~, ~) + + handles.inflation = get(handles.ui_slider_inflation, 'Value'); + + str1 = sprintf('%.4f',handles.inflation); + set(handles.ui_edit_inflation_label, 'String', str1); + + % Just in case the inflation label was in the error state, reset + set(handles.ui_edit_inflation_label, 'BackgroundColor', 'White', 'FontWeight', 'Normal'); + set(handles.ui_text_inf_err_print, 'Visible', 'Off') + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_observation, 'Enable', 'On') + set(handles.ui_edit_obs_error_sd, 'Enable', 'On') + set(handles.ui_button_create_new_ens, 'Enable', 'On') + set(handles.ui_button_update_ens, 'Enable', 'On') + + + % Plot the updated distribution + update_obs_dist_plot(); + + end + +%% ----------------------------------------------------------------------------- + + function button_update_ens_Callback (~, ~) + + % Turn Off any old points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_inf_up_ens, 'Visible', 'Off'); + set(handles.h_inf_ens_member, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + set(handles.h_inf_lines, 'Visible', 'Off'); + set(handles.h_inf_axis, 'Visible', 'Off'); + + ensemble = handles.ens_members; + + % Plot the updated distribution + update_obs_dist_plot(); + + % Figure out which filter option is currently selected + val = get(handles.ui_radio_button_group,'SelectedObject'); + filter_type = get(val,'String'); + + switch filter_type + + case 'EAKF' + % Plot the continuous prior distribution + handles.h_prior_pdf = plot_gaussian(mean(ensemble), std(ensemble), 1); + set(handles.h_prior_pdf, 'Color', atts.green, 'Linewidth', 2.0); + + [obs_increments, ~] = ... + obs_increment_eakf(ensemble, handles.observation, handles.obs_error_sd^2); + + % Add on increments to get new ensemble + new_ensemble = ensemble + obs_increments; + + % Plot the continuous posterior distribution + handles.h_post_pdf = plot_gaussian(mean(new_ensemble), std(new_ensemble), 1); + set(handles.h_post_pdf, 'Color', atts.blue, 'Linewidth', 2.0); + + case 'Gamma' + % Get the shape and scale of the likelihood + lgamma_mean = handles.observation; + lgamma_sd = handles.obs_error_sd; + lgamma_var = lgamma_sd .^ 2; + lgamma_shape = lgamma_mean.^2 / lgamma_var; + lgamma_scale = lgamma_var / lgamma_mean; + + % Get the shape and scale of the prior + pgamma_params = gamfit(ensemble); + pgamma_shape = pgamma_params(1); pgamma_scale = pgamma_params(2); + % Plot the prior + handles.h_prior_pdf = plot_gamma(pgamma_shape, pgamma_scale); + set(handles.h_prior_pdf, 'Color', atts.green, 'Linewidth', 2.0); + + % Compute the quantiles of the prior ensemble members + prior_q = gamcdf(ensemble, pgamma_shape, pgamma_scale); + + % Compute the posterior + agamma_shape = pgamma_shape + lgamma_shape - 1; + agamma_scale = pgamma_scale * lgamma_scale / (pgamma_scale + lgamma_scale); + + % Plot the posterior (analysis) + handles.h_post_pdf = plot_gamma(agamma_shape, agamma_scale); + set(handles.h_post_pdf, 'Color', atts.blue, 'Linewidth', 2.0); + + % Get the posterior ensemble + new_ensemble = gaminv(prior_q, agamma_shape, agamma_scale); + case 'BNRHF' + bounded_left = true; + [obs_increments, err, xp_prior, yp_prior, xp_post, yp_post] = ... + obs_increment_rhf(ensemble, handles.observation, handles.obs_error_sd^2, bounded_left); + handles.h_prior_pdf = plot(xp_prior, yp_prior, 'linewidth', 2, 'color', atts.green); + handles.h_post_pdf = plot(xp_post, yp_post, 'linewidth', 2, 'color', atts.blue); + + % Add on increments to get new ensemble + new_ensemble = ensemble + obs_increments; + end + + + y(1:size(ensemble, 2)) = -0.1; + handles.h_update_ens = plot(new_ensemble, y, '*', 'MarkerSize', 16, 'Color', atts.blue); + + % 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'); + end + + % Add in a label of the updated mean and sd + new_mean = mean(new_ensemble); + new_sd = std(new_ensemble); + + % Update mean and sd of old posterior + str1 = sprintf('Posterior Mean = %.4f',new_mean); + set(handles.ui_text_post_mean, 'String', str1, 'Visible', 'on'); + + str1 = sprintf('Posterior SD = %.4f',new_sd); + set(handles.ui_text_post_sd, 'String', str1, 'Visible', 'on'); + + % If the checkbox isn't set, return now + if(not(get(handles.ui_checkbox_inflation, 'Value'))) + return + end + + % Inflate the prior ensemble + inf_ens = zeros(1,handles.ens_size); + handles.prior_mean = mean(handles.ens_members(1:handles.ens_size)); + switch filter_type + + case 'EAKF' + for i = 1: handles.ens_size + inf_ens(i) = (handles.ens_members(i) - handles.prior_mean) ... + * sqrt(handles.inflation) + handles.prior_mean; + end + + case 'Gamma' + inf_ens = inflate_gamma(handles.ens_members, handles.ens_size, handles.inflation); + + case 'BNRHF' + inf_ens = inflate_bnrh(handles.ens_members, handles.ens_size, handles.inflation, ... + true, false, 0, -99); + end + + % Plot the inflated ensemble members + y = -0.2; + for i = 1:handles.ens_size + handles.h_inf_ens_member(i) = plot(inf_ens(i), y, '*', 'MarkerSize', 16, 'Color', atts.green,'LineWidth',2.0); + end + + % Update mean and sd of old posterior + handles.inf_prior_mean = mean(inf_ens(1:handles.ens_size)); + handles.inf_prior_sd = std(inf_ens(1:handles.ens_size)); + + str1 = sprintf('Inflated = %.4f',handles.inf_prior_mean); + set(handles.ui_text_inflated_prior_mean,'String',str1,'Visible','on'); + str1 = sprintf('Inflated = %.4f',handles.inf_prior_sd); + set(handles.ui_text_inflated_prior_sd, 'String',str1,'Visible','on'); + + % Get the update for the inflated ensemble + switch filter_type + + case 'EAKF' + [obs_increments, ~] = ... + obs_increment_eakf(inf_ens, handles.observation, handles.obs_error_sd^2); + % Add on increments to get new ensemble + new_ensemble = inf_ens + obs_increments; + + case 'Gamma' + % Get the shape and scale of the inflated prior + pgamma_params = gamfit(inf_ens); + pgamma_shape = pgamma_params(1); pgamma_scale = pgamma_params(2); + + % Compute the quantiles of the prior ensemble members + prior_q = gamcdf(inf_ens, pgamma_shape, pgamma_scale); + + % Compute the posterior: Likelihood is still the same + agamma_shape = pgamma_shape + lgamma_shape - 1; + agamma_scale = pgamma_scale * lgamma_scale / (pgamma_scale + lgamma_scale); + + % Get the posterior ensemble + new_ensemble = gaminv(prior_q, agamma_shape, agamma_scale); + + case 'BNRHF' + [obs_increments, ~] = ... + obs_increment_rhf(inf_ens, handles.observation, handles.obs_error_sd^2, bounded_left); + % Add on increments to get new ensemble + new_ensemble = inf_ens + obs_increments; + end + + y(1:size(ensemble, 2)) = -0.3; + handles.h_inf_up_ens = plot(new_ensemble, y, '*', 'MarkerSize', 16, 'Color', atts.blue); + + % 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'); + + end + + % Set a basic plotting domain range that includes mean +/- 3 obs SDs + % Plus all inflated members + 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([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); + handles.h_inf_axis = plot([xlower xupper], [-0.2 -0.2], 'k', 'Linewidth', 1.7); + + % Update mean and sd of old posterior + handles.update_inf_mean = mean(new_ensemble(1:handles.ens_size)); + handles.update_inf_sd = std (new_ensemble(1:handles.ens_size)); + + str1 = sprintf('Inflated = %.4f',handles.update_inf_mean); + set(handles.ui_text_inflated_post_mean, 'String', str1, 'Visible','on'); + + str1 = sprintf('Inflated = %.4f',handles.update_inf_sd); + set(handles.ui_text_inflated_post_sd, 'String', str1, 'Visible', 'on'); + + end + +%% ----------------------------------------------------------------------------- + + function clear_ui_labels() + + % Turns Off all labels except for the prior mean and SD + set(handles.ui_text_post_sd, 'Visible', 'Off'); + set(handles.ui_text_post_mean, 'Visible', 'Off'); + set(handles.ui_text_inflated_prior_mean, 'Visible', 'Off'); + set(handles.ui_text_inflated_prior_sd, 'Visible', 'Off'); + set(handles.ui_text_inflated_post_sd, 'Visible', 'Off'); + set(handles.ui_text_inflated_post_mean, 'Visible', 'Off'); + + end + +%% ----------------------------------------------------------------------------- + + function edit_inflation_Callback(~, ~) + + % Turn Off any old updated points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_inf_up_ens, 'Visible', 'Off'); + set(handles.h_inf_ens_member, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + set(handles.h_inf_lines, 'Visible', 'Off'); + set(handles.h_inf_axis, 'Visible', 'Off'); + + % Enable things that an error might have turned Off + set(handles.ui_edit_observation, 'Enable', 'on') + set(handles.ui_edit_obs_error_sd, 'Enable', 'on') + set(handles.ui_button_create_new_ens, 'Enable', 'on') + + % Only enable the update ensemble pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_update_ens, 'Enable', 'on'); + + end + + % Get the value of the inflation + inf_value = str2double(get(handles.ui_edit_inflation_label, 'String')); + + if( isfinite(inf_value) && (inf_value >= 1) && (inf_value <= 5)) + inflation = inf_value; + + else + set(handles.ui_edit_inflation_label, 'String', '??','FontWeight','Bold', ... + 'BackgroundColor', atts.red); + set(handles.ui_text_inf_err_print,'Visible','On') + + fprintf('ERROR: Inflation value must be between 1 and 5.\n') + fprintf('ERROR: Inflation value must be between 1 and 5.\n') + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_observation, 'Enable', 'Off') + set(handles.ui_edit_obs_error_sd, 'Enable', 'Off') + set(handles.ui_button_create_new_ens, 'Enable', 'Off') + set(handles.ui_button_update_ens, 'Enable', 'Off') + + return + + end + + % Update the value in global storage + handles.inflation = inflation; + set(handles.ui_edit_inflation_label, 'BackgroundColor', 'White', 'FontWeight', 'Normal'); + set(handles.ui_slider_inflation,'Value', handles.inflation); + set(handles.ui_text_inf_err_print,'Visible','Off') + + + % Plot the updated distribution + update_obs_dist_plot() + + % Set a basic plotting domain range that includes mean +/- 3 obs SDs + 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]); + + hold on + + plot([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); + + end + +%% ----------------------------------------------------------------------------- + + function edit_observation_Callback(~, ~) + + % Turn Off any old updated points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_inf_up_ens, 'Visible', 'Off'); + set(handles.h_inf_ens_member, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + set(handles.h_inf_lines, 'Visible', 'Off'); + set(handles.h_inf_axis, 'Visible', 'Off'); + + % Enable things that an error might have turned Off + set(handles.ui_edit_obs_error_sd, 'Enable', 'on') + set(handles.ui_edit_inflation_label, 'Enable', 'on') + set(handles.ui_button_create_new_ens, 'Enable', 'on') + + % Only enable the update ensemble pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_update_ens, 'Enable', 'on'); + + end + + % Get the value of the observation + if(isfinite( str2double(get(handles.ui_edit_observation, 'String')))) + observation = str2double(get(handles.ui_edit_observation, 'String')); + + else + set(handles.ui_edit_observation, 'String', '??','FontWeight','Bold', ... + 'BackgroundColor', atts.red); + set(handles.ui_text_obs_err_print,'Visible','On') + + fprintf('ERROR: Observation value must be numeric.\n') + fprintf('ERROR: Observation value must be numeric.\n') + + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_obs_error_sd, 'Enable', 'Off') + set(handles.ui_edit_inflation_label, 'Enable', 'Off') + set(handles.ui_button_create_new_ens, 'Enable', 'Off') + set(handles.ui_button_update_ens, 'Enable', 'Off') + + return + + end + + % Update the global storage + handles.observation = observation; + set(handles.ui_edit_observation, 'BackgroundColor', 'White','FontWeight', 'Normal'); + set(handles.ui_text_obs_err_print,'Visible','Off') + + % Plot the updated distribution + update_obs_dist_plot(); + + % Move the observation asterisk + set(handles.h_obs_ast, 'Visible', 'Off'); + 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 + 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]); + + hold on + plot([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); + + end + +%% ----------------------------------------------------------------------------- + + function edit_obs_error_sd_Callback(~, ~) + + % Turn Off any old updated points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_inf_up_ens, 'Visible', 'Off'); + set(handles.h_inf_ens_member, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + set(handles.h_inf_lines, 'Visible', 'Off'); + set(handles.h_inf_axis, 'Visible', 'Off'); + + % Enable things that an error might have turned Off + set(handles.ui_edit_observation, 'Enable', 'on') + set(handles.ui_edit_inflation_label, 'Enable', 'on') + set(handles.ui_button_create_new_ens, 'Enable', 'on') + + % Only enable the update ensemble pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_update_ens, 'Enable', 'on'); + end + + % Get the value of the observation error sd + obs_error_value = str2double(get(handles.ui_edit_obs_error_sd, 'String')); + + if(isfinite(obs_error_value) && (obs_error_value > 0)) + obs_error_sd = obs_error_value; + + else + + set(handles.ui_edit_obs_error_sd, 'String', '??','FontWeight','Bold', ... + 'BackgroundColor', atts.red); + set(handles.ui_text_obs_sd_err_print,'Visible','On') + + fprintf('ERROR: Obs. Error SD value must be numeric.\n') + fprintf('ERROR: Obs. Error SD value must be numeric.\n') + + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_observation, 'Enable', 'Off') + set(handles.ui_edit_inflation_label, 'Enable', 'Off') + set(handles.ui_button_create_new_ens, 'Enable', 'Off') + set(handles.ui_button_update_ens, 'Enable', 'Off') + + return + + end + + % Update the value in global storage + handles.obs_error_sd = obs_error_sd; + set(handles.ui_edit_obs_error_sd, 'BackgroundColor', 'White', 'FontWeight', 'Normal'); + set(handles.ui_text_obs_sd_err_print,'Visible','Off') + + % Plot the updated distribution + update_obs_dist_plot(); + + % Set a basic plotting domain range that includes mean +/- 3 obs SDs + 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]); + + hold on + + plot([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); + + end + +%% ----------------------------------------------------------------------------- + + function update_obs_dist_plot + + % Plot the observation distribution + % Figure out which filter option is currently selected + val = get(handles.ui_radio_button_group,'SelectedObject'); + filter_type = get(val,'String'); + + % Turn the old obsevation likelihood plot off + set(handles.h_obs_plot, 'Visible', 'Off'); + + % Gamma prior and likelihood + if(strcmp(filter_type, 'Gamma')) + gmean = handles.observation; + gsd = handles.obs_error_sd; + gvar = gsd .^ 2; + gamma_shape = gmean.^2 / gvar; + gamma_scale = gvar / gmean; + + handles.h_obs_plot = plot_gamma(gamma_shape, gamma_scale); + else + handles.h_obs_plot = plot_gaussian(handles.observation, handles.obs_error_sd, 1); + end + set(handles.h_obs_plot, 'Color', atts.red, 'Linestyle', '--', 'Linewidth', 1.7); + + end + +%% ----------------------------------------------------------------------------- + +end + +% +% $URL$ +% $Revision$ +% $Date$ diff --git a/guide/DART_LAB/matlab/gaussian_product.m b/guide/DART_LAB/matlab/gaussian_product.m index 3f8c4c048c..66f67c9bf1 100644 --- a/guide/DART_LAB/matlab/gaussian_product.m +++ b/guide/DART_LAB/matlab/gaussian_product.m @@ -3,15 +3,16 @@ % % 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) +% Gaussian for the Prior (green) and the Observation (red) % and click on 'Plot Posterior'. % -% 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. +% 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.m oned_model_inf.m oned_ensemble.m -% twod_ensemble.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% See also: bounded_oned_ensemble.m oned_cycle.m oned_ensemble.m oned_model.m +% oned_model_inf.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at @@ -19,6 +20,8 @@ % % DART $Id$ +help gaussian_product + atts = stylesheet; % get the default fonts and colors % Specify the figure size in pixels. After that, all positions are @@ -248,6 +251,7 @@ %%Plots an initial graph with the default values g_prod_plot(handles); +set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); %% --------------------------------------------------------------------- @@ -257,6 +261,7 @@ function plotGraph_Callback(~,~) %handles. This is done through the function definition. [prior_mean, prior_sd, obs_mean, obs_err_sd, is_err] = g_prod_plot(handles); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); % If there is an error, zero out the posterior text values % don't try to do posterior computation @@ -300,21 +305,25 @@ function plotGraph_Callback(~,~) function edit_prior_mean_Callback(~, ~) g_prod_plot(handles); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); reset_Posterior(); end function edit_prior_sd_Callback(~, ~) g_prod_plot(handles); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); reset_Posterior(); end function edit_observation_Callback(~, ~) g_prod_plot(handles); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); reset_Posterior(); end function edit_obs_error_sd_Callback(~, ~) g_prod_plot(handles); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); reset_Posterior(); end diff --git a/guide/DART_LAB/matlab/oned_cycle.m b/guide/DART_LAB/matlab/oned_cycle.m new file mode 100644 index 0000000000..aa596d240d --- /dev/null +++ b/guide/DART_LAB/matlab/oned_cycle.m @@ -0,0 +1,843 @@ +function oned_cycle +%% ONED_CYCLE explores the cycling Kalman Filter for a single variable and allows +% comparison of the continuous KF to three ensemble filter algorithms. +% +% The system being explored has a true value that is always zero (marked by +% the black 'X' on the plot axis). The observation error standard deviation (SD) +% is specified by the 'Obs. Error SD' entry in the red box. A single observation +% is generated for each assimilation cycle as a random draw from Normal(0, SD). +% +% The forecast model is a linear error growth with growth rate (G) specified by the +% 'Growth Rate' entry in the red box. The difference equation is +% x(t+1) = Gx(t) where each integer value of t is the time of an assimilation. +% +% The model state is a Gaussian distribution. The initial distribution has mean +% 1 and standard deviation 2. Pressing the 'Cycle DA' button +% advances the model state to get the prior distribution (green) at the next +% assimilation time and computes a Kalman filter update at that time to get a +% posterior distribution (blue). The observation value at this time is plotted +% as a red asterisk and the likelihood as a red curve. The prior and +% posterior mean and standard deviation are plotted below the axis. +% +% With the growth rate set to the default value of 1, cycling is equivalent +% to assimilation a sequence of observations at a single time. +% +% The continuous Kalman filter can also be compared to three types of ensemble +% filters assimilating the same observations. The Ensemble Adjustment +% Kalman filter (EAKF), perturbed obs ensemble Kalman filter (EnKF), and the +% rank histogram filter (RHF) can be selected with the pushbuttons at the +% bottom right. +% +% Push on the 'Create New Ensemble' button to activate the interactive +% mechanism to generate a prior ensemble estimate. +% This is done by 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. As you add ensemble members, +% the ensemble prior mean and prior standard deviation are updated in the upper +% left corner of the plot. +% +% Once an ensemble is created, the 'Cycle DA' button will also advance the +% ensemble to the next time and perform the selected ensemble assimilation. +% The prior and posterior distributions that are fit to the ensemble are +% plotted by dashed curves for the EAKF and RHF. Prior and posterior normal +% fits to the EnKF ensemble are also plotted. The ensemble prior and posterior +% sample mean and standard deviation are plotted adjacent to their Kalman +% filter counterparts. +% +% See also: bounded_oned_ensemble.m gaussian_product.m oned_ensemble.m oned_model.m +% oned_model_inf.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +help oned_cycle + +atts = stylesheet; % get the default fonts and colors + +% Set random number seed to same value to generate known sequences +rng('default') + +% Initialize the basic structure we will be using to hold everything. +handles.ens_size = 0; +handles.ens_members = 0; +handles.h_obs_plot = []; +handles.h_update_ens = []; +handles.h_prior_pdf = []; +handles.h_post_pdf = []; +handles.h_prior_cont_pdf = []; +handles.h_post_cont_pdf = []; +handles.h_ens_member = []; +handles.h_obs_ast = []; +handles.h_update_lines = []; + +handles.prior_cont_mean = 1; +handles.prior_cont_sd = 2; + +%% ----------------------------------------------------------------------------- + +% Specify the figure size in pixels. After that, all positions are +% specified as fractions (units=Normalized). That way, the objects +% scale proportionally as the figure gets resized. +figXmin = 450; % The horizontal position of the entire figure, in pixels +figYmin = 250; % The vertical position of the entire figure, in pixels +figWidth = 670; % The width of the entire figure, in pixels +figHeight = 500; % The height of the entire figure, in pixels + +handles.figure1 = figure('Position', [figXmin figYmin figWidth figHeight], ... + 'Units', 'Pixels', ... + 'Name', 'oned_cycle', ... + 'Color', atts.background); + +%% ----------------------------------------------------------------------------- +%Position has units of normalized and therefore the components must be +%Fractions of figure. Also, the actual text has FontUnits which are +%normalized, so the Font Size is a fraction of the text box/ edit box that +%contains it. Making the font size and box size normalized allows the text +%size and box size to change proportionately if the user resizes the +%figure window. This is true for all text/edit boxes and buttons. + +%% ----------------------------------------------------------------------------- +% Set up a parent container so we can move the one container around instead of +% trying to manipulate the positions of all the components. + +handles.observation_panel = uipanel('BackgroundColor',atts.red, ... + 'BorderType','none', ... + 'Units', 'Normalized', ... + 'Position',[0.66 0.744 0.325 0.2]); + +handles.ui_text_growth_rate = uicontrol(handles.observation_panel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.035 0.600 0.600 0.275], ... + 'String', 'Growth Rate' , ... + 'BackgroundColor', atts.red, ... + 'ForegroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'FontWeight', 'Bold', ... + 'HorizontalAlignment', 'center'); + +handles.ui_edit_growth_rate= uicontrol(handles.observation_panel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', ... + 'Position', [0.675 0.562 0.300 0.373], ... + 'String', '1', ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'Callback', @edit_growth_rate_Callback); +handles.growth_rate = str2double(get(handles.ui_edit_growth_rate,'String')); + +handles.ui_text_obs_error_sd = uicontrol(handles.observation_panel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.035 0.100 0.600 0.275], ... + 'String', 'Obs. Error SD', ... + 'BackgroundColor', atts.red,... + 'ForegroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'center'); + +handles.ui_edit_obs_error_sd = uicontrol(handles.observation_panel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', ... + 'Position', [0.675 0.062 0.300 0.373], ... + 'String', '1', ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.6, ... + 'Callback', @edit_obs_error_sd_Callback); +handles.obs_error_sd = str2double(get(handles.ui_edit_obs_error_sd,'String')); + +% The true observation value can be initialized to 0 for plotting +handles.observation = 0; + +%% ----------------------------------------------------------------------------- +% Try to center the boxes with what is on top. + +box = get(handles.observation_panel,'Position'); +center = box(1) + box(3)/2.0; +mywid = 0.250; +myleft = center - mywid/2.0; + +handles.ui_button_create_new_ens = uicontrol('Style', 'pushbutton', ... + 'Units', 'Normalized', ... + 'Position', [myleft 0.3 mywid 0.075], ... + 'String', 'Create New Ensemble', ... + 'BackgroundColor', 'White', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'Callback', @button_create_new_ens_Callback); + +mywid = 0.200; +myleft = center - mywid/2.0; +handles.ui_button_cycle_da = uicontrol('style', 'pushbutton', ... + 'Units', 'Normalized', ... + 'Enable','On', ... + 'Position', [myleft 0.6 mywid 0.075], ... + 'String', 'Cycle DA' , ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.4, ... + 'Callback', @button_cycle_da_Callback); + +%% ----------------------------------------------------------------------------- + +% doesn't seem to do a great job of centering ... but the distribute is nice. + +hlist = [handles.observation_panel, ... + handles.ui_button_create_new_ens, ... + handles.ui_button_cycle_da]; + +align(hlist,'Center','Distribute') + +%% ----------------------------------------------------------------------------- + +handles.ui_radio_button_group = uibuttongroup ('BackgroundColor', atts.background, ... + 'BorderType', 'none', ... + 'Position',[0.776 0.013 0.200 0.256]); + +handles.ui_radio_button_eakf = uicontrol(handles.ui_radio_button_group, ... + 'Style', 'radio button', ... + 'Units', 'Normalized', ... + 'Position', [0.007 0.667 670 0.253], ... + 'String', 'EAKF', ... + 'BackgroundColor', atts.background, ... + 'Foreground', 'Black', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontWeight','Bold', ... + 'FontSize', 0.5, ... + 'HandleVisibility', 'On'); + +handles.ui_radio_button_enkf = uicontrol(handles.ui_radio_button_group, ... + 'Style', 'radio button', ... + 'Units', 'Normalized', ... + 'Position', [0.007 0.347 670 0.253], ... + 'String', 'EnKF', ... + 'BackgroundColor', atts.background, ... + 'Foreground', 'Black', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontWeight','Bold', ... + 'FontSize', 0.5, ... + 'HandleVisibility', 'Off'); + +handles.ui_radio_button_rhf = uicontrol(handles.ui_radio_button_group, ... + 'Style', 'radio button', ... + 'Units', 'Normalized', ... + 'Position', [0.007 0.06 670 0.253], ... + 'String', 'RHF', ... + 'BackgroundColor', atts.background, ... + 'Foreground', 'Black', ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontWeight','Bold', ... + 'FontSize', 0.5, ... + 'HandleVisibility', 'Off'); + +%% ----------------------------------------------------------------------------- + +handles.axes = axes ('Units', 'Normalized', ... + 'Position', [30/figWidth 30/figWidth 0.6000 0.9000], ... + 'FontName', atts.fontname, ... + 'FontSize', atts.fontsize, ... + 'Color','White'); + +% This section specifies the annotation for the values on the main graph. +% By specifying them this way, we can turn them On/Off at will and specify +% that they scale when the object is resized. +% FIXME ... these should be part of some uipanel that is on the graphic. + +handles.ui_text_prior_ens_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.152 0.273 0.065], ... + 'String', 'Prior Ens Mean = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'right'); + +handles.ui_text_prior_cont_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.122 0.273 0.065], ... + 'String', 'Prior KF Mean = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'right'); + +handles.ui_text_prior_ens_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.062 0.273 0.065], ... + 'String', 'Prior Ens SD = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'right'); + +handles.ui_text_prior_cont_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.052 0.032 0.273 0.065], ... + 'String', 'Prior KF SD = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.green, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment', 'right'); + +handles.ui_text_post_ens_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.152 0.300 0.065], ... + 'String', 'Posterior Ens Mean = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_post_cont_mean = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.122 0.300 0.065], ... + 'String', 'Posterior KF Mean = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_post_ens_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.062 0.300 0.065], ... + 'String', 'Posterior Ens SD = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +handles.ui_text_post_cont_sd = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [0.340 0.032 0.300 0.065], ... + 'String', 'Posterior KF SD = ', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.blue, ... + 'FontUnits', 'normalized', ... + 'FontName', atts.fontname, ... + 'FontSize', 0.4, ... + 'FontWeight','Bold', ... + 'HorizontalAlignment','right'); + +% justify all the strings, but don't turn them on yet: +set(handles.ui_text_prior_ens_mean, 'String', 'Prior Ens Mean = ', 'visible', 'off'); +set(handles.ui_text_prior_ens_sd, 'String', 'Prior Ens SD = ', 'visible', 'off'); + +set(handles.ui_text_post_ens_mean, 'String', 'Posterior Ens Mean = ', 'visible', 'off'); +set(handles.ui_text_post_cont_mean, 'String', 'Posterior KF Mean = ', 'visible', 'off'); +set(handles.ui_text_post_ens_sd, 'String', 'Posterior Ens SD = ', 'visible', 'off'); +set(handles.ui_text_post_cont_sd, 'String', 'Posterior KF SD = ', 'visible', 'off'); + +% Initial values of the continuous prior distribution +str1 = sprintf('Prior KF Mean = %.4f', handles.prior_cont_mean); +set(handles.ui_text_prior_cont_mean, 'String', str1, 'Visible', 'on'); +str1 = sprintf('Prior KF SD = %.4f', handles.prior_cont_sd); +set(handles.ui_text_prior_cont_sd, 'String', str1, 'Visible', 'on'); + +% These are error messages that may or may not be present. +% We create them and then turn them off. When we need to, IF we need to, +% they can be turned back on. + +handles.ui_text_growth_rate_err_print = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [40/figWidth 260/figHeight 400/figWidth 40/figHeight], ... + 'String', 'ERROR: Growth Rate must be numeric.', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.red, ... + 'FontSize', atts.fontsize, ... + 'FontWeight', 'Bold', ... + 'FontName', atts.fontname, ... + 'Visible', 'Off'); + +handles.ui_text_obs_sd_err_print = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [40/figWidth 260/figHeight 400/figWidth 40/figHeight], ... + 'String', 'ERROR: Obs. Error SD value must be numeric.', ... + 'BackgroundColor', 'White', ... + 'ForegroundColor', atts.red, ... + 'FontSize', atts.fontsize, ... + 'FontWeight', 'Bold', ... + 'FontName', atts.fontname, ... + 'Visible','Off'); + +% Set a basic plotting domain range that includes mean +/- 3 obs SDs +xlower = - 3*handles.obs_error_sd; +xupper = + 3*handles.obs_error_sd; +ylower = -0.4; +yupper = 1.0; +axis([xlower xupper ylower yupper]); +hold on + +%set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]); +set(gca, 'FontSize', atts.fontsize) +set(gca, 'XGrid', 'on') +set(gca, 'color', [0.8, 0.8, 0.8]); + +% Add on a long black zero axis +plot([-1000 1000], [0 0], 'k', 'Linewidth', 1.7); + +% Plot a marker at the truth value +plot(0, 0, 'x', 'color', 'k', 'markersize', 16, 'linewidth', 4); + +set(gca, 'FontUnits', 'Normalized'); +title('oned_cycle','Interpreter','none') + +% Format the information and error messages for ensemble creation +h_click = text(0, 0.6, {'Click inside grey graphics box to create member', ... + '(only X value is used)'}, 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center', ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_err_text = text(0, 0.35, 'An ensemble has to have at least 2 members.', ... + 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center','Color', atts.red, ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_finish = text(0, 0.35, 'Click outside of grey to finish', ... + 'Fontsize', atts.fontsize, 'HorizontalAlignment', 'center', ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); + + +%% ----------------------------------------------------------------------------- + + function button_create_new_ens_Callback(~,~) + + % Disable the update ensemble button and all other active buttons + set(handles.ui_button_cycle_da, 'Enable', 'Off'); + set(handles.ui_edit_growth_rate, 'Enable', 'Off'); + set(handles.ui_edit_obs_error_sd, 'Enable', 'Off'); + + % Clear out any old ensemble members if they exist + set(handles.h_ens_member, 'Visible', 'Off'); + + set(handles.h_update_lines, 'Visible', 'Off'); + + % Turn Off any old update points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_prior_cont_pdf, 'Visible', 'Off'); + set(handles.h_post_cont_pdf, 'Visible', 'Off'); + + % Turn off legend + legend('hide'); + + clear_ui_labels; + + hold on + + % Set a basic plotting domain range that includes mean +/- 3 obs SDs + 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]); + + % Messages are centered in the middle. + xmid = (xupper + xlower) / 2.0; + set(h_click, 'Position', [xmid, 0.6, 0], 'Visible', 'On'); + set(h_err_text, 'Position', [xmid, 0.35, 0], 'Visible', 'On'); + set(h_finish, 'Position', [xmid, 0.35, 0], 'Visible', 'Off'); + + ens_size = 0; + + while ens_size < 100 + [xt, yt] = ginput(1); + + if(xt >= xlower && xt <= xupper && yt >= ylower && yt <= yupper) + ens_size = ens_size + 1; + x(ens_size) = xt; %#ok + y(ens_size) = 0; %#ok + handles.h_ens_member(ens_size) = ... + plot(x(ens_size), y(ens_size), '*', 'MarkerSize', 16, 'Color', atts.green,'LineWidth',2.0); + + % Display the prior ens mean and sd + prior_ens_mean = mean(x); + prior_ens_sd = std(x); + str1 = sprintf('Prior Ens Mean = %.4f', prior_ens_mean); + set(handles.ui_text_prior_ens_mean, 'String', str1, 'Visible', 'on'); + str1 = sprintf('Prior Ens SD = %.4f', prior_ens_sd); + set(handles.ui_text_prior_ens_sd, 'String', str1, 'Visible', 'on'); + + elseif (ens_size < 2) + set(h_err_text,'FontWeight','bold') + + else + break; + + 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 + + % Ensemble created, compute mean and sd, clean up and return + % Set the global gui storage + handles.ens_size = ens_size; + handles.ens_members = x; + + % Turn Off the data entry messages + set(h_click, 'Visible', 'Off'); + set(h_finish, 'Visible', 'Off'); + + % Enable the cycle DA button + set(handles.ui_button_cycle_da, 'Enable', 'On'); + set(handles.ui_edit_growth_rate, 'Enable', 'On'); + set(handles.ui_edit_obs_error_sd, 'Enable', 'On'); + + end + + +%% ----------------------------------------------------------------------------- + + function button_cycle_da_Callback (~, ~) + + % Turn Off any old points + set(handles.h_ens_member, 'Visible', 'Off'); + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_prior_cont_pdf, 'Visible', 'Off'); + set(handles.h_post_cont_pdf, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + + % The true value is zero in this example; obs error sd comes from the box + % Generate a random observation + handles.observation = randn * handles.obs_error_sd; + + % Move the observation asterisk + set(handles.h_obs_ast, 'Visible', 'Off'); + handles.h_obs_ast = plot(handles.observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0); + hold on; + + % Do the assimilation for the continuous normal distribution + post_cont_var = 1 / (1 / handles.prior_cont_sd^2 + 1 / handles.obs_error_sd^2); + post_cont_sd = sqrt(post_cont_var); + post_cont_mean = post_cont_var * (handles.prior_cont_mean / handles.prior_cont_sd^2 + ... + handles.observation / handles.obs_error_sd^2); + + % Plot the prior and posterior continuous normal (KF) PDFs + handles.h_prior_cont_pdf = plot_gaussian(handles.prior_cont_mean, handles.prior_cont_sd, 1); + set(handles.h_prior_cont_pdf, 'Color', atts.green, 'Linestyle', '-', 'Linewidth', 1.7, 'Visible', 'on'); + handles.h_leg(2) = handles.h_prior_cont_pdf; + handles.h_post_cont_pdf = plot_gaussian(post_cont_mean, post_cont_sd, 1); + set(handles.h_post_cont_pdf, 'Color', atts.blue, 'Linestyle', '-', 'Linewidth', 1.7); + handles.h_leg(3) = handles.h_post_cont_pdf; + + % Plot the updated distribution + set(handles.h_obs_plot, 'Visible', 'Off'); + handles.h_obs_plot = plot_gaussian(handles.observation, handles.obs_error_sd, 1); + set(handles.h_obs_plot, 'Color', atts.red, 'Linestyle', '-', 'Linewidth', 1.7); + handles.h_leg(1) = handles.h_obs_plot; + + % Display the prior cont mean and sd + str1 = sprintf('Prior KF Mean = %.4f', handles.prior_cont_mean); + set(handles.ui_text_prior_cont_mean, 'String', str1, 'Visible', 'on'); + str1 = sprintf('Prior KF SD = %.4f', handles.prior_cont_sd); + set(handles.ui_text_prior_cont_sd, 'String', str1, 'Visible', 'on'); + + % Update mean and sd of old cont posterior + str1 = sprintf('Posterior KF Mean = %.4f', post_cont_mean); + set(handles.ui_text_post_cont_mean, 'String', str1, 'Visible', 'on'); + str1 = sprintf('Posterior KF SD = %.4f', post_cont_sd); + set(handles.ui_text_post_cont_sd, 'String', str1, 'Visible', 'on'); + + % Update the continuous for the next cycle + handles.prior_cont_mean = post_cont_mean * handles.growth_rate; + handles.prior_cont_sd = post_cont_sd * handles.growth_rate; + + + % Only plot this three entry legend if no ensemble has been created + ens_exists = handles.ens_size > 1; + if(~ens_exists) legend(handles.h_leg, 'Likelihood', 'Prior', 'Posterior'); end + + % If ensemble has been created, update an plot it also + if(ens_exists) + % Get the prior ensmeble members + ensemble = handles.ens_members; + + % Figure out which filter option is currently selected + val = get(handles.ui_radio_button_group,'SelectedObject'); + filter_type = get(val,'String'); + + switch filter_type + + case 'EAKF' + [obs_increments, ~] = ... + obs_increment_eakf(ensemble, handles.observation, handles.obs_error_sd^2); + + handles.h_prior_pdf = plot_gaussian(mean(ensemble), std(ensemble), 1); + set(handles.h_prior_pdf, 'linewidth', 2, 'color', atts.green, 'Linestyle', '--'); + case 'EnKF' + [obs_increments, ~] = ... + obs_increment_enkf(ensemble, handles.observation, handles.obs_error_sd^2); + % There is no prior distribution to plot for the EnKF, it's not a QCEF + % However, here we can plot the fit to the ensemble anyway to show how it compares + handles.h_prior_pdf = plot_gaussian(mean(ensemble), std(ensemble), 1); + set(handles.h_prior_pdf, 'linewidth', 2, 'color', atts.green, 'Linestyle', '--'); + case 'RHF' + bounded_left = false; + [obs_increments, err, xp_prior, yp_prior, xp_post, yp_post] = ... + obs_increment_rhf(ensemble, handles.observation, handles.obs_error_sd^2, bounded_left); + handles.h_prior_pdf = plot(xp_prior, yp_prior, 'linewidth', 2, 'color', atts.green, 'Linestyle', '--'); + handles.h_post_pdf = plot(xp_post, yp_post, 'linewidth', 2, 'color', atts.blue, 'Linestyle', '--'); + end + + % Add on increments to get new ensemble + new_ensemble = ensemble + obs_increments; + + % Plot the prior and posterior ensemble members + ens_size = size(ensemble, 2); + for i = 1:ens_size + handles.h_ens_member(i) = ... + plot(handles.ens_members(i), 0, '*', 'MarkerSize', 16, 'Color', atts.green); + end + y(1:size(ensemble, 2)) = -0.1; + handles.h_update_ens = plot(new_ensemble, y, '*', 'MarkerSize', 16, 'Color', atts.blue); + + % 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'); + end + + % Add in a label of the updated mean and sd + new_mean = mean(new_ensemble); + new_sd = std(new_ensemble); + + % Display the prior ens mean and sd + prior_ens_mean = mean(handles.ens_members); + prior_ens_sd = std(handles.ens_members); + str1 = sprintf('Prior Ens Mean = %.4f', prior_ens_mean); + set(handles.ui_text_prior_ens_mean, 'String', str1, 'Visible', 'on'); + str1 = sprintf('Prior Ens SD = %.4f', prior_ens_sd); + set(handles.ui_text_prior_ens_sd, 'String', str1, 'Visible', 'on'); + + % Update mean and sd of old ens posterior + str1 = sprintf('Posterior Ens Mean = %.4f',new_mean); + set(handles.ui_text_post_ens_mean, 'String', str1, 'Visible', 'on'); + str1 = sprintf('Posterior Ens SD = %.4f',new_sd); + set(handles.ui_text_post_ens_sd, 'String', str1, 'Visible', 'on'); + + % Plot the posterior sample continuous for the EAKF or ENKF + if(strcmp(filter_type, 'EAKF') | strcmp(filter_type, 'EnKF')) + handles.h_post_pdf = plot_gaussian(new_mean, new_sd, 1); + set(handles.h_post_pdf, 'linewidth', 2, 'color', atts.blue, 'Linestyle', '--'); + end + + handles.h_leg(4) = handles.h_prior_pdf; + handles.h_leg(5) = handles.h_post_pdf; + + legend(handles.h_leg, 'Likelihood', 'Cont. Prior', 'Cont. Post', 'Ens. Prior', 'Ens. Post'); + + % Update the ensemble for next cycle + handles.ens_members = new_ensemble * handles.growth_rate; + end + + end + +%% ----------------------------------------------------------------------------- + + function clear_ui_labels() + + % Turns Off all labels except for the prior ens mean and SD + set(handles.ui_text_post_ens_sd, 'Visible', 'Off'); + set(handles.ui_text_post_ens_mean, 'Visible', 'Off'); + + % Turns Off all labels except for the prior ens mean and SD + set(handles.ui_text_post_cont_sd, 'Visible', 'Off'); + set(handles.ui_text_post_cont_mean, 'Visible', 'Off'); + + end + +%% ----------------------------------------------------------------------------- + + + function edit_growth_rate_Callback(~, ~) + + % Turn Off any old updated points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_prior_cont_pdf, 'Visible', 'Off'); + set(handles.h_post_cont_pdf, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + + % Enable things that an error might have turned Off + set(handles.ui_edit_obs_error_sd, 'Enable', 'on') + set(handles.ui_button_create_new_ens, 'Enable', 'on') + + % Only enable the cycle DA pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_cycle_da, 'Enable', 'on'); + + end + + % Get the value of the growth_rate + if(isfinite( str2double(get(handles.ui_edit_growth_rate, 'String')))) + growth_rate = str2double(get(handles.ui_edit_growth_rate, 'String')); + + else + set(handles.ui_edit_growth_rate, 'String', '??','FontWeight','Bold', ... + 'BackgroundColor', atts.red); + set(handles.ui_text_growth_rate_err_print,'Visible','On') + + fprintf('ERROR: Growth Rate value must be numeric.\n') + fprintf('ERROR: Growth Rate value must be numeric.\n') + + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_obs_error_sd, 'Enable', 'Off') + set(handles.ui_button_create_new_ens, 'Enable', 'Off') + set(handles.ui_button_cycle_da, 'Enable', 'Off') + + return + + end + + % Update the global storage + handles.growth_rate = growth_rate; + set(handles.ui_edit_growth_rate, 'BackgroundColor', 'White','FontWeight', 'Normal'); + set(handles.ui_text_growth_rate_err_print,'Visible','Off') + + end + +%% ----------------------------------------------------------------------------- + + function edit_obs_error_sd_Callback(~, ~) + + % Turn Off any old updated points + set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); + set(handles.h_prior_cont_pdf, 'Visible', 'Off'); + set(handles.h_post_cont_pdf, 'Visible', 'Off'); + + % Remove mean and sd of old posterior + clear_ui_labels; + + % And the lines in between + set(handles.h_update_lines, 'Visible', 'Off'); + + % Enable things that an error might have turned Off + set(handles.ui_edit_growth_rate, 'Enable', 'on') + set(handles.ui_button_create_new_ens, 'Enable', 'on') + + % Only enable the update ensemble pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_cycle_da, 'Enable', 'on'); + end + + % Get the value of the observation error sd + obs_error_value = str2double(get(handles.ui_edit_obs_error_sd, 'String')); + + if(isfinite(obs_error_value) && (obs_error_value > 0)) + obs_error_sd = obs_error_value; + + else + + set(handles.ui_edit_obs_error_sd, 'String', '??','FontWeight','Bold', ... + 'BackgroundColor', atts.red); + set(handles.ui_text_obs_sd_err_print,'Visible','On') + + fprintf('ERROR: Obs. Error SD value must be numeric.\n') + fprintf('ERROR: Obs. Error SD value must be numeric.\n') + + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_growth_rate, 'Enable', 'Off') + set(handles.ui_button_create_new_ens, 'Enable', 'Off') + set(handles.ui_button_cycle_da, 'Enable', 'Off') + + return + + end + + % Update the value in global storage + handles.obs_error_sd = obs_error_sd; + set(handles.ui_edit_obs_error_sd, 'BackgroundColor', 'White', 'FontWeight', 'Normal'); + set(handles.ui_text_obs_sd_err_print,'Visible','Off') + + + % Plot the updated distribution + set(handles.h_obs_plot, 'Visible', 'Off'); + handles.h_obs_plot = plot_gaussian(handles.observation, handles.obs_error_sd, 1); + set(handles.h_obs_plot, 'Color', atts.red, 'Linestyle', '--', 'Linewidth', 1.7); + + % Set a basic plotting domain range that includes mean +/- 3 obs SDs + 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', atts.red, 'Linestyle', '--', 'Linewidth', 1.7); + + set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]); + + hold on + + plot([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); + + end + +%% ----------------------------------------------------------------------------- + +end + +% +% $URL$ +% $Revision$ +% $Date$ diff --git a/guide/DART_LAB/matlab/oned_ensemble.m b/guide/DART_LAB/matlab/oned_ensemble.m index 6f2e4daf29..91d45183f5 100644 --- a/guide/DART_LAB/matlab/oned_ensemble.m +++ b/guide/DART_LAB/matlab/oned_ensemble.m @@ -1,36 +1,52 @@ function oned_ensemble -%% ONED_ENSEMBLE explore the details of ensemble data assimilation for a scalar. +%% ONED_ENSEMBLE explores different ensemble filter assimilation algorithms for a scalar. % % Push 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 -% 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 -% have all the ensemble members you want, click in the grey area of -% the window outside of the white axis plot. +% mechanism to generate a prior ensemble estimate of the observation. +% This is done by 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. As you add ensemble members, +% the ensemble prior mean and prior standard deviation are updated in the upper +% left corner of the plot. % -% 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. +% Once you have created a prior ensemble (green asterisks), click the +% 'Update Ensemble' button. With the default settings, this will apply the EAKF +% algorithm to produce a posterior ensemble (blue asterisks) just below the axis. +% 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. +% Two other ensemble filter variants, the EnKF (sometimes referred to as +% the perturbed observations ensemble Kalman filter) and the rank histogram +% filter (RHF) can be selected with the pushbuttons at the lower right. +% Selecting one of these and pressing 'Update Ensemble' will produce +% the posterior ensemble and posterior statistics using the selected +% filter algorithm. % -% Checking the 'Show Inflation' box will also apply inflation to the -% prior before doing the update and will print the mean and standard +% For the EAKF and RHF, the continuous prior distribution (green) +% that is fit to the prior ensemble and the posterior continuous distribution +% (blue) from which the QCEFF algorithm determines the posterior ensemble +% members are also plotted. The EnKF is not a QCEFF filter and does +% not directly make use of a continuous prior and posterior distribution fit. +% +% Checking the 'Apply Inflation' box will also apply inflation to the +% prior ensemble 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 -% axis for the uninflated ensemble. +% axis for the uninflated ensemble. The value of the covariance inflation +% applied can be adjusted using the slider or by typing in a value in the +% 'Inflation Amount' box. % -% The 'EAKF' is a stochastic algorithm so repeated updates can be done -% for the same prior and observation. +% The 'EnKF' is a stochastic algorithm so repeated updates can be done +% for the same prior and observation by repeatedly pressing 'Update Ensemble'. +% The 'EnKF' version implemented here adjusts the mean of the perturbed +% observations so that they are equal to the actual observation. This means +% that the posterior mean is the same for each repeat. % -% change the Observation Error SD, lay down an ensemble pretty far away -% from the observation - have fun with it. +% The mean and standard deviation of the observation likelihood +% can be changed in the red box. % -% See also: gaussian_product.m oned_model.m oned_model_inf.m -% twod_ensemble.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% See also: bounded_oned_ensemble.m gaussian_product.m oned_cycle.m oned_model.m +% oned_model_inf.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at @@ -50,6 +66,8 @@ handles.ens_members = 0; handles.h_obs_plot = []; handles.h_update_ens = []; +handles.h_prior_pdf = []; +handles.h_post_pdf = []; handles.h_ens_member = []; handles.h_obs_ast = []; handles.h_update_lines = []; @@ -467,10 +485,24 @@ set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]); set(gca, 'FontSize', atts.fontsize) set(gca, 'XGrid', 'on') +set(gca, 'color', [0.8, 0.8, 0.8]); plot([xlower xupper], [0 0], 'k', 'Linewidth', 1.7); +set(gca, 'FontUnits', 'Normalized'); title('oned_ensemble','Interpreter','none') +% Format the information and error messages for ensemble creation +h_click = text(0, 0.6, {'Click inside grey graphics box to create member', ... + '(only X value is used)'}, 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center', ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_err_text = text(0, -0.15, 'An ensemble has to have at least 2 members.', ... + 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center','Color', atts.red, ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); +h_finish = text(0, -0.15, 'Click outside of grey to finish', ... + 'Fontsize', atts.fontsize, 'HorizontalAlignment', 'center', ... + 'FontUnits', 'Normalized', 'Visible', 'Off'); + + %% ----------------------------------------------------------------------------- function button_create_new_ens_Callback(~,~) @@ -491,6 +523,8 @@ function button_create_new_ens_Callback(~,~) % Turn Off any old update points set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); set(handles.h_inf_up_ens, 'Visible', 'Off'); set(handles.h_inf_ens_member, 'Visible', 'Off'); @@ -509,14 +543,9 @@ function button_create_new_ens_Callback(~,~) % Messages are centered in the middle. xmid = (xupper + xlower) / 2.0; - h_click = text(xmid, 0.6, {'Click inside graphics box to create member', ... - '(only X value is used)'}, 'FontSize', atts.fontsize, 'HorizontalAlignment', 'center'); - - h_err_text = text(xmid, -0.15, 'An ensemble has to have at least 2 members.', ... - 'FontSize', atts.fontsize, 'Visible', 'on', 'HorizontalAlignment', 'center','Color', atts.red); - - h_finish = text(xmid, -0.15, 'Click outside of plot to finish', ... - 'Fontsize', atts.fontsize, 'Visible', 'Off', 'HorizontalAlignment', 'center'); + set(h_click, 'Position', [xmid, 0.6, 0], 'Visible', 'On'); + set(h_err_text, 'Position', [xmid, -0.15, 0], 'Visible', 'On'); + set(h_finish, 'Position', [xmid, -0.15, 0], 'Visible', 'Off'); ens_size = 0; @@ -625,6 +654,8 @@ function button_update_ens_Callback (~, ~) % Turn Off any old points set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); set(handles.h_inf_up_ens, 'Visible', 'Off'); set(handles.h_inf_ens_member, 'Visible', 'Off'); @@ -647,18 +678,24 @@ function button_update_ens_Callback (~, ~) case 'EAKF' [obs_increments, ~] = ... obs_increment_eakf(ensemble, handles.observation, handles.obs_error_sd^2); + + handles.h_prior_pdf = plot_gaussian(mean(ensemble), std(ensemble), 1); + set(handles.h_prior_pdf, 'linewidth', 2, 'color', atts.green); case 'EnKF' [obs_increments, ~] = ... obs_increment_enkf(ensemble, handles.observation, handles.obs_error_sd^2); + % There is no prior distribution to plot for the EnKF, it's not a QCEF case 'RHF' - [obs_increments, ~] = ... - obs_increment_rhf(ensemble, handles.observation, handles.obs_error_sd^2); + bounded_left = false; + [obs_increments, err, xp_prior, yp_prior, xp_post, yp_post] = ... + obs_increment_rhf(ensemble, handles.observation, handles.obs_error_sd^2, bounded_left); + handles.h_prior_pdf = plot(xp_prior, yp_prior, 'linewidth', 2, 'color', atts.green); + handles.h_post_pdf = plot(xp_post, yp_post, 'linewidth', 2, 'color', atts.blue); end % Add on increments to get new ensemble new_ensemble = ensemble + obs_increments; - - y(1:size(ensemble)) = -0.1; + y(1:size(ensemble,2)) = -0.1; handles.h_update_ens = plot(new_ensemble, y, '*', 'MarkerSize', 16, 'Color', atts.blue); % Plot lines connecting the prior and posterior ensemble members @@ -678,6 +715,12 @@ function button_update_ens_Callback (~, ~) str1 = sprintf('Posterior SD = %.4f',new_sd); set(handles.ui_text_post_sd, 'String', str1, 'Visible', 'on'); + + % Plot the posterior continuous for the EAKF + if(strcmp(filter_type, 'EAKF')) + handles.h_post_pdf = plot_gaussian(new_mean, new_sd, 1); + set(handles.h_post_pdf, 'linewidth', 2, 'color', atts.blue); + end % If the checkbox isn't set, return now if(not(get(handles.ui_checkbox_inflation, 'Value'))) @@ -698,9 +741,10 @@ function button_update_ens_Callback (~, ~) end % Update mean and sd of old posterior + handles.inf_prior_mean = mean(inf_ens(1:handles.ens_size)); handles.inf_prior_sd = std(inf_ens(1:handles.ens_size)); - str1 = sprintf('Inflated = %.4f',handles.prior_mean); + str1 = sprintf('Inflated = %.4f',handles.inf_prior_mean); set(handles.ui_text_inflated_prior_mean,'String',str1,'Visible','on'); str1 = sprintf('Inflated = %.4f',handles.inf_prior_sd); set(handles.ui_text_inflated_prior_sd, 'String',str1,'Visible','on'); @@ -716,13 +760,13 @@ function button_update_ens_Callback (~, ~) 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); + obs_increment_rhf(inf_ens, handles.observation, handles.obs_error_sd^2, bounded_left); end % Add on increments to get new ensemble new_ensemble = inf_ens + obs_increments; - y(1:size(ensemble)) = -0.3; + y(1:size(ensemble, 2)) = -0.3; handles.h_inf_up_ens = plot(new_ensemble, y, '*', 'MarkerSize', 16, 'Color', atts.blue); % Plot lines connecting the prior and posterior ensemble members @@ -777,6 +821,8 @@ function edit_inflation_Callback(~, ~) % Turn Off any old updated points set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); set(handles.h_inf_up_ens, 'Visible', 'Off'); set(handles.h_inf_ens_member, 'Visible', 'Off'); @@ -858,6 +904,8 @@ function edit_observation_Callback(~, ~) % Turn Off any old updated points set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); set(handles.h_inf_up_ens, 'Visible', 'Off'); set(handles.h_inf_ens_member, 'Visible', 'Off'); @@ -937,6 +985,8 @@ function edit_obs_error_sd_Callback(~, ~) % Turn Off any old updated points set(handles.h_update_ens, 'Visible', 'Off'); + set(handles.h_prior_pdf, 'Visible', 'Off'); + set(handles.h_post_pdf, 'Visible', 'Off'); set(handles.h_inf_up_ens, 'Visible', 'Off'); set(handles.h_inf_ens_member, 'Visible', 'Off'); diff --git a/guide/DART_LAB/matlab/oned_model.m b/guide/DART_LAB/matlab/oned_model.m index 203b95c3b8..2a146b385c 100644 --- a/guide/DART_LAB/matlab/oned_model.m +++ b/guide/DART_LAB/matlab/oned_model.m @@ -1,42 +1,61 @@ function oned_model -%% ONED_MODEL simple ensemble data assimilation example. +%% ONED_MODEL explores cycling ensemble data assimilation for a single variable +% model. % -% There are no input arguments. Simply invoke by typing the name. +% This application explores cycling data assimilation for a system that +% has a true value that is always zero. A single observation is generated +% for each assimilation cycle as a random draw from Normal(0, 1). % -% ONED_MODEL demonstrates the simplest possible case of ensemble data -% assimilation. It is possible to explore assimilation algorithms, -% 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 -% process is ready to repeat. +% The forecast model is a linear error growth with growth rate 2, a +% systematic bias term, and a nonlinear quadratic error term. +% The difference equation is: +% x(t+1) = x(t) + x(t) + bias + a * x(t) * |x(t)|. The values of the +% model_bias and the nonlinear coefficient a are initially set to 0 +% but can be changed using two of the boxes on the right side of the +% display. % -% ONED_MODEL opens a gui control window that plots -% the most recent prior, posterior, and observation, -% time sequences of the assimilation, the RMS error, -% spread and kurtosis, and prior and posterior rank histograms. +% The button in the upper left, initially labeled 'Advance Model' can +% be used to step through alternating model advance and assimilation +% steps; the button will be relabeled as "Assimilate Obs" to do the +% assimilation. Pressing the 'Start Auto Run' will automatically do +% alternating advance and assimilation steps. The button will be relabeled +% as 'Pause Auto Run' which can be pressed to pause the cycling. % -% The top button alternates between "Advance Model" and "Assimilate" to -% single-step the model. The "Start Auto Run" button is useful to watch -% the system evolve and generate estimates from many assimilation cycles. +% The ensemble size can be selected using the box on the right side. +% A constant multiplicative inflation value has a default of 1 (no +% inflation) but can be changed with the box on the right. % -% 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 -% Truth is not (in general) the same as the observation! +% The Ensemble Adjustment Kalman filter (EAKF), perturbed obs ensemble +% Kalman filter (EnKF), and the rank histogram filter (RHF) can be +% selected with the pushbuttons at the bottom right. % -% 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 -% in a flat rank histogram, given enough samples. +% Five panels of information are displayed about the assimilation. +% The upper left panel shows the prior ensemble members (green +% asterisks), the posterior ensemble members (blue asterisks), the +% observation (red asterisk), the likelihood distribution (red dashed +% curve), and the truth (black) for the current assimilation cycle +% plotted on a horizontal axis. % -% See also: gaussian_product.m oned_model_inf.m oned_ensemble.m -% twod_ensemble.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% The upper right panel shows the time evolution of the ensemble on +% a vertical axis. The truth (zero axis, black dashed), the prior +% ensemble (green asterisks), posterior ensemble (blue asterisks), and +% observation (red asterisk) are plotted, along with green line +% segments connecting the posterior at the previous time to the prior +% at the next time. +% +% The second panel in the left column shows the time evolution of +% ensemble mean error (blue) and ensemble spread (red) and the third +% panel shows the ensemble kurtosis. +% +% The two bottom panels show the prior (left) and posterior (right) +% rank histograms with the most recent time's entry highlighted in +% yellow. The 'Clear Hist' button on the lower right resets the rank +% histograms. The 'Reset' button restores the application to its +% initial state. +% +%% See also: bounded_oned_ensemble.m gaussian_product.m oned_cycle.m oned_ensemble.m +% oned_model_inf.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at @@ -44,11 +63,6 @@ % % DART $Id$ -%>@todo FIXME really should set plot limits based on model bias and inflation values -%> could set state limits to +/- 7 -%> observation likelihood graphic (handles.axes) can have a static range -%> based on these values as well. - help oned_model atts = stylesheet; % get the default fonts and colors @@ -316,6 +330,7 @@ 'FontSize' , atts.fontsize, ... 'FontName' , atts.fontname, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_model_bias_err_print = uicontrol('Style', 'text', ... @@ -327,6 +342,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_nonlin_err_print = uicontrol('Style', 'text', ... @@ -338,6 +354,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_ens_size_err_print = uicontrol('Style', 'text', ... @@ -349,6 +366,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); %% ----------------------------------------------------------------------------- @@ -675,7 +693,7 @@ function step_ahead(~, ~) % plot the model evolution handles.time_step = handles.time_step + 1; - plot(handles.time_step - 0.1, ens_new, '*', 'MarkerSize', 6, 'Color', atts.green); + plot(handles.time_step - 0.1, ens_new, '*', 'MarkerSize', 10, 'Color', atts.green); % Load up to plot all segments at once, more time efficient than previous loop bx(1:2, 1:handles.ens_size) = 0; @@ -705,7 +723,7 @@ function step_ahead(~, ~) handles.spread = prior_spread; L = legend([h_e h_s], 'Error', 'Spread', 'Location', 'NorthEast'); - set(L,'FontName', atts.fontname, 'FontSize', atts.fontsize,'Box','on'); + set(L, 'Box','on'); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -731,7 +749,7 @@ function step_ahead(~, ~) handles.kurtosis = prior_kurtosis; - %% Update the prior rank histogram figure + % Update the prior rank histogram figure axes(handles.h_prior_rank_histogram); ens_rank = get_ens_rank(ens_new, 0); @@ -745,8 +763,14 @@ function step_ahead(~, ~) B = bar(temp_rank,'stacked'); B(1).FaceColor= atts.blue ; B(1).EdgeColor= 'k'; B(2).FaceColor= atts.yellow ; B(2).EdgeColor= 'k'; - - %% Plot the figure window for this update + + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + title ('Prior Rank Histogram') + + + % Plot the figure window for this update axes(handles.axes); cla; @@ -793,19 +817,20 @@ function step_ahead(~, ~) base_x = ens_axis(4) - text_width; end - text(base_x, -0.1, ['Truth in Prior Bin ', num2str(ens_rank)], ... + hsss= text(base_x, -0.1, ['Truth in Prior Bin ', num2str(ens_rank)], ... 'FontSize', 14, 'FontWeight', 'Bold','FontName', atts.fontname); + set(hsss, 'FontUnits', 'Normalized', 'Fontsize', 0.05); % Draw a line from the label string to the truth h = line([base_x + text_width / 8, 0], [-0.08, -0.03]); set(h, 'Color', 'k'); % Label this plot - xlabel('State' ,'FontName', atts.fontname, 'FontSize', atts.fontsize); - title('Latest Ensemble Prior','FontName', atts.fontname, 'FontSize', atts.fontsize); + xlabel('State', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.05); + title('Latest Ensemble Prior','FontName', atts.fontname, ... + 'FontUnits', 'Normalized', 'FontSize', 0.05); L = legend([hg_prior hg_truth],'Prior','Truth'); - set(L, 'FontName', atts.fontname, 'FontSize', atts.fontsize); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -849,6 +874,20 @@ function step_ahead(~, ~) axlims(1) = handles.time_step; axlims(2) = handles.time_step + 10; axis(axlims) + + % There are normalization issues with this legend when it is resized + % from large to small. No way to independently set legend normalized + % font size. + % plot the invisible stuff and capture a nice handle array for later. + h_truth = plot(-1, 0, 'k--', 'Visible', 'on'); + h_obs = plot(-1, 0, 'r*' , 'Visible', 'on', 'MarkerSize', 10); + h_prior = plot(-1, 0, 'g*-', 'Visible', 'on', 'MarkerSize', 10, ... + 'Color', atts.green); + h_posterior = plot(-1, 0, 'b*' , 'Visible', 'on', 'MarkerSize', 10); + h_evolution_handles = [h_truth h_obs h_prior h_posterior]; + + L = legend(h_evolution_handles, 'Truth', 'Observation', 'Prior', 'Posterior'); + set(L, 'Box', 'on', 'Position',[0.821 0.770 0.118 0.148]) % Want the lower y limit to stay 0 for error spread axes(handles.h_err_spread_evolution); @@ -891,13 +930,13 @@ function step_ahead(~, ~) obs_increment_rhf(ens, observation, obs_error_sd^2); end - %% Plot the evolution of the state + % Plot the evolution of the state new_ens = ens + obs_increments; axes(handles.h_state_evolution); - plot(handles.time_step + 0.1, new_ens, 'b*', 'MarkerSize', 6); + plot(handles.time_step + 0.1, new_ens, 'b*', 'MarkerSize', 10); handles.ens = new_ens; - %% Update the rank data + % Update the rank data axes(handles.h_post_rank_histogram); ens_rank = get_ens_rank(handles.ens, 0); @@ -905,35 +944,44 @@ function step_ahead(~, ~) temp_rank(:, 1) = handles.post_rank(1:handles.ens_size + 1); temp_rank(:, 2) = 0; temp_rank(ens_rank, 2) = 1; + hold off - B = bar(temp_rank, 'stacked'); B(1).FaceColor= atts.blue ; B(1).EdgeColor= 'k'; B(2).FaceColor= atts.yellow ; B(2).EdgeColor= 'k'; + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency', 'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank', 'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + title ('Posterior Rank Histogram') % Update the permanent storage of the rank values handles.post_rank(ens_rank) = handles.post_rank(ens_rank) + 1; - %% Plot the segment for the updated error + % Plot the segment for the updated error axes(handles.h_err_spread_evolution); post_error = calculate_rmse(new_ens, 0.0); - h = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... + h_e = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... [handles.error, post_error]); - set(h,'Color', atts.blue, 'LineWidth', 2.0); + set(h_e,'Color', atts.blue, 'LineWidth', 2.0); handles.error = post_error; - %% Plot the segment for the updated spread + % Plot the segment for the updated spread axes(handles.h_err_spread_evolution); post_spread = std(new_ens); - h = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... + h_s = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... [handles.spread, post_spread]); - set(h, 'Color', atts.red, 'LineWidth', 2.0); + set(h_s, 'Color', atts.red, 'LineWidth', 2.0); handles.spread = post_spread; + + % Recreate legend when axes shift in time + if( mod(handles.time_step, 10) == 0) + L = legend([h_e h_s], 'Error', 'Spread', 'Location', 'NorthEast'); + end %% Plot the segment for the updated kurtosis axes(handles.h_kurtosis_evolution); @@ -1014,8 +1062,8 @@ function step_ahead(~, ~) base_x = ens_axis(2) - text_width; end - text(base_x, -0.18, ['Truth in Posterior Bin ', num2str(ens_rank)], ... - 'FontSize', 14, 'FontWeight', 'Bold','FontName', atts.fontname); + httt = text(base_x, -0.18, ['Truth in Posterior Bin ', num2str(ens_rank)]); + set(httt, 'FontUnits', 'Normalized', 'Fontsize', 0.05); % Draw a line from the label string to the truth plot([base_x + text_width/8, 0], [-0.16, -0.13], 'k'); @@ -1035,13 +1083,12 @@ function step_ahead(~, ~) plot(ens_axis(1:2), [-0.1 -0.1], 'k', 'LineWidth', 2); % Label this plot - xlabel('State','FontName', atts.fontname,'FontSize', atts.fontsize); - title('Latest Ensemble Prior, Likelihood, Posterior', ... - 'FontName', atts.fontname, 'FontSize', atts.fontsize,'FontWeight', 'Bold'); + xlabel('State', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.05); + title('Latest Ensemble Prior, Likelihood, Posterior','FontName', atts.fontname, ... + 'FontUnits', 'Normalized', 'FontSize', 0.05); % Put on legend L = legend([hg_prior hg_post hg_like], 'Prior', 'Posterior', 'Likelihood', 'Location','NorthEast'); - set(L, 'FontName', atts.fontname, 'FontSize', atts.fontsize); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1126,9 +1173,13 @@ function turn_on_controls () hg_like = line([0 1], [0 0.1]); set(hg_like, 'LineWidth', 2, 'Color', atts.red, 'LineStyle','--', 'Visible', 'off'); - + + % Set original horizontal axes + axis([-7 7 -Inf Inf]) + + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); + L = legend([hg_prior hg_post hg_like],'Prior','Posterior','Likelihood'); - set(L, 'FontName', atts.fontname, 'FontSize', atts.fontsize,'Box','on'); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1139,9 +1190,7 @@ function turn_on_controls () end hold on; - % Set original horizontal axes - axis([-7 7 -Inf Inf]) - + end %% ----------------------------------------------------------------------------- @@ -1164,7 +1213,7 @@ function turn_on_controls () x(1:handles.ens_size) = handles.time_step + 0.1; - plot(x, handles.ens, 'b*', 'MarkerSize', 6); + plot(x, handles.ens, 'b*', 'MarkerSize', 10); hold on str1 = '$x_{t+1} = x_t + (x_t+$model bias$) + a{\cdot}x_t{\cdot}{\mid}x_t{\mid}$'; str2 = '\hspace{1.5mm} observation is a draw from $\mathcal{N}(0,1)$'; @@ -1175,21 +1224,17 @@ function turn_on_controls () plot([1 100000], [0 0], 'k--'); % plot the invisible stuff and capture a nice handle array for later. - h_truth = plot(1, 0, 'k--', 'Visible', 'on'); - h_obs = plot(1, 0, 'r*' , 'Visible', 'on', 'MarkerSize', 10); - h_prior = plot(1, 0, 'g*-', 'Visible', 'on', 'MarkerSize', 6, 'Color', atts.green); - h_posterior = plot(1, 0, 'b*' , 'Visible', 'on', 'MarkerSize', 6); + h_truth = plot(-1, 0, 'k--', 'Visible', 'on'); + h_obs = plot(-1, 0, 'r*' , 'Visible', 'on', 'MarkerSize', 10); + h_prior = plot(-1, 0, 'g*-', 'Visible', 'on', 'MarkerSize', 10, 'Color', atts.green); + h_posterior = plot(-1, 0, 'b*' , 'Visible', 'on', 'MarkerSize', 10); h_evolution_handles = [h_truth h_obs h_prior h_posterior]; % Want the y axis limits to take care of themselves set(gca, 'YLimMode', 'Auto','XTickLabel',[],'XGrid','on'); - ylabel('State','FontName', atts.fontname,'FontSize', atts.fontsize); L = legend(h_evolution_handles, 'Truth', 'Observation', 'Prior', 'Posterior'); - set(L,'FontName', atts.fontname, ... - 'FontSize', 12, ... - 'Box', 'on', ... - 'Position',[0.821 0.770 0.118 0.148]) + set(L, 'Box', 'on', 'Position',[0.821 0.770 0.118 0.148]) if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1200,6 +1245,8 @@ function turn_on_controls () end axis([1 10 -Inf Inf]); + set(gca, 'FontUnits', 'Normalized', 'FontSize', 0.2) + ylabel('State', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.15); end @@ -1220,10 +1267,11 @@ function turn_on_controls () 'Color', 'White'); end - ylabel('Error, Spread','FontName', atts.fontname,'FontSize', atts.fontsize); axis([1 10 0 Inf]); set(gca,'XTickLabel',[],'XGrid','on') + set(gca, 'FontUnits', 'Normalized', 'FontSize', 0.1) hold on; + ylabel('Error, Spread','FontName', atts.fontname, 'FontUnits', 'Normalized','FontSize', 0.15); end @@ -1244,9 +1292,10 @@ function turn_on_controls () end axis([1 10 0 Inf]); - ylabel('Kurtosis', 'FontName', atts.fontname, 'FontSize', atts.fontsize); - xlabel('Timestep', 'FontName', atts.fontname, 'FontSize', atts.fontsize); set(gca,'XGrid', 'on') + set(gca, 'FontUnits', 'Normalized', 'FontSize', 0.1) + ylabel('Kurtosis', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.15); + xlabel('Timestep', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.15); hold on; end @@ -1264,11 +1313,12 @@ function set_prior_histogram() handles.h_prior_rank_histogram = axes('Position',[0.050 0.075 0.333 0.208]); end - ylabel('Frequency' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - xlabel('Rank' ,'FontName', atts.fontname,'FontSize', atts.fontsize); title ('Prior Rank Histogram','FontName', atts.fontname,'FontSize', atts.fontsize); axis([0 handles.ens_size+2 -Inf Inf]) set(handles.h_prior_rank_histogram,'XTick',1:(handles.ens_size+1)); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); hold on end @@ -1286,11 +1336,12 @@ function set_posterior_histogram() handles.h_post_rank_histogram = axes('Position',[0.43 0.075 0.333 0.208]); end - ylabel('Frequency' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - xlabel('Rank' ,'FontName', atts.fontname,'FontSize', atts.fontsize); title ('Posterior Rank Histogram','FontName', atts.fontname,'FontSize', atts.fontsize); axis([0 handles.ens_size+2 -Inf Inf]) set(handles.h_post_rank_histogram,'XTick',1:(handles.ens_size+1)); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); hold on; end diff --git a/guide/DART_LAB/matlab/oned_model_inf.m b/guide/DART_LAB/matlab/oned_model_inf.m index 44b216c451..860f509424 100644 --- a/guide/DART_LAB/matlab/oned_model_inf.m +++ b/guide/DART_LAB/matlab/oned_model_inf.m @@ -1,42 +1,31 @@ function oned_model_inf -%% ONED_MODEL_INF simple ensemble data assimilation example. +%% ONED_MODEL_INF extends the oned_model application to allow adaptive inflation. +% +% See the help for oned_model.m. Additional inflation options are added +% by a group of controls on the right. A toggle allows switching between +% 'Fixed Inflation' and 'Adaptive Inflation'. The adaptive inflation is +% controlled by five boxes. The minimum value of the inflation is set by +% 'Inf Min' and the maximum by 'Inf Max'. The inflation damping is set by +% 'Inf Damp'. The initial value of the inflation standard deviation is +% set by 'Inf Std Init' while the minimum standard deviation is set +% by 'Inf Std Min'. % -% There are no input arguments. Simply invoke by typing the name. +% All display panels are the same as in oned_model_mod.m except for the +% third panel from the top in the right column which displays the time +% evolution of the inflation value. The value is indicated by a blue dot, +% its standard deviation by a bracketed range, and the values of the +% mean (lambda) and standard deviation (sigma) are also printed in the +% upper left of the panel. +% The inflation algorithm is described in: +% El Gharamti, M., 2018: Enhanced Adaptive Inflation Algorithm for +% Ensemble Filters. Mon. Wea. Rev., 146, 623–640. % -% ONED_MODEL_INF demonstrates the simplest possible case of ensemble data -% assimilation. It is possible to explore assimilation algorithms, -% 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 -% process is ready to repeat. +% A log file 'oned_model_inf.log' is created when this application is run. +% It records parameters of the assimilation for each experiment run. % -% ONED_MODEL_INF opens a gui control window that plots -% the most recent prior, posterior, and observation, -% time sequences of the assimilation, the RMS error, -% spread and inflation, and prior and posterior rank histograms. -% -% The top button alternates between "Advance Model" and "Assimilate" to -% single-step the model. The "Start Auto Run" button is useful to watch -% the system evolve and generate estimates from many assimilation cycles. -% -% 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 -% Truth is not (in general) the same as the observation! -% -% 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 -% in a flat rank histogram, given enough samples. -% -% See also: gaussian_product.m oned_model_inf.m oned_ensemble.m -% twod_ensemble.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +%% See also: bounded_oned_ensemble.m gaussian_product.m oned_cycle.m oned_ensemble.m +% oned_model.m run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at @@ -454,6 +443,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_nonlin_err_print = uicontrol('Style', 'text', ... @@ -465,6 +455,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_ens_size_err_print = uicontrol('Style', 'text', ... @@ -476,6 +467,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_inf_err_print = uicontrol('Style', 'text', ... @@ -487,6 +479,7 @@ 'FontSize' , atts.fontsize, ... 'FontName' , atts.fontname, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_inf_damp_err_print = uicontrol('Style', 'text', ... @@ -498,6 +491,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_inf_min_err_print = uicontrol('Style', 'text', ... @@ -509,6 +503,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_inf_max_err_print = uicontrol('Style', 'text', ... @@ -520,6 +515,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_inf_std_err_print = uicontrol('Style', 'text', ... @@ -531,6 +527,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); handles.ui_text_inf_std_min_err_print = uicontrol('Style', 'text', ... @@ -542,6 +539,7 @@ 'FontName' , atts.fontname, ... 'FontSize' , atts.fontsize, ... 'FontWeight' , 'Bold', ... + 'FontUnits' , 'Normalized', ... 'Visible' , 'Off'); %% ----------------------------------------------------------------------------- @@ -1038,7 +1036,7 @@ function ClearStats_Callback(~, ~) axes(handles.h_err_spread_evolution); L = legend('Error','Spread','Location', 'NorthWest'); - set(L,'FontName', atts.fontname, 'FontSize', atts.fontsize, 'EdgeColor', 'w'); + set(L, 'EdgeColor', 'w'); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1219,8 +1217,7 @@ function step_ahead(~, ~) % plot the model evolution handles.time_step = handles.time_step + 1; - h_evolution.prior = plot(handles.time_step - 0.1, ens_new, '*', ... - 'MarkerSize', 6, 'Color', atts.green); + plot(handles.time_step - 0.1, ens_new, '*', 'MarkerSize', 10, 'Color', atts.green); % Load up to plot all segments at once, more time efficient than previous loop bx(1:2, 1:handles.ens_size) = 0; @@ -1251,7 +1248,25 @@ function step_ahead(~, ~) handles.error_hist(handles.time_step) = handles.error; handles.spread_hist(handles.time_step) = handles.spread; + + L = legend([h_e h_s], 'Error', 'Spread', 'Location', 'NorthEast'); + if verLessThan('matlab','R2017a') + % Convince Matlab to not autoupdate the legend with each new line. + % Before 2017a, this was the default behavior, so do nothing. + % We do not want to add the bias line to the legend, for example. + else + L.AutoUpdate = 'off'; + end + + axlims = axis; + axlims(3) = 0.0; + axis(axlims) + + + + + % Update the prior rank histogram figure axes(handles.h_prior_rank_histogram); @@ -1266,9 +1281,11 @@ function step_ahead(~, ~) B = bar(temp_rank,'stacked'); B(1).FaceColor= atts.blue ; B(1).EdgeColor= 'k'; B(2).FaceColor= atts.yellow ; B(2).EdgeColor= 'k'; - ylabel('Frequency' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - xlabel('Rank' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - title ('Prior Rank Histogram','FontName', atts.fontname,'FontSize', atts.fontsize); + + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + title('Prior Rank Histogram') % Plot the figure window for this update axes(handles.axes); @@ -1284,7 +1301,7 @@ function step_ahead(~, ~) % Put on a black axis line using data limits plot([xmin xmax], [0, 0], 'k', 'Linewidth', 2); - hold on + hold on; ens_axis = [xmin xmax -0.2 y_max + 0.02]; grid on @@ -1310,12 +1327,28 @@ function step_ahead(~, ~) base_x = ens_axis(4) - text_width; end - text(base_x, -0.1, ['Truth in Prior Bin ', num2str(ens_rank)], ... + hsss = text(base_x, -0.1, ['Truth in Prior Bin ', num2str(ens_rank)], ... 'FontSize', 14, 'FontWeight', 'Bold','FontName', atts.fontname); + set(hsss, 'FontUnits', 'Normalized', 'Fontsize', 0.05); % Draw a line from the label string to the truth h = line([base_x + text_width / 8, 0], [-0.08, -0.03]); set(h, 'Color', 'k'); + + % Label this plot + xlabel('State', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.05); + title('Latest Ensemble Prior','FontName', atts.fontname, ... + 'FontUnits', 'Normalized', 'FontSize', 0.05); + + L = legend([hg_prior hg_truth],'Prior','Truth'); + + if verLessThan('matlab','R2017a') + % Convince Matlab to not autoupdate the legend with each new line. + % Before 2017a, this was the default behavior, so do nothing. + % We do not want to add the bias line to the legend, for example. + else + L.AutoUpdate = 'off'; + end % Update the permanent storage of the rank values handles.prior_rank(ens_rank) = handles.prior_rank(ens_rank) + 1; @@ -1352,7 +1385,21 @@ function step_ahead(~, ~) axlims(1) = handles.time_step; axlims(2) = handles.time_step + 10; axis(axlims) - + + % There are normalization issues with this legend when it is resized + % from large to small. No way to independently set legend normalized + % font size. + % plot the invisible stuff and capture a nice handle array for later. + h_truth = plot(-1, 0, 'k--', 'Visible', 'on'); + h_obs = plot(-1, 0, 'r*' , 'Visible', 'on', 'MarkerSize', 10); + h_prior = plot(-1, 0, 'g*-', 'Visible', 'on', 'MarkerSize', 10, ... + 'Color', atts.green); + h_posterior = plot(-1, 0, 'b*' , 'Visible', 'on', 'MarkerSize', 10); + h_evolution_handles = [h_truth h_obs h_prior h_posterior]; + + L = legend(h_evolution_handles, 'Truth', 'Observation', 'Prior', 'Posterior'); + set(L, 'Box', 'on', 'Position',[0.821 0.770 0.118 0.148]) + % Want the lower y limit to stay 0 for error spread axes(handles.h_err_spread_evolution); cla @@ -1416,10 +1463,11 @@ function step_ahead(~, ~) [obs_increments, ~] = ... obs_increment_rhf(ens, observation, obs_error_sd^2); end - + + % Plot the evolution of the state new_ens = ens + obs_increments; axes(handles.h_state_evolution); - plot(handles.time_step + 0.1, new_ens, 'b*', 'MarkerSize', 6); + plot(handles.time_step + 0.1, new_ens, 'b*', 'MarkerSize', 10); handles.ens = new_ens; % Update the rank data @@ -1435,34 +1483,39 @@ function step_ahead(~, ~) B = bar(temp_rank, 'stacked'); B(1).FaceColor= atts.blue ; B(1).EdgeColor= 'k'; B(2).FaceColor= atts.yellow ; B(2).EdgeColor= 'k'; - ylabel('Frequency' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - xlabel('Rank' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - title ('Posterior Rank Histogram','FontName', atts.fontname,'FontSize', atts.fontsize); - + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency', 'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank', 'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + title('Posterior Rank Histogram') + % Update the permanent storage of the rank values handles.post_rank(ens_rank) = handles.post_rank(ens_rank) + 1; - %% Plot the segment for the updated error + % Plot the segment for the updated error axes(handles.h_err_spread_evolution); post_error = calculate_rmse(new_ens, 0.0); - h = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... + h_e = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... [handles.error, post_error]); - set(h,'Color', atts.blue, 'LineWidth', 2.0); + set(h_e,'Color', atts.blue, 'LineWidth', 2.0); handles.error = post_error; % Plot the segment for the updated spread - post_spread = std(new_ens); - axes(handles.h_err_spread_evolution); - h = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... + post_spread = std(new_ens); + h_s = plot([handles.time_step - 0.1, handles.time_step + 0.1], ... [handles.spread, post_spread]); - set(h, 'Color', atts.red, 'LineWidth', 2.0); + set(h_s, 'Color', atts.red, 'LineWidth', 2.0); handles.spread = post_spread; + + % Recreate legend when axes shift in time + if( mod(handles.time_step, 10) == 0) + L = legend([h_e h_s], 'Error', 'Spread', 'Location', 'NorthEast'); + end time_range = [ handles.time_last_change, handles.time_step ]; show_rms_on_plot(handles.error_hist, handles.spread_hist, time_range); @@ -1478,19 +1531,6 @@ function step_ahead(~, ~) g = errorbar(handles.time_step + 0.1, post_inflation, handles.adap_inf_Std,'-.ob','MarkerSize',8,... 'MarkerEdgeColor',atts.blue,'MarkerFaceColor',atts.blue); - L = legend( g, [ '\lambda= ' sprintf('%.4f', post_inflation) ... - ', \sigma= ' sprintf('%.4f', handles.adap_inf_Std) ], ... - 'Location', 'NorthWest'); - set(L, 'FontName', atts.fontname, 'FontSize', 14, 'EdgeColor', 'w'); - - if verLessThan('matlab','R2017a') - % Convince Matlab to not autoupdate the legend with each new line. - % Before 2017a, this was the default behavior, so do nothing. - % We do not want to add the bias line to the legend, for example. - else - L.AutoUpdate = 'off'; - end - %% Plot the figure for this update axes(handles.axes); cla; @@ -1500,7 +1540,7 @@ function step_ahead(~, ~) % Plot the observation likelihood [hg_like, ~, ylims] = plot_gaussian(observation, obs_error_sd, 1.0); set(hg_like, 'Color', atts.red, 'LineWidth', 2, 'LineStyle', '--'); - hold on + hold on; % Want axes to encompass likely values for plotted obs_likelihood % The observed value will be between -4 and 4 with very high probability, @@ -1509,6 +1549,7 @@ function step_ahead(~, ~) xmin = -10; xmax = 10; ens_axis = [xmin xmax -0.2 ylims(2)+0.02]; + axis(ens_axis); % Put on a black axis line using data limits plot([xmin xmax], [0, 0], 'k', 'Linewidth', 2); @@ -1518,14 +1559,14 @@ function step_ahead(~, ~) tick_half = 0.015; for n_tick = 1:handles.ens_size - plot([ens(n_tick), ens(n_tick)], ... + hg_prior = plot([ens(n_tick), ens(n_tick)], ... [-tick_half, tick_half], 'Color', atts.green, ... 'LineWidth', 2); end % Plot the posterior ensemble members in blue for n_tick = 1:handles.ens_size - plot([new_ens(n_tick), new_ens(n_tick)], ... + hg_post = plot([new_ens(n_tick), new_ens(n_tick)], ... [-0.1 - tick_half, -0.1 + tick_half], 'Color', atts.blue, ... 'LineWidth', 2); end @@ -1549,8 +1590,8 @@ function step_ahead(~, ~) base_x = ens_axis(2) - text_width; end - text(base_x, -0.18, ['Truth in Posterior Bin ', num2str(ens_rank)], ... - 'FontSize', 14, 'FontWeight', 'Bold','FontName', atts.fontname); + httt = text(base_x, -0.18, ['Truth in Posterior Bin ', num2str(ens_rank)]); + set(httt, 'FontUnits', 'Normalized', 'Fontsize', 0.05); % Draw a line from the label string to the truth plot([base_x + text_width/8, 0], [-0.16, -0.13], 'k'); @@ -1564,11 +1605,27 @@ function step_ahead(~, ~) set(gca, 'YTick', [0 0.1 0.2 0.3 0.4], ... 'XLim' , [def_axis(1) def_axis(2)], ... 'YLim' , [def_axis(3) def_axis(4)]); - grid on + grid on; % Plot an additional axis plot(ens_axis(1:2), [-0.1 -0.1], 'k', 'LineWidth', 2); - + + % Label this plot + xlabel('State', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.05); + title('Latest Ensemble Prior, Likelihood, Posterior','FontName', atts.fontname, ... + 'FontUnits', 'Normalized', 'FontSize', 0.05); + + % Put on legend + L = legend([hg_prior hg_post hg_like], 'Prior', 'Posterior', 'Likelihood', 'Location','NorthEast'); + + if verLessThan('matlab','R2017a') + % Convince Matlab to not autoupdate the legend with each new line. + % Before 2017a, this was the default behavior, so do nothing. + % We do not want to add the bias line to the legend, for example. + else + L.AutoUpdate = 'off'; + end + end end @@ -1662,7 +1719,7 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) str2 = sprintf('Spread: %.2f', mean(prior_aes_vals_new) ); L = legend( str1, str2, 'Location', 'NorthWest'); - set(L, 'EdgeColor', 'w', 'FontSize', atts.fontsize) + set(L, 'EdgeColor', 'w'); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1700,10 +1757,17 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) set(hg_post, 'LineWidth', 2, 'Color', atts.blue, 'Visible', 'on'); hg_like = line([0 0], [0 0]); - set(hg_like, 'LineWidth', 2, 'Color', atts.red, 'LineStyle','--', 'Visible', 'on'); + set(hg_like, 'LineWidth', 2, 'Color', atts.red, 'LineStyle','--', 'Visible', 'off'); + % Axis Limits + y_max = 1 / (sqrt(2 * pi) * handles.obs_error_sd); + xmin = -10; + xmax = 10; + axis([xmin xmax -0.2 y_max + 0.02]) + + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.05); + L = legend([hg_prior hg_post hg_like],'Prior','Posterior','Likelihood'); - set(L,'FontName', atts.fontname, 'FontSize', atts.fontsize, 'Box', 'on'); if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1712,17 +1776,8 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) else L.AutoUpdate = 'off'; end - xlabel('State', 'FontName', atts.fontname, 'FontSize', atts.fontsize); - title('Latest Ensemble Prior, Likelihood, Posterior', ... - 'FontName', atts.fontname, 'FontSize', atts.fontsize,'FontWeight', 'Bold'); - - % Axis Limits - y_max = 1 / (sqrt(2 * pi) * handles.obs_error_sd); - xmin = -10; - xmax = 10; hold on; - axis([xmin xmax -0.2 y_max + 0.02]) end @@ -1746,7 +1801,7 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) x(1:handles.ens_size) = handles.time_step + 0.1; - plot(x, handles.ens, 'b*', 'MarkerSize', 6); + plot(x, handles.ens, 'b*', 'MarkerSize', 10); hold on str1 = '$x_{t+1} = x_t + (x_t+$model bias$) + a{\cdot}x_t{\cdot}{\mid}x_t{\mid}$'; str2 = '\hspace{1.5mm} observation is a draw from $\mathcal{N}(0,1)$'; @@ -1757,21 +1812,17 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) plot([1 100000], [0 0], 'k--'); % plot the invisible stuff and capture a nice handle array for later. - h_truth = plot(1, 0, 'k--', 'Visible', 'on'); - h_obs = plot(1, 0, 'r*' , 'Visible', 'on', 'MarkerSize', 10); - h_prior = plot(1, 0, 'g*-', 'Visible', 'on', 'MarkerSize', 6, 'Color', atts.green); - h_posterior = plot(1, 0, 'b*' , 'Visible', 'on', 'MarkerSize', 6); + h_truth = plot(-1, 0, 'k--', 'Visible', 'on'); + h_obs = plot(-1, 0, 'r*' , 'Visible', 'on', 'MarkerSize', 10); + h_prior = plot(-1, 0, 'g*-', 'Visible', 'on', 'MarkerSize', 10, 'Color', atts.green); + h_posterior = plot(-1, 0, 'b*' , 'Visible', 'on', 'MarkerSize', 10); h_evolution_handles = [h_truth h_obs h_prior h_posterior]; % Want the y axis limits to take care of themselves set(gca, 'YLimMode', 'Auto','XTickLabel',[],'XGrid','on'); - ylabel('State','FontName', atts.fontname,'FontSize', atts.fontsize); L = legend(h_evolution_handles, 'Truth', 'Observation', 'Prior', 'Posterior'); - set(L, 'FontName', atts.fontname, ... - 'FontSize', 14, ... - 'Box', 'on', ... - 'Position',[0.821 0.770 0.118 0.148]) + set(L, 'Box', 'on', 'Position',[0.821 0.770 0.118 0.148]) if verLessThan('matlab','R2017a') % Convince Matlab to not autoupdate the legend with each new line. @@ -1780,7 +1831,10 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) else L.AutoUpdate = 'off'; end + axis([1 10 -10 10]); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.2) + ylabel('State','FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.15); end @@ -1806,22 +1860,12 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) h_s = line([0 0], [0 0]); set(h_s, 'LineWidth', 2, 'Color', atts.red, 'Visible', 'on'); - - ylabel('Error, Spread','FontName', atts.fontname,'FontSize', atts.fontsize); + axis([1 10 0 10]); set(gca,'XTickLabel',[],'XGrid','on') - hold on - - L = legend([h_e h_s], 'Error', 'Spread', 'Location', 'NorthWest'); - set(L, 'FontName', atts.fontname, 'FontSize', atts.fontsize, 'EdgeColor', 'w'); - - if verLessThan('matlab','R2017a') - % Convince Matlab to not autoupdate the legend with each new line. - % Before 2017a, this was the default behavior, so do nothing. - % We do not want to add the bias line to the legend, for example. - else - L.AutoUpdate = 'off'; - end + set(gca, 'FontUnits', 'Normalized', 'FontSize', 0.1) + hold on; + ylabel('Error, Spread','FontName', atts.fontname, 'FontUnits', 'Normalized','FontSize', 0.15); title('Averaging Over Steps(1:n)'); @@ -1845,11 +1889,11 @@ function show_rms_on_plot(prior_rms_vals, prior_aes_vals, ranges) plot([1 100000], [1 1], 'k:'); axis([1 10 0. 3]); - - ylabel('Inflation', 'FontName', atts.fontname, 'FontSize', atts.fontsize); - xlabel('Timestep', 'FontName', atts.fontname, 'FontSize', atts.fontsize); set(gca,'XGrid', 'on') - hold on + set(gca, 'FontUnits', 'Normalized', 'FontSize', 0.1) + ylabel('Inflation', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.15); + xlabel('Timestep', 'FontName', atts.fontname, 'FontUnits', 'Normalized', 'FontSize', 0.15); + hold on; end @@ -1865,12 +1909,13 @@ function set_prior_histogram() else handles.h_prior_rank_histogram = axes('Position',[0.050 0.075 0.333 0.208]); end - - ylabel('Frequency' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - xlabel('Rank' ,'FontName', atts.fontname,'FontSize', atts.fontsize); + title ('Prior Rank Histogram','FontName', atts.fontname,'FontSize', atts.fontsize); axis([0 handles.ens_size+2 -Inf Inf]) set(handles.h_prior_rank_histogram,'XTick',1:(handles.ens_size+1)); + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); hold on end @@ -1887,14 +1932,16 @@ function set_posterior_histogram() else handles.h_post_rank_histogram = axes('Position',[0.43 0.075 0.333 0.208]); end - - ylabel('Frequency' ,'FontName', atts.fontname,'FontSize', atts.fontsize); - xlabel('Rank' ,'FontName', atts.fontname,'FontSize', atts.fontsize); + + title ('Posterior Rank Histogram','FontName', atts.fontname,'FontSize', atts.fontsize); axis([0 handles.ens_size+2 -Inf Inf]) set(handles.h_post_rank_histogram,'XTick',1:(handles.ens_size+1)); - hold on - + set(gca, 'FontUnits', 'Normalized', 'Fontsize', 0.1); + ylabel('Frequency' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + xlabel('Rank' ,'FontName', atts.fontname,'FontUnits', 'Normalized','FontSize', 0.125); + hold on; + end %% ----------------------------------------------------------------------------- diff --git a/guide/DART_LAB/matlab/private/bnrh_cdf.m b/guide/DART_LAB/matlab/private/bnrh_cdf.m new file mode 100644 index 0000000000..5a0a7363d4 --- /dev/null +++ b/guide/DART_LAB/matlab/private/bnrh_cdf.m @@ -0,0 +1,114 @@ +function [sort_x, quantiles, tail_amp_left, tail_mean_left, tail_sd_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, ... + do_uniform_tail_left, do_uniform_tail_right] = ... + bnrh_cdf(x, ens_size, bounded_below, bounded_above, lower_bound, upper_bound) + + +%% bnrh_cdf Computes the quantiles for a bnrh distribution +% This is modified directly from the Fortran version and is not currently +% using matlab efficiently and clearly. + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% Computes all information about a rank histogram cdf given the ensemble and bounds + +% Get ensemble sd for the normal tails +tail_sd_left = std(x); +tail_sd_right = tail_sd_left; + +% Don't know what to do if sd is 0; tail_sds returned are illegal value to indicate this +if(tail_sd_left <= 0.0) + tail_sd_left = -99; + tail_sd_right = -99; + return +end + +% Sort. For now, don't worry about efficiency, but may need to somehow pass previous +% sorting indexes and use a sort that is faster for nearly sorted data. Profiling can guide the need +[sort_x, sort_index] = sort(x); + +% Fail if lower bound is larger than smallest ensemble member +if(bounded_below) + % Do in two ifs in case the bound is not defined + if(sort_x(1) < lower_bound) + %write(errstring, *) 'Smallest ensemble member less than lower bound', & + %sort_x(1), lower_bound + %call error_handler(E_ERR, 'bnrh_cdf', errstring, source) + % Just die, could be more graceful + stop + end +end + +% Fail if upper bound is smaller than the largest ensemble member +if(bounded_above) + if(sort_x(ens_size) > upper_bound) + %write(errstring, *) 'Largest ensemble member greater than upper bound', & + %sort_x(ens_size), upper_bound + %call error_handler(E_ERR, 'bnrh_cdf', errstring, source) + stop + end +end + +% The ensemble size array q contains the sorted quantiles corresponding to the sorted ensemble sort_x +q = ens_quantiles(sort_x, ens_size, bounded_below, bounded_above, lower_bound, upper_bound); +% The quantiles array has the unsorted quantiles corresponding to the unsorted input ensemble, x +for i = 1:ens_size + indx = sort_index(i); + quantiles(indx) = q(i); +end + +% Compute the characteristics of tails + +% For unit normal, find distance from mean to where cdf is 1/(ens_size+1) (del_q_. +% Saved to avoid redundant computation for repeated calls with same ensemble size +del_q = 1.0 / (ens_size + 1.8); + +% This will be negative, want it to be a distance so make it positive +dist_for_unit_sd = -1.0 * norminv(del_q, 0, 1); + +% Find a mean so that 1 / (ens_size + 1) probability is in outer regions +tail_mean_left = sort_x(1) + dist_for_unit_sd * tail_sd_left; +tail_mean_right = sort_x(ens_size) - dist_for_unit_sd * tail_sd_right; + +% If the distribution is bounded, still want 1 / (ens_size + 1) (del_q) in outer regions +% Put an amplitude term (greater than 1) in front of the tail normals +% Amplitude is 1 if there are no bounds, so start with that +tail_amp_left = 1.0; +tail_amp_right = 1.0; + +% Switch to uniform for cases where bound and outermost ensemble have close quantiles +% Default: not close +uniform_threshold = 0.01; +do_uniform_tail_left = false; +if(bounded_below) + % Compute the CDF at the bounds + bound_quantile = normcdf(lower_bound, tail_mean_left, tail_sd_left); + % Note that due to roundoff it is possible for del_q - quantile to be slightly negative + if((del_q - bound_quantile) / del_q < uniform_threshold) + % If bound and ensemble member are too close, do uniform approximation + do_uniform_tail_left = true; + else + % Compute the left tail amplitude + tail_amp_left = del_q / (del_q - bound_quantile); + end +end + +% Default: not close +do_uniform_tail_right = false; +if(bounded_above) + % Compute the CDF at the bounds + bound_quantile = normcdf(upper_bound, tail_mean_right, tail_sd_right); + % Note that due to roundoff it is possible for the numerator to be slightly negative + if((bound_quantile - (1.0 - del_q)) / del_q < uniform_threshold) + % If bound and ensemble member are too close, do uniform approximation + do_uniform_tail_right = true; + else + % Compute the right tail amplitude + tail_amp_right = del_q / (del_q - (1.0 - bound_quantile)); + end +end + diff --git a/guide/DART_LAB/matlab/private/bnrh_cdf_initialized.m b/guide/DART_LAB/matlab/private/bnrh_cdf_initialized.m new file mode 100644 index 0000000000..42c4028053 --- /dev/null +++ b/guide/DART_LAB/matlab/private/bnrh_cdf_initialized.m @@ -0,0 +1,100 @@ +function [quantile] = bnrh_cdf_initialized(x, ens_size, sort_ens, ... + bounded_below, bounded_above, lower_bound, upper_bound, ... + tail_amp_left, tail_mean_left, tail_sd_left, do_uniform_tail_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, do_uniform_tail_right, q) + + +%% bnrh_cdf_initialed Computes the quantiles for a bnrh distribution for a different +% ensemble. +% This is modified directly from the Fortran version and is not currently +% using matlab efficiently and clearly. + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% Quantile increment between ensemble members for bnrh +del_q = 1.0 / (ens_size + 1.0); + +if(x < sort_ens(1)) + % In the left tail + % Do an error check to make sure ensemble member isn't outside bounds, may be redundant + if(bounded_below & x < lower_bound) + stop + %write(errstring, *) 'Ensemble member less than lower bound', x, lower_bound + %call error_handler(E_ERR, 'bnrh_cdf_initialized', errstring, source) + % This error can occur due to roundoff in increment generation from BNRHF + % See discussion in function fix_bounds. + end + + if(do_uniform_tail_left) + % Uniform approximation for left tail; Note that denominator cannot be 0 but could be small + quantile = (x - lower_bound) / (sort_ens(1) - lower_bound) * del_q; + else + % It's a normal tail + if(bounded_below) + quantile = tail_amp_left * (normcdf(x, tail_mean_left, tail_sd_left) - ... + normcdf(lower_bound, tail_mean_left, tail_sd_left)); + else % Unbounded, tail normal goes all the way down to quantile 0, amplitude is 1 + quantile = (normal_cdf(x, tail_mean_left, tail_sd_left) / ... + normal_cdf(sort_ens(1), tail_mean_left, tail_sd_left)) * del_q; + end + % Make sure it doesn't sneak past the quantile of the smallest ensemble member due to round-off + quantile = min(quantile, q(1)); + end +elseif(x == sort_ens(1)) + % This takes care of cases where there are multiple bnrh values at the bdry or at first ensemble + quantile = q(1); +elseif(x > sort_ens(ens_size)) + % In the right tail + % Do an error check to make sure ensemble member isn't outside bounds, may be redundant + if(bounded_above & x > upper_bound) + stop + %write(errstring, *) 'Ensemble member greater than upper bound first check(see code)', x, upper_bound + %call error_handler(E_ERR, 'bnrh_cdf_initialized', errstring, source) + % This error can occur due to roundoff in increment generation from bounded BNRHF + % See discussion in function fix_bounds + end + + if(do_uniform_tail_right) + % Uniform approximation for right tail + % The division here could be a concern. However, if sort_ens(ens_size) == upper_bound, then + % x cannot be > sort_ens(ens_size). + quantile = ens_size *del_q + ... + (x - sort_ens(ens_size)) / (upper_bound - sort_ens(ens_size)) * del_q; + else + % It's a normal tail + q_at_largest_ens = normcdf(sort_ens(ens_size), tail_mean_right, tail_sd_right); + % Want to avoid quantiles exceeding 1 due to numerical issues. Do fraction of the normal part + if(bounded_above) + upper_q = tail_amp_right * normcdf(upper_bound, tail_mean_right, tail_sd_right); + fract = (tail_amp_right * normcdf(x, tail_mean_right, tail_sd_right) - ... + tail_amp_right * q_at_largest_ens) / (upper_q - tail_amp_right * q_at_largest_ens); + else + % Normal goes all the way to infinity, amplitude is 1, q at infinity is 1 + fract = (normcdf(x, tail_mean_right, tail_sd_right) - q_at_largest_ens) / ... + (1.0 - q_at_largest_ens); + end + + quantile = ens_size * del_q + fract * del_q; + quantile = min(quantile, 1.0); + end + +else + % In an interior bin + for j = 1:ens_size - 1 + if(x < sort_ens(j+1)) + % The division here could be a concern. + % However, sort_ens(j)< x < sort_ens(j+1) so the two cannot be equal + quantile = j * del_q + ... + ((x - sort_ens(j)) / (sort_ens(j+1) - sort_ens(j))) * del_q; + break + elseif(x == sort_ens(j+1)) + quantile = q(j+1); + break + end + end +end + diff --git a/guide/DART_LAB/matlab/private/ens_quantiles.m b/guide/DART_LAB/matlab/private/ens_quantiles.m new file mode 100644 index 0000000000..36450f41a7 --- /dev/null +++ b/guide/DART_LAB/matlab/private/ens_quantiles.m @@ -0,0 +1,83 @@ +function [q] = ... + ens_quantiles(sorted_ens, ens_size, bounded_below, bounded_above, lower_bound, upper_bound) + + +%% ensemble_quantiles Computes the quantiles of a prior ensemble for a bnrh distribution +% This is modified directly from the Fortran version and is not currently +% using matlab efficiently and clearly. + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% This is directly converted from fortran code +% Given sorted ensemble which may have members identical to the bounds or may contain +% duplicates, compute the quantiles for each member in an bounded normal rh distribution + +% Get number of ensemble members that are duplicates of the lower bound +lower_dups = 0; +if(bounded_below) + for i = 1:ens_size + if(sorted_ens(i) == lower_bound) + lower_dups = lower_dups + 1; + else + break; + end + end +end + +% Get number of ensemble members that are duplicates of the upper bound +upper_dups = 0; +if(bounded_above) + for i = ens_size:-1:1 + if(sorted_ens(i) == upper_bound) then + upper_dups = upper_dups + 1 + else + break; + end + end +end + +% If there are duplicate ensemble members away from the boundaries need to revise quantiles +% Make sure not to count duplicates already handled at the boundaries +% Outer loop determines if a series of duplicates starts at sorted index i +d_start = lower_dups + 1; +d_end = ens_size - upper_dups; + +% Get start, length, and end of each series of duplicates away from the bounds +series_num = 1; +series_start(series_num) = d_start; +series_length(series_num) = 1; +for i = d_start + 1:d_end + if(sorted_ens(i) == sorted_ens(i - 1)) + series_length(series_num) = series_length(series_num) + 1; + else + series_end(series_num) = i-1; + series_num = series_num + 1; + series_start(series_num) = i; + series_length(series_num) = 1; + end +end + +% Off the end, finish up the last series +series_end(series_num) = d_end; + +% Now get the value of the quantile for the exact ensemble members +% Start with the lower bound duplicates +for i = 1:lower_dups + q(i) = lower_dups / (2.0 * (ens_size + 1.0)); +end + +% Top bound duplicates next +for i = ens_size - upper_dups + 1:ens_size + q(i) = 1.0 - upper_dups / (2.0 * (ens_size + 1.0)); +end + +% Do the interior series +for i = 1:series_num + for j = series_start(i):series_end(i) + q(j) = series_start(i) / (ens_size + 1.0) + (series_length(i) - 1.0) / (2.0 * (ens_size + 1.0)); + end +end diff --git a/guide/DART_LAB/matlab/private/gamma_ppi_update.m b/guide/DART_LAB/matlab/private/gamma_ppi_update.m new file mode 100644 index 0000000000..41ae71c039 --- /dev/null +++ b/guide/DART_LAB/matlab/private/gamma_ppi_update.m @@ -0,0 +1,50 @@ +function [post_state, prior_obs_ppi, post_obs_ppi, prior_state_ppi, post_state_ppi] = ... + gamma_ppi_update(prior_obs, prior_state, post_obs) + +%% gamma_ppi_update Computes increments for unobserved variable with gamma transform +% Need to discuss the available options eventually +% For now this implements the default options + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% Do the transform for the normally distributed observed ensemble (just a normalizaion) +prior_obs_mean = mean(prior_obs); +prior_obs_sd = std(prior_obs); +% Note that both are transformed with the prior statistics as per QCEFF algorithm +prior_obs_ppi = (prior_obs - prior_obs_mean) / prior_obs_sd + prior_obs_mean; +post_obs_ppi = (post_obs - prior_obs_mean) / prior_obs_sd + prior_obs_mean; + +% Get the obs_increments in the transformed space +obs_increments = post_obs_ppi - prior_obs_ppi; + +% Get the shape and scale of the prior state ensemble +pgamma_params = gamfit(prior_state); +pgamma_shape = pgamma_params(1); pgamma_scale = pgamma_params(2); + +% Compute the quantiles of the prior state members +prior_state_q = gamcdf(prior_state, pgamma_shape, pgamma_scale); + +% Do a probit transform on the quantiles +prior_state_ppi = norminv(prior_state_q, 0, 1); + +% Update the unobserved variable +covar = cov(prior_obs_ppi, prior_state_ppi); +state_inc = obs_increments * covar(1, 2) / covar(1, 1); +post_state_ppi = prior_state_ppi + state_inc; + +% Transform back to gamma cdf space +post_state_q = normcdf(post_state_ppi, 0, 1); + +% Transform back to regular space +post_state = gaminv(post_state_q, pgamma_shape, pgamma_scale); + +%----------------------------------------------- + +% +% $URL$ +% $Revision$ +% $Date$ diff --git a/guide/DART_LAB/matlab/private/inflate_bnrh.m b/guide/DART_LAB/matlab/private/inflate_bnrh.m new file mode 100644 index 0000000000..deee174619 --- /dev/null +++ b/guide/DART_LAB/matlab/private/inflate_bnrh.m @@ -0,0 +1,39 @@ +function [inf_ens] = inflate_bnrh(ens, ens_size, var_inf, ... + bounded_below, bounded_above, lower_bound, upper_bound) + +%% inflate_bnrh Inflates an ensemble in a bnrh/probit transformed space + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% The sd inflation +inf = sqrt(var_inf); + +% Compute the quantiles of the prior ensemble members +[sort_x, prior_q, tail_amp_left, tail_mean_left, tail_sd_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, ... + do_uniform_tail_left, do_uniform_tail_right] = ... + bnrh_cdf(ens, ens_size, bounded_below, bounded_above, lower_bound, upper_bound); + +% Do a probit transform on the quantiles +probit_ens = norminv(prior_q, 0, 1); + +% Inflate in probit space +probit_inf_ens = zeros(1, ens_size); +probit_mean = mean(probit_ens); +for i = 1:ens_size + probit_inf_ens(i) = (probit_ens(i) - probit_mean) * inf + probit_mean; +end + +% Transform back to gamma cdf space +inf_prior_q = normcdf(probit_inf_ens, 0, 1); + +% Transform back to regular space +inf_ens = inv_bnrh_cdf(inf_prior_q, ens_size, sort_x, ... + bounded_below, bounded_above, lower_bound, upper_bound, ... + tail_amp_left, tail_mean_left, tail_sd_left, do_uniform_tail_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, do_uniform_tail_right); + diff --git a/guide/DART_LAB/matlab/private/inflate_gamma.m b/guide/DART_LAB/matlab/private/inflate_gamma.m new file mode 100644 index 0000000000..84a9a59058 --- /dev/null +++ b/guide/DART_LAB/matlab/private/inflate_gamma.m @@ -0,0 +1,37 @@ +function [inf_ens] = inflate_gamma(ens, ens_size, var_inf) + +%% inflate_gamma Inflates an ensemble in a gamma/probit transformed space + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% The sd inflation +inf = sqrt(var_inf); + +% Get the shape and scale of the prior +pgamma_params = gamfit(ens); +pgamma_shape = pgamma_params(1); pgamma_scale = pgamma_params(2); + +% Compute the quantiles of the prior ensemble members +prior_q = gamcdf(ens, pgamma_shape, pgamma_scale); + +% Do a probit transform on the quantiles +probit_ens = norminv(prior_q, 0, 1); + +% Inflate in probit space +probit_inf_ens = zeros(1, ens_size); +probit_mean = mean(probit_ens); +for i = 1:ens_size + probit_inf_ens(i) = (probit_ens(i) - probit_mean) * inf + probit_mean; +end + +% Transform back to gamma cdf space +inf_prior_q = normcdf(probit_inf_ens, 0, 1); + +% Transform back to regular space +inf_ens = gaminv(inf_prior_q, pgamma_shape, pgamma_scale); + + diff --git a/guide/DART_LAB/matlab/private/inv_bnrh_cdf.m b/guide/DART_LAB/matlab/private/inv_bnrh_cdf.m new file mode 100644 index 0000000000..c9372f3795 --- /dev/null +++ b/guide/DART_LAB/matlab/private/inv_bnrh_cdf.m @@ -0,0 +1,96 @@ +function [x] = inv_bnrh_cdf(quantiles, ens_size, sort_ens, ... + bounded_below, bounded_above, lower_bound, upper_bound, ... + tail_amp_left, tail_mean_left, tail_sd_left, do_uniform_tail_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, do_uniform_tail_right) + +%% inv_bnrh_cdf Computes the inverse cdf for a bnrh distribution +% This is modified directly from the Fortran version and is not currently +% using matlab efficiently and clearly. + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +% Quantile increment between ensemble members for bnrh +del_q = 1.0 / (ens_size + 1.0); + +for i = 1:ens_size + q(i) = i * del_q; +end + +% Loop through each ensemble member to find posterior state +for i = 1:ens_size + curr_q = quantiles(i); + % Which region is this quantile in? + % BNRH quantiles are uniform; finding region for this quantile is trivial + region = floor(curr_q * (ens_size + 1.0)); + % Careful about numerical issues moving outside of region [0 ens_size] + if(region < 0) region = 0; end + if(region > ens_size) region = ens_size; end + + if(region == 0) + % Lower tail + if(bounded_below & do_uniform_tail_left) + % Lower tail uniform + upper_state = sort_ens(1); + x(i) = lower_bound + (curr_q / q(1)) * (upper_state - lower_bound); + else + % Find the mass at the lower bound (which could be unbounded) + if(bounded_below) + lower_mass = tail_amp_left * ... + normcdf(lower_bound, tail_mean_left, tail_sd_left); + else + lower_mass = 0.0; + end + % Find the mass at the upper bound (ensemble member 1) + upper_mass = tail_amp_left * ... + normcdf(sort_ens(1), tail_mean_left, tail_sd_left); + % What fraction of this mass difference should we go? + fract = curr_q / q(1); + target_mass = lower_mass + fract * (upper_mass - lower_mass); + x(i) = weighted_norm_inv(tail_amp_left, tail_mean_left, ... + tail_sd_left, target_mass); + end + + elseif(region == ens_size) + % Upper tail + if(bounded_above & do_uniform_tail_right) + % Upper tail is uniform + lower_state = sort_ens(ens_size); + upper_state = upper_bound; + x(i) = lower_state + (curr_q - q(ens_size)) * ... + (upper_state - lower_state) / (1.0 - q(ens_size)); + else + % Upper tail is (bounded) normal + % Find the mass at the upper bound (which could be unbounded) + if(bounded_above) then + upper_mass = tail_amp_right * ... + normcdf(upper_bound, tail_mean_right, tail_sd_right); + else + upper_mass = 1.0; + end + % Find the mass at the lower edge of the region (ensemble member n) + lower_mass = tail_amp_right * ... + normcdf(sort_ens(ens_size), tail_mean_right, tail_sd_right); + % What fraction of the last interval do we need to move + fract = (curr_q - q(ens_size)) / (1.0 - q(ens_size)); + target_mass = lower_mass + fract * (upper_mass - lower_mass); + x(i) = weighted_norm_inv(tail_amp_right, tail_mean_right, ... + tail_sd_right, target_mass); + end + + else + % Interior region; get the quantiles of the region boundary + lower_q = q(region); + upper_q = q(region + 1); + x(i) = sort_ens(region) + ((curr_q - lower_q) / (upper_q - lower_q)) * ... + (sort_ens(region + 1) - sort_ens(region)); + end + + % Imprecision can lead to x being slightly out of bounds, fix it to bounds + % Not needed for matlab test scripting + %call check_bounds(x(i), curr_q, bounded_below, lower_bound, & + %bounded_above, upper_bound, 'inf_bnrh_cdf') +end diff --git a/guide/DART_LAB/matlab/private/obs_increment_rhf.m b/guide/DART_LAB/matlab/private/obs_increment_rhf.m index 9ded3c55d6..f686775b29 100644 --- a/guide/DART_LAB/matlab/private/obs_increment_rhf.m +++ b/guide/DART_LAB/matlab/private/obs_increment_rhf.m @@ -1,4 +1,6 @@ -function [obs_increments, err] = obs_increment_rhf(ensemble, observation, obs_error_var) +function [obs_increments, err, xp_prior, yp_prior, xp_post, yp_post] = ... + obs_increment_rhf(ensemble, observation, obs_error_var, bounded_left) + %% obs_increment_rhf Computes increments for a rank histogram filter % Need to discuss the available options eventually % For now this implements the default options @@ -9,6 +11,13 @@ % % DART $Id$ +arguments + ensemble (1, :) double + observation (1, 1) double + obs_error_var (1, 1) double + bounded_left (1, 1) logical = false +end + % Set error return to default successful err = 0; @@ -42,10 +51,28 @@ left_var = prior_var; left_sd = prior_sd; +% If the distribution is bounded on the left, the amplitude is greater than one +if(bounded_left) + cdf_at_ens = 1 / (ens_size + 1); + cdf_at_bound = normcdf(0, left_mean, left_sd); + left_prior_amp = cdf_at_ens / (cdf_at_ens - cdf_at_bound); +else + left_prior_amp = 1; +end + + % 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_prior_amp = 1; + +% Define the mass in each of the bins, uniform for prior +mass(1:ens_size + 1) = 1 / (ens_size + 1); + +% Get the points to plot the prior +[xp_prior, yp_prior] = plot_rhf_pdf(x, ens_size, mass, left_mean, left_sd, left_prior_amp, ... + right_mean, right_sd, right_prior_amp, bounded_left); % Eventually want to support options, for now gaussian_likelihood_tails = false; @@ -95,7 +122,7 @@ % Get the weight for the final normalized tail gaussians % This is the same as left_amp=(ens_size + 1)*nmass(1) -left_amp = prod_weight_left / mass_sum; +left_amp = left_prior_amp * prod_weight_left / mass_sum; % This is the same as right_amp=(ens_size + 1)*nmass(ens_size+1) right_amp = prod_weight_right / mass_sum; @@ -118,9 +145,14 @@ % If it is in the inner or outer range have to use normal if(umass < cumul_mass(2)) % It is in the left tail - % Get position of x in weighted gaussian where the cdf has value umass - new_ens(i) = weighted_norm_inv(left_amp, new_mean_left, ... - new_sd_left, umass); + if(bounded_left) + % Bounded, how much mass is in the tail normal to the left of the bound + lower_mass = left_amp * normcdf(0, new_mean_left, new_sd_left); + new_ens(i) = weighted_norm_inv(left_amp, new_mean_left, new_sd_left, umass + lower_mass); + else + % Get position of x in weighted gaussian where the cdf has value umass + new_ens(i) = weighted_norm_inv(left_amp, new_mean_left, new_sd_left, umass); + end elseif (umass > cumul_mass(ens_size + 1)) % It's in the right tail % Get the position of x in weighted gaussian where the cdf has value umass @@ -159,25 +191,12 @@ obs_increments(e_ind(i)) = new_ens(i) - x(i); end +% Get the points to plot the posterior PDF +[xp_post, yp_post] = plot_rhf_pdf(x, ens_size, nmass, left_mean, left_sd, left_amp, ... + right_mean, right_sd, right_amp, bounded_left); %----------------------------------------------- - -function [x] = weighted_norm_inv(alpha, mean, sd, p) - -% Find the value of x for which the cdf of a N(mean, sd) -% multiplied times alpha has value p. - -% Can search in a standard normal, then multiply by sd at end and add mean -% Divide p by alpha to get the right place for weighted normal -np = p / alpha; - -% Find spot in standard normal -x = norm_inv(np); - -% Add in the mean and normalize by sd -x = mean + x * sd; - % % $URL$ % $Revision$ diff --git a/guide/DART_LAB/matlab/private/plot_gamma.m b/guide/DART_LAB/matlab/private/plot_gamma.m new file mode 100644 index 0000000000..b289766d57 --- /dev/null +++ b/guide/DART_LAB/matlab/private/plot_gamma.m @@ -0,0 +1,30 @@ +function [plot_handle] = plot_gamma(shape, scale) +%% plot_gamma Plot gamma over 5 standard deviations + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + + +% Get the mean and variance for this shape and scale +[gmean, gvar] = gamstat(shape, scale); +gsd = sqrt(gvar); + +% Get points to a range where pdf is small +x_min = 0; +x_max = gmean + 5*gsd; + +% Number of points is 1001 +num_points = 1001; +interval = x_max / num_points; +x = 0:interval:x_max; +y = gampdf(x, shape, scale); + +plot_handle = plot(x, y); + +% +% $URL$ +% $Revision$ +% $Date$ diff --git a/guide/DART_LAB/matlab/private/plot_rhf_pdf.m b/guide/DART_LAB/matlab/private/plot_rhf_pdf.m new file mode 100644 index 0000000000..21e5039f81 --- /dev/null +++ b/guide/DART_LAB/matlab/private/plot_rhf_pdf.m @@ -0,0 +1,56 @@ +function [xp, yp] = plot_rhf_pdf(x, ens_size, mass, left_mean, left_sd, left_amp, ... + right_mean, right_sd, right_amp, bounded_left) + +% Creates the x and y point vectors to plot a rhf PDF +% Bounds of plotting domain are given by xlow and xhigh +% x is the ensemble values +% mass is the amount of probability mass in each of the ens_size+1 bins +% left_mean, left_sd, and left_amp are the left tail PDF, similar for right + +% Use the tail distributions to get 5 standard deviations on the plots +if(bounded_left) + xlow = 0.0; +else + xlow = left_mean - 5.0 * left_sd; +end + +xhigh = right_mean + 5.0 * right_sd; + +% Point spacing +point_del = (xhigh - xlow) / 1000; +% Plot the left tail +xptsb = xlow:point_del:x(1); +yptsb = left_amp * normpdf(xptsb, left_mean, left_sd); + +% Do the interior bins +for i = 1:ens_size - 1 + % Do the vertical line on the right + height(i+1) = mass(i+1) / (x(i + 1) - x(i)); + + xpts(3*i - 2) = x(i); + xpts(3*i - 1) = x(i); + xpts(3*i) = x(i + 1); + + if(i == 1) + ypts(3*i -2) = left_amp * normpdf(x(1), left_mean, left_sd); + else + ypts(3*i -2) = height(i); + end + ypts(3*i - 1) = height(i + 1); + ypts(3*i) = height(i+1); +end + +% Plot the vertical on right of last interior bin +xptsv(1:2) = x(ens_size); +yptsv(1) = height(ens_size); +yptsv(2) = right_amp * normpdf(x(ens_size), right_mean, right_sd); + +% Plot the right tail +xptsa = x(ens_size):point_del:xhigh; +yptsa = right_amp * normpdf(xptsa, right_mean, right_sd); + +% Concatenate the different pieces +xp = [xptsb xpts xptsv xptsa]; +yp = [yptsb ypts yptsv yptsa]; + +end diff --git a/guide/DART_LAB/matlab/private/ppi_update.m b/guide/DART_LAB/matlab/private/ppi_update.m new file mode 100644 index 0000000000..b1d2a87b98 --- /dev/null +++ b/guide/DART_LAB/matlab/private/ppi_update.m @@ -0,0 +1,124 @@ +function [post_state, prior_obs_ppi, post_obs_ppi, prior_state_ppi, post_state_ppi] = ... + ppi_update(prior_obs, prior_state, post_obs, state_dist_type, obs_dist_type) + +%% ppi_update Computes increments for unobserved variable with gamma transform +% Need to discuss the available options eventually +% For now this implements the default options + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +ens_size = size(prior_state, 2); + +% Convert the prior_obs and posterior_obs to PPI space +% First get quantiles +switch(obs_dist_type) + case 'Normal' + % This is just a scaling and removing the mean but done explicitly here + prior_obs_mean = mean(prior_obs); + prior_obs_sd = std(prior_obs); +% Note that both are transformed with the prior statistics as per QCEFF algorithm + prior_obs_q = normcdf(prior_obs, prior_obs_mean, prior_obs_sd); + post_obs_q = normcdf(post_obs, prior_obs_mean, prior_obs_sd); + + case 'RHF' + % Get the quantiles for this ensemble for an unbounded BNRH distribution + % Everything except the + [sort_x, prior_obs_q, tail_amp_left, tail_mean_left, tail_sd_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, ... + do_uniform_tail_left, do_uniform_tail_right] = bnrh_cdf(prior_obs, ... + ens_size, false, false, -99, -99); + + % Need to use the prior_obs distribution for the quantile transform + for i = 1:ens_size + post_obs_q(i) = bnrh_cdf_initialized(post_obs(i), ens_size, sort_x, ... + false, false, -99, -99, ... + tail_amp_left, tail_mean_left, tail_sd_left, do_uniform_tail_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, do_uniform_tail_right, prior_obs_q); + end + +end + +% Do the probit transform +prior_obs_ppi = norminv(prior_obs_q, 0, 1); +post_obs_ppi = norminv(post_obs_q, 0, 1); + +% Get the obs_increments in the transformed space +obs_increments = post_obs_ppi - prior_obs_ppi; + + +% Get the quantiles for the appropriate distribution +switch(state_dist_type) + case 'Normal' + % Get the mean and std + pmean = mean(prior_state); + psd = std(prior_state); + % Get quantiles + prior_state_q = normcdf(prior_state, pmean, psd); + + case 'Gamma' + % Get the shape and scale of the prior state ensemble + pgamma_params = gamfit(prior_state); + pgamma_shape = pgamma_params(1); pgamma_scale = pgamma_params(2); + + % Compute the quantiles of the prior state members + prior_state_q = gamcdf(prior_state, pgamma_shape, pgamma_scale); + + case 'RHF' + % Get the quantiles for this ensemble for an unbounded BNRH distribution + [sort_x, prior_state_q, tail_amp_left, tail_mean_left, tail_sd_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, ... + do_uniform_tail_left, do_uniform_tail_right] = bnrh_cdf(prior_state, ... + ens_size, false, false, -99, -99); + + case 'BNRH' + % Get the quantiles for this ensemble for an unbounded BNRH distribution + [sort_x, prior_state_q, tail_amp_left, tail_mean_left, tail_sd_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, ... + do_uniform_tail_left, do_uniform_tail_right] = bnrh_cdf(prior_state, ... + ens_size, true, false, 0, -99); + +end + +% Do a probit transform on the quantiles +prior_state_ppi = norminv(prior_state_q, 0, 1); + +% Update the unobserved variable +covar = cov(prior_obs_ppi, prior_state_ppi); +state_inc = obs_increments * covar(1, 2) / covar(1, 1); +post_state_ppi = prior_state_ppi + state_inc; + +% Transform back to gamma cdf space +post_state_q = normcdf(post_state_ppi, 0, 1); + +% Transform back to regular space +switch(state_dist_type) + case 'Normal' + post_state = norminv(post_state_q, pmean, psd); + + case 'Gamma' + post_state = gaminv(post_state_q, pgamma_shape, pgamma_scale); + + case 'RHF' + post_state = inv_bnrh_cdf(post_state_q, ens_size, sort_x, ... + false, false, -99, -99, ... + tail_amp_left, tail_mean_left, tail_sd_left, do_uniform_tail_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, do_uniform_tail_right); + + case 'BNRH' + post_state = inv_bnrh_cdf(post_state_q, ens_size, sort_x, ... + true, false, 0, -99, ... + tail_amp_left, tail_mean_left, tail_sd_left, do_uniform_tail_left, ... + tail_amp_right, tail_mean_right, tail_sd_right, do_uniform_tail_right); +end + + +%----------------------------------------------- + +% +% $URL$ +% $Revision$ +% $Date$ diff --git a/guide/DART_LAB/matlab/private/weighted_norm_inv.m b/guide/DART_LAB/matlab/private/weighted_norm_inv.m new file mode 100644 index 0000000000..35664294de --- /dev/null +++ b/guide/DART_LAB/matlab/private/weighted_norm_inv.m @@ -0,0 +1,17 @@ + +function [x] = weighted_norm_inv(alpha, mean, sd, p) + +% Find the value of x for which the cdf of a N(mean, sd) +% multiplied times alpha has value p. + +% Can search in a standard normal, then multiply by sd at end and add mean +% Divide p by alpha to get the right place for weighted normal +np = p / alpha; + +% Find spot in standard normal +x = norm_inv(np); + +% Add in the mean and normalize by sd +x = mean + x * sd; + + diff --git a/guide/DART_LAB/matlab/run_lorenz_63.m b/guide/DART_LAB/matlab/run_lorenz_63.m index ff8ba2b42b..7be8323be6 100644 --- a/guide/DART_LAB/matlab/run_lorenz_63.m +++ b/guide/DART_LAB/matlab/run_lorenz_63.m @@ -7,9 +7,18 @@ % % This is another 'perfect_model' experiment. As the model is % advanced, observations are taken from the true state and are -% assimilated by the ensemble members. The true trajectory, the -% observations, the observation increments, and the Prior and Posterior -% ensemble states are displayed. +% assimilated by the ensemble members. The three model components +% are observed at each assimilation cycle with a specified +% observation error distribution that is Normal(0, 1). +% +% See the documentation for oned_model_mod.m for a description of +% using the 'Advance Model' button and 'Start Auto Run' button +% to do cycling assimilation. +% +% No assimilation, the Ensemble Adjustment Kalman filter (EAKF), +% perturbed obs ensemble Kalman filter (EnKF), and the rank +% histogram filter (RHF) can be selected with the pushbuttons +% at the right. The 'Reset' button restarts the application. % % To provide context and highlight the details of the assimilation, % two views are presented. The larger view is a detailed view of @@ -19,6 +28,17 @@ % on the Z-axis although they can be rotated when the assimilation is % stopped. % +% The detailed view includes the truth (black asterisk) and its +% recent trajectory (black line). It also contains the recent +% trajectory of the 20 ensemble members (green lines). If +% assimilation is being done, it also displays the observation +% (red asterisk) and the increments to the prior ensemble members +% resulting from the assimilation (red lines). +% +% The global panel plots a time history of the true trajectory +% (black line with asterisks at start and end) and the most +% recent observation (red asterisk) when assimilation is being done. +% % After you get the feel for a few single steps through the process % (by repeatedly pressing the 'Advance/Assimilate' button), select % 'No Assimilation' and start a free run. After a while some of the @@ -30,8 +50,9 @@ % particularly difficult - especially in the highly nonlinear % 'saddle' region. % -% See also: gaussian_product.m oned_model.m oned_model_inf.m oned_ensemble.m -% twod_ensemble.m run_lorenz_96.m run_lorenz_96_inf.m +%% See also: bounded_oned_ensemble.m gaussian_product.m oned_cycle.m oned_ensemble.m +% oned_model.m oned_model_inf.m run_lorenz_96.m run_lorenz_96_inf.m +% twod_ensemble.m twod_ppi_ensemble.m %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at diff --git a/guide/DART_LAB/matlab/twod_ensemble.m b/guide/DART_LAB/matlab/twod_ensemble.m index 3b732aba72..3aeb6f074e 100644 --- a/guide/DART_LAB/matlab/twod_ensemble.m +++ b/guide/DART_LAB/matlab/twod_ensemble.m @@ -1,26 +1,46 @@ function twod_ensemble -%% TWOD_ENSEMBLE demonstrates a powerful aspect of ensemble data assimilation. -% The ensemble makes it possible to estimate the impact of observations -% on unobserved state variables. +%% TWOD_ENSEMBLE explores how and unobserved variable is updated by an +% observed variable in an ensemble filter. % % Click on the 'Create New Ensemble' button to activate the interactive -% observation generation mechanism and lay down a set of ensemble -% samples of an unobserved variable (vertical axis) and an observed -% variable (horizontal axis). The ensemble members are created by -% left clicking in the central portion of the figure window. -% Start out small, say 6 or so. -% In this case, some H() operator would generate the Observed Quantity. -% The Unobserved State Variable could simply be some portion of the -% model state. +% ensemble generation mechanism. Ensemble members are created by +% left clicking in the central portion of the plot on the left side +% of the panel. The horizontal axis is for an observed quantity, +% and the vertical axis for an unobserved quantity. When you have all +% the ensemble members you want, click outside of the grey plot area. +% As you add ensemble members, marginal plots for the observed and +% unobserved axes are created and the sample correlation is displayed. % -% After creating the ensemble, the correlation between the Observed -% Quantity and the Unobserved State Variable is calculated. -% Select an assimilation algorithm and click 'Update Ensemble'. -% The increments are shown as red lines, the new Posterior estimates -% are in blue. +% Once you have created a prior ensemble (green asterisks), click the +% 'Update Ensemble' button. With the default settings, this will apply the EAKF +% algorithm to produce a posterior ensemble (blue asterisks) in the central +% bivariate plot and in the marginal plots. The prior and posteriors are +% connected by cyan segments. % -% See also: gaussian_product.m oned_model.m oned_model_inf.m oned_ensemble.m -% run_lorenz_63.m run_lorenz_96.m run_lorenz_96_inf.m +% Two other ensemble filter variants, the EnKF (sometimes referred to as +% the perturbed observations ensemble Kalman filter) and the rank histogram +% filter (RHF) can be selected with the pushbuttons at the lower right. +% Selecting one of these and pressing 'Update Ensemble' will produce +% the posterior ensemble and posterior statistics using the selected +% filter algorithm. The EnKF is a stochastic algorithm so repeated +% updates can be done for the same prior and observation by repeatedly +% pressing 'Update Ensemble'. +% +% The marginal distributions of the observed quantity and the +% unobserved quantity are shown in the rectangular plots abutting the +% square joint distribution. The prior and posterior ensembles and +% the increments between them are highlighted in all three plots after +% the 'Update Ensemble' button is pushed. +% +% The plot at the lower right shows more detail of the update for the +% observed variable including the likelihood in red. +% +% The mean and standard deviation of the likelihood can be +% changed in the red box. +% +% See also: bounded_oned_ensemble.m gaussian_product.m oned_cycle.m oned_ensemble.m +% oned_model.m oned_model_inf.m run_lorenz_63.m run_lorenz_96.m +% run_lorenz_96_inf.m twod_ppi_ensemble.m %% DART software - Copyright UCAR. This open source software is provided % by UCAR, "as is", without charge, subject to all terms of use at @@ -230,7 +250,8 @@ 'FontName', atts.fontname, ... 'FontSize', atts.fontsize); axis([0 1 0 10]); -ylabel('Unobserved State Variable', 'Fontsize', atts.fontsize); +set(gca, 'FontUnits', 'Normalized'); +ylabel('Unobserved State Variable', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); hold on %% Get a subplot for the joint @@ -241,28 +262,30 @@ axis([0 10 0 10]); grid on hold on -title('Joint Distribution'); - -h_click = text(5, 9, {'Click inside Joint Distribution plot','to create member.'}, ... - 'FontSize', 16, 'HorizontalAlignment','center'); -h_finish = text(5, 8, 'Click outside of plot to finish.', ... - 'FontSize', 16, 'HorizontalAlignment','center'); +title('Joint Distribution', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); +set(gca, 'color', [0.8, 0.8, 0.8]) +set(gca, 'FontUnits', 'Normalized') + +h_click = text(5, 9, {'Click inside Joint Distribution plot','grey area to create member.'}, ... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); +h_finish = text(5, 8, 'Click outside of grey to finish.', ... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); h_err_text = text(5, 4, 'An ensemble must have at least 2 members.', 'Color', atts.red, ... - 'FontSize', 16, 'HorizontalAlignment','center'); + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); handles.h_correl = text(5, 9, ' ', 'Color', atts.green, 'FontWeight', 'Bold', ... - 'FontSize', 16, 'HorizontalAlignment','center'); + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); %% Create a subplot for the observed variable marginal handles.h_obMarginal = axes('Position', plotposition(3,:), ... - 'FontName', atts.fontname, ... - 'FontSize', atts.fontsize); + 'FontName', atts.fontname, 'FontSize', atts.fontsize); % May want to mess with vertical axis for prior density axis([0 10 0 1]); hold on -xlabel('Observed Quantity'); +xlabel('Observed Quantity', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); +set(gca, 'FontUnits', 'Normalized'); handles.h_obs_marg = plot(observation, 0, '*', 'MarkerSize', 16, 'LineWidth', 2.0); set(handles.h_obs_marg,'Color',atts.red) @@ -286,9 +309,10 @@ set(handles.h_marg_obs_plot, 'Color', atts.red, ... 'Linestyle', '--', ... 'LineWidth', 2); -xlabel('Observed Quantity', 'FontSize', atts.fontsize); -ylabel('Observation Likelihood', 'FontSize', atts.fontsize); -title('Marginal Distribution of Observation', 'FontSize', atts.fontsize); +set(gca, 'FontUnits', 'Normalized'); +xlabel('Observed Quantity', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); +ylabel('Observation Likelihood', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); +title('Marginal Distribution of Observation', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); hold on % Plot an asterisk for the observed value @@ -330,6 +354,11 @@ function create_ensemble_Callback(~,~) % Clear out the old best fit line set(handles.h_best_fit, 'Visible', 'off'); set(handles.h_correl, 'Visible', 'off'); + + % Turn on the data entry messages + set(h_click, 'Visible', 'on'); + set(h_finish, 'Visible', 'on'); + set(h_err_text, 'Visible', 'on'); % Work in the joint distribution plot axes(handles.h_joint); @@ -408,10 +437,11 @@ function create_ensemble_Callback(~,~) prior_mean = mean(x, 2); prior_cov = cov(x(1, :), x(2, :)); slope = prior_cov(1, 2) / var(x(1, :)); + intercept = prior_mean(2) - slope * prior_mean(1); - best_x = [0 10]; - best_y(1) = prior_mean(2) - (prior_mean(1)) * slope; - best_y(2) = best_y(1) + 10 * slope; + best_x = [-1000 1000]; + best_y = slope * best_x + intercept; + handles.h_best_fit = plot(best_x, best_y, 'g', 'LineWidth', 2.0); set(handles.h_best_fit, 'Color', atts.green); diff --git a/guide/DART_LAB/matlab/twod_ppi_ensemble.m b/guide/DART_LAB/matlab/twod_ppi_ensemble.m new file mode 100644 index 0000000000..db89d37a7c --- /dev/null +++ b/guide/DART_LAB/matlab/twod_ppi_ensemble.m @@ -0,0 +1,950 @@ +function twod_ppi_ensemble +%% TWOD_PPI_ENSEMBLE demonstrates capabilities of a variety of QCEFF filters +% for computing the impact of an observed variable on an unobserved +% variable. +% +% Click on the 'Create New Ensemble' button to activate the interactive +% ensemble generation mechanism. Ensemble members are created by +% left clicking in the central portion of the plot on the left side +% of the panel. The horizontal axis is for an observed quantity, +% and the vertical axis for an unobserved quantity. The unobserved +% quantity is nonnegative so ensemble members must not be below +% the axis. When you have all the ensemble members you want, click +% outside of the grey plot area. As you add ensemble members, marginal +% plots for the observed and unobserved axes are created and the +% sample correlation is displayed. +% +% The push buttons at the bottom of the joint distribution plot select +% two different continuous distribution choices for a QCEFF for the +% observed variable, normal and RHF. The push buttons at the right +% select the continuous distribution for the unobserved variable +% and include normal, Gamma, an unbounded BNRH (RHF) and a bounded +% BNRH. +% +% Clicking on 'Update Ensemble' results in an updated joint ensemble +% along with updated marginals for both the observed and unobserved +% variables. The plot to the right is the joint distribution of the +% prior and posterior ensemble in the probit probability integral +% (PPI) transform space. +% +% The plot at the lower right shows more detail of the update for the +% observed variable including the likelihood in red. +% +% The mean and standard deviation of the likelihood can be +% changed in the red box. +% +% See also: bounded_oned_ensemble.m gaussian_product.m oned_cycle.m oned_ensemble.m +% oned_model.m oned_model_inf.m run_lorenz_63.m run_lorenz_96.m +% run_lorenz_96_inf.m twod_ensemble.m + +%% DART software - Copyright UCAR. This open source software is provided +% by UCAR, "as is", without charge, subject to all terms of use at +% http://www.image.ucar.edu/DAReS/DART/DART_download +% +% DART $Id$ + +help twod_ppi_ensemble; + +atts = stylesheet; % get the default fonts and color + +figureXmin = 0; +figureYmin = 0; +figureWidth = 1100; handles.figureWidth = figureWidth; +figureHeight = 800; handles.figureHeight = figureHeight; +handles.figureNormalization = 1000; + +%% Create a figure layout +handles.figure1 = figure('Name', 'twod_ppi_ensemble', ... + 'Position', [figureXmin figureYmin figureWidth figureHeight], ... + 'Color', atts.background); + +% Position has units of normalized and therefore the components must be +% Fractions of figure. Also, the actual text has FontUnits which are +% normalized, so the Font Size is a fraction of the text box/ edit box that +% contains it. Making the font size and box size normalized allows the text +% size and box size to change proportionately if the user resizes the +% figure window. This is true for all text/edit boxes and buttons. + +%% Create a button that lets you create a new ensemble. +handles.ui_button_create_ensemble = uicontrol('Style', 'pushbutton', ... + 'Units', 'Normalized', ... + 'Position', [gpx(0.2) gpy(0.75) gpx(0.2) gpy(0.04)] , ... + 'String', 'Create New Ensemble', ... + 'BackgroundColor', 'White', ... + 'Callback', @create_ensemble_Callback, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.38); + +%% Create a Button that updates the ensemble. +% This button is disabled at first because there is no ensemble to update. +handles.ui_button_update_ensemble = uicontrol('Style', 'pushbutton', ... + 'Units', 'Normalized', ... + 'Position', [gpx(0.42) gpy(0.75) gpx(0.2) gpy(0.04)] , ... + 'String', 'Update Ensemble', ... + 'BackgroundColor', 'White', ... + 'Callback', @update_ensemble_Callback, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.38, ... + 'Enable', 'Off'); + +%% Create Radio Button group with choices for the state variable PPI distribution + +handles.ui_assimilation_buttons = uibuttongroup('Visible','off', ... + 'Position', [gpx(0.01) gpy(0.55) gpx(0.15) gpy(0.11)] , ... + 'BorderType', 'none', ... + 'BackgroundColor', atts.background, ... + 'SelectionChangeFcn', @Assimilation_selection); + +set(handles.ui_assimilation_buttons, 'title', 'PPI DISTRIBUTION', 'fontsize', 16, ... + 'FontUnits', 'Normalized'); + +% Positions are relative to the Button Group +radioXmin = 0.05; +radioWidth = 0.9; +nudge = 0.02; +radioHeight = (1 - 5*nudge)/4; + +radioYmin = 1.0 - radioHeight - nudge; +handles.ui_radio_Normal = uicontrol(handles.ui_assimilation_buttons, ... + 'Style','radiobutton', ... + 'Units' , 'Normalized', ... + 'Position',[radioXmin radioYmin radioWidth radioHeight], ... + 'String','Normal', ... + 'BackgroundColor' , atts.background, ... + 'Foreground' , 'Black', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'Normalized' , ... + 'FontWeight', 'Bold' , ... + 'FontSize', 0.8, ... + 'HandleVisibility','off'); + +radioYmin = radioYmin - radioHeight - nudge; +handles.ui_radio_Gamma = uicontrol(handles.ui_assimilation_buttons, ... + 'Style','radiobutton', ... + 'String','Gamma', ... + 'Units' , 'Normalized', ... + 'Position',[radioXmin radioYmin radioWidth radioHeight], ... + 'BackgroundColor' , atts.background, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'Normalized' , ... + 'FontWeight', 'Bold' , ... + 'FontSize' , 0.8, ... + 'HandleVisibility','off'); + +radioYmin = radioYmin - radioHeight - nudge; +handles.ui_radio_RHF = uicontrol(handles.ui_assimilation_buttons, ... + 'Style','radiobutton', ... + 'String','RHF', ... + 'Units' , 'Normalized', ... + 'Position',[radioXmin radioYmin radioWidth radioHeight], ... + 'BackgroundColor' , atts.background, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'Normalized' , ... + 'FontWeight', 'Bold' , ... + 'FontSize' , 0.8, ... + 'HandleVisibility','off'); + +radioYmin = radioYmin - radioHeight - nudge; +handles.ui_radio_BNRH = uicontrol(handles.ui_assimilation_buttons, ... + 'Style','radiobutton', ... + 'String','BNRH', ... + 'Units' , 'Normalized', ... + 'Position',[radioXmin radioYmin radioWidth radioHeight], ... + 'BackgroundColor' , atts.background, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'Normalized' , ... + 'FontWeight', 'Bold' , ... + 'FontSize' , 0.8, ... + 'HandleVisibility','off'); + +% Make the uibuttongroup visible after creating child objects. +set(handles.ui_assimilation_buttons, 'Visible','On'); +selected = get(handles.ui_assimilation_buttons,'SelectedObject'); +handles.state_dist_type = get(selected,'String'); + +%% Create Radio Button group with choices for the obs variable PPI distribution + +handles.ui_obs_ppi_buttons = uibuttongroup('Visible','off', ... + 'Position', [gpx(0.33) gpy(0.1) gpx(0.16) gpy(0.11)] , ... + 'BorderType', 'none', ... + 'BackgroundColor', atts.background, ... + 'SelectionChangeFcn', @obs_ppi_selection); + +set(handles.ui_obs_ppi_buttons, 'title', 'PPI DISTRIBUTION', 'fontsize', 16, ... + 'FontUnits', 'Normalized'); + +% Positions are relative to the Button Group +radioXmin = 0.05; +radioWidth = 0.9; +nudge = 0.02; +radioHeight = (1 - 5*nudge)/4; + +radioYmin = 1.0 - radioHeight - nudge; +handles.ui_radio_obs_normal = uicontrol(handles.ui_obs_ppi_buttons, ... + 'Style','radiobutton', ... + 'Units' , 'Normalized', ... + 'Position',[radioXmin radioYmin radioWidth radioHeight], ... + 'String','Normal', ... + 'BackgroundColor' , atts.background, ... + 'Foreground' , 'Black', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'Normalized' , ... + 'FontWeight', 'Bold' , ... + 'FontSize', 0.8, ... + 'HandleVisibility','off'); + +radioYmin = radioYmin - radioHeight - nudge; +handles.ui_radio_obs_RHF = uicontrol(handles.ui_obs_ppi_buttons, ... + 'Style','radiobutton', ... + 'String','RHF', ... + 'Units' , 'Normalized', ... + 'Position',[radioXmin radioYmin radioWidth radioHeight], ... + 'BackgroundColor' , atts.background, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'Normalized' , ... + 'FontWeight', 'Bold' , ... + 'FontSize' , 0.8, ... + 'HandleVisibility','off'); + +% Make the uibuttongroup visible after creating child objects. +set(handles.ui_obs_ppi_buttons, 'Visible','On'); +selected = get(handles.ui_obs_ppi_buttons,'SelectedObject'); +handles.obs_dist_type = get(selected,'String'); + +%% Create a red panel with two text boxes and two edit boxes next to them + +handles.ObservationPanel = uipanel('BackgroundColor', atts.red, ... + 'BorderType','none', ... + 'Units', 'Normalized', ... + 'Position',[gpx(0.7) gpy(0.20) gpx(0.2), gpy(0.06)]); + +% Positions are relative to the panel +nudge = 0.05; +Xmin = 0.05; +dX1 = 0.65; +dX2 = 1.0 - dX1 - 2*nudge; +Height = (1.0 - 3*nudge)/2; + +Ymin = 1.0 - Height - nudge; +handles.ui_text_observation = uicontrol(handles.ObservationPanel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [Xmin Ymin-nudge dX1 Height] , ... + 'String', 'Observation' , ... + 'BackgroundColor', atts.red, ... + 'ForegroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontWeight', 'Bold', ... + 'FontSize', 0.50); + +handles.ui_edit_observation = uicontrol(handles.ObservationPanel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', ... + 'Position', [Xmin+dX1 Ymin dX2 Height] , ... + 'String', '5' , ... + 'Callback', @edit_observation_Callback, ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.50); + +Ymin = Ymin - Height - nudge; +handles.ui_text_obs_error_sd = uicontrol(handles.ObservationPanel, ... + 'Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [Xmin Ymin-nudge dX1 Height] , ... + 'String', 'Obs. Error SD' , ... + 'BackgroundColor', atts.red, ... + 'ForegroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontWeight', 'Bold', ... + 'FontSize', 0.50); + +handles.ui_edit_obs_error_sd = uicontrol(handles.ObservationPanel, ... + 'Style', 'edit', ... + 'Units', 'Normalized', .... + 'Position', [Xmin+dX1 Ymin dX2 Height] , ... + 'String', '1' , ... + 'Callback', @edit_obs_error_sd_Callback, ... + 'BackgroundColor', 'White', ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontSize', 0.50); + +%% Creates an error text inside the graph that is invisible at first but +% appears if a user enters an incorrect number for the observation or the +% observation error sd +handles.ui_text_error = uicontrol('Style', 'text', ... + 'Units', 'Normalized', ... + 'Position', [gpx(0.275), gpy(0.4), gpx(0.35), gpy(0.15)], ... + 'String', 'Error' , ... + 'BackgroundColor', 'white', ... + 'ForegroundColor', atts.red, ... + 'FontName', atts.fontname, ... + 'FontUnits', 'normalized', ... + 'FontWeight', 'Bold', ... + 'FontSize', 0.15, ... + 'Visible', 'Off'); + +%Set up all the handles variables and create the graphs +initialize(); + +%% Creates the graph of the observation and the observation error sd +% well as the ensemble estimates of the observation. + +h_observation = get(handles.ui_edit_observation); +h_obs_error_sd = get(handles.ui_edit_obs_error_sd); +observation = str2double(h_observation.String); +obs_error_sd = str2double(h_obs_error_sd.String); + +plotposition = [ gpx(0.20) gpy(0.35) gpx(0.07), gpy(0.35) ; ... + gpx(0.275) gpy(0.35) gpx(0.35) gpy(0.35); ... + gpx(0.275) gpy(0.275) gpx(0.35) gpy(0.07); + gpx(0.7) gpy(0.35) gpx(0.35) gpy(0.35)]; + + +%% Create the axis for the unobserved marginal graphic + +% Set parameters for the initial x and y value ranges +handles.xmin = -10; +handles.xmax = 10; +handles.ymin = -2; +handles.ymax = 10; + +handles.h_unobMarginal = axes('Position', plotposition(1,:), ... + 'FontName', atts.fontname, ... + 'FontSize', atts.fontsize); +axis([0 1 handles.ymin handles.ymax]); +set(gca, 'FontUnits', 'Normalized'); +ylabel('Unobserved State Variable', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); +hold on + +%% Get a subplot for the joint prior + +handles.h_joint = axes('Position', plotposition(2,:), ... + 'FontName', atts.fontname, ... + 'FontSize', atts.fontsize); + +axis([handles.xmin, handles.xmax, handles.ymin, handles.ymax]); +grid on +hold on +title('Joint Distribution', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); +set(gca, 'color', [0.8, 0.8, 0.8]) +set(gca, 'FontUnits', 'Normalized'); + +% Plot long axes at x and y = 0 +ax = [-1000, 1000]; +ay = [0, 0]; +plot(ax, ay, 'k', 'linewidth', 2); +ax = [0 0]; +ay = [-1000 1000]; +plot(ax, ay, 'k', 'linewidth', 2); + +h_click = text(0, 8, {'Click inside Joint Distribution plot','grey area to create member.'}, ... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); +h_finish = text(0, 6, 'Click outside of grey to finish.',... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); +h_err_text = text(0, 4, 'An ensemble must have at least 2 members.', 'Color', atts.red, ... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); + +handles.h_correl = text(5, 9, ' ', 'Color', atts.green, 'FontWeight', 'Bold', ... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'Visible', 'off'); +% Prepare error message for negative y in case it's needed +h_neg_y = text(0, 3, 'Unobserved variable must be nonnegative', ... + 'FontSize', 16, 'HorizontalAlignment','center', 'FontUnits', 'Normalized', 'visible', 'off'); + + +%% Get a subplot for the ppi joint prior + +handles.h_ppi_joint = axes('Position', plotposition(4,:), ... + 'FontName', atts.fontname, ... + 'FontSize', atts.fontsize); +axis([-3 3 -3 3]); +grid on +hold on +title('Joint PPI SpaceDistribution', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); +set(gca, 'FontUnits', 'Normalized'); + +xlabel('PPI Transformed Observed', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); +ylabel('PPI Transformed Unobserved', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); + + +%% Create a subplot for the observed variable marginal + +handles.h_obMarginal = axes('Position', plotposition(3,:), ... + 'FontName', atts.fontname, 'FontSize', atts.fontsize); + +% May want to mess with vertical axis for prior density +axis([handles.xmin handles.xmax 0 1]); +hold on +xlabel('Observed Quantity', 'Fontsize', atts.fontsize, 'FontUnits', 'Normalized'); +set(gca, 'FontUnits', 'Normalized'); +handles.h_obs_marg = plot(observation, 0, '*', 'MarkerSize', 16, 'LineWidth', 2.0); +set(handles.h_obs_marg,'Color',atts.red) + +% Turn off the unwanted tick labels +set(handles.h_unobMarginal, 'Xticklabel', []); +set(handles.h_joint , 'Xticklabel', [], 'Yticklabel', []); +set(handles.h_obMarginal , 'Yticklabel', []); + +%This graph is the graph of the observation +handles.h_obs_likelihood = axes( ... + 'Position', [gpx(0.63) gpy(0.035), gpx(0.35), gpy(0.125)], ... + 'FontName', atts.fontname, ... + 'FontSize', atts.fontsize); + +handles.h_marg_obs_plot = plot_gaussian(observation, obs_error_sd, 1); + +axis([handles.xmin handles.xmax -handles.y_max/5 1]); + +% Set the ticks +set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]); +set(handles.h_marg_obs_plot, 'Color', atts.red, ... + 'Linestyle', '--', ... + 'LineWidth', 2); +set(gca, 'FontUnits', 'Normalized'); +xlabel('Observed Quantity', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); +ylabel('Observation Likelihood', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); +title('Marginal Distribution of Observation', 'FontSize', atts.fontsize, 'FontUnits', 'Normalized'); +hold on + +% Plot an asterisk for the observed value +handles.h_obs_ast = plot(observation, 0, 'r*', 'MarkerSize', 16, 'LineWidth', 2.0); +set(handles.h_obs_ast,'Color',atts.red) + +% Plot an axis; +plot([-1000 1000], [0 0], 'k', 'LineWidth', 2); + +%% ---------------------------------------------------------------------------- + + function create_ensemble_Callback(~,~) + % Allows the user to create a new ensemble in the left axes + + % Disable the update ensemble button and all other active buttons + set(handles.ui_button_create_ensemble, 'Enable', 'Off'); + set(handles.ui_button_update_ensemble, 'Enable', 'Off'); + set(handles.ui_edit_observation, 'Enable', 'Off'); + set(handles.ui_edit_obs_error_sd, 'Enable', 'Off'); + + % Clear out any old ensemble members if they exist + axes(handles.h_joint) + for i = 1:handles.ens_size + set(handles.h_ens_member(i), 'Visible', 'off'); + set(handles.h_gui_marg(i), 'Visible', 'off'); + set(handles.h_unobs(i), 'Visible', 'off'); + set(handles.h_marg(i), 'Visible', 'off'); + end + + % Turn off any posterior old plotting + set(handles.h_update_ens, 'Visible', 'off'); + set(handles.h_marg_update, 'Visible', 'off'); + set(handles.h_marg_inc, 'Visible', 'off'); + set(handles.h_marg_state, 'Visible', 'off'); + set(handles.h_state_inc, 'Visible', 'off'); + set(handles.h_joint_update, 'Visible', 'off'); + set(handles.h_joint_inc, 'Visible', 'off'); + set(handles.h_ppi_prior, 'Visible', 'off'); + set(handles.h_ppi_post, 'Visible', 'off'); + set(handles.h_joint_ppi, 'Visible', 'off'); + + % Clear out the old best fit line + set(handles.h_best_fit, 'Visible', 'off'); + set(handles.h_correl, 'Visible', 'off'); + + % Turn on the data entry messages + set(h_click, 'Visible', 'on'); + set(h_finish, 'Visible', 'on'); + set(h_err_text, 'Visible', 'on'); + + % Work in the joint distribution plot + axes(handles.h_joint); + hold on + + % Need to guarantee at least 2 ensemble members + ens_size = 0; + + while ens_size < 1000 + [xt, yt] = ginput(1); + gca; + + % Make sure negative y message is off + set(h_neg_y, 'visible', 'off'); + + % Make sure that the click was in the correct set of axes + % Terminate by clicking outside of graph range + if(xt < handles.xmin || xt > handles.xmax || yt < handles.ymin || ... + yt > handles.ymax || gca ~= handles.h_joint) + axes(handles.h_joint); %#ok + break; + elseif(yt<0 & yt >= handles.ymin) + % Still in grey shaded part of plot but value of y is negative + set(h_neg_y, 'visible', 'on'); + else + + ens_size = ens_size + 1; + x(1, ens_size) = xt; %#ok + x(2, ens_size) = yt; %#ok + + axes(handles.h_joint); %#ok + handles.h_ens_member(ens_size) = ... + plot(x(1, ens_size), x(2, ens_size), '*', ... + 'MarkerSize', 16, 'Color', atts.green, 'LineWidth',2.0); + + % Plot the marginal for the unobserved state variable + %>@ TODO POSSIBLE IMPROVEMENT ... annotate new marginal mean, sd + axes(handles.h_unobMarginal); %#ok + handles.h_unobs(ens_size) = ... + plot(0, x(2, ens_size), '*', 'MarkerSize', 16, 'Color', atts.green, 'LineWidth',2.0); + + % Plot the marginal for the observed quantity + axes(handles.h_obMarginal); %#ok + handles.h_marg(ens_size) = ... + plot(x(1, ens_size), 0, '*', 'MarkerSize', 16, 'Color', atts.green, 'LineWidth',2.0); + + % Plot the marginal in the gui frame + axes(handles.h_obs_likelihood); %#ok + handles.h_gui_marg(ens_size) = ... + plot(x(1, ens_size), 0, '*', 'MarkerSize', 16, 'Color', atts.green, 'LineWidth',2.0); + + % Then switch back to axes(handles.h_joint) + axes(handles.h_joint); %#ok + + if (ens_size < 2) + continue + end + + % Clear out the error message if it's been made visible + set(h_err_text, 'Visible', 'off'); + set(h_click, 'Visible', 'off'); + + prior_correl = corrcoef(x(1, :), x(2, :)); + str1 = sprintf('Correlation = %f', prior_correl(1,2)); + set(handles.h_correl,'String', str1, 'Visible', 'on') + + end + end + + % it is possible that they click outside the box before completing a viable + % ensemble ... in this case, just return and let them start over. + if (ens_size > 0) + + % Turn off the data entry messages + set(h_finish, 'Visible', 'off'); + + %% Ensemble created, compute mean and sd, clean up and return + % Set the global gui storage + handles.ens_size = ens_size; + handles.ens_members = x; + + % Plot the best fit line on the ensemble + prior_mean = mean(x, 2); + prior_cov = cov(x(1, :), x(2, :)); + slope = prior_cov(1, 2) / var(x(1, :)); + intercept = prior_mean(2) - slope * prior_mean(1); + + best_x = [-1000 1000]; + best_y = slope * best_x + intercept; + + handles.h_best_fit = plot(best_x, best_y, 'g', 'LineWidth', 2.0); + set(handles.h_best_fit, 'Color', atts.green); + + end + + % Enable the update ensemble button + set(handles.ui_button_create_ensemble, 'Enable', 'On'); + set(handles.ui_button_update_ensemble, 'Enable', 'On'); + set(handles.ui_edit_observation, 'Enable', 'On'); + set(handles.ui_edit_obs_error_sd, 'Enable', 'On'); + + % Reset focus to the menu gui window + axes(handles.h_obs_likelihood); + + end + +%% ------------------------------------------------------------------------- + + function update_ensemble_Callback(~,~) + % Uses the assimilation to update the ensemble according to the + % observation, and then plots it on the main graph on the left, the + % two marginals, and the right observation graph + + axes(handles.h_obs_likelihood); + % Turn off any old points + set(handles.h_update_ens, 'Visible', 'off'); + set(handles.h_marg_update, 'Visible', 'off'); + set(handles.h_marg_inc, 'Visible', 'off'); + set(handles.h_marg_state, 'Visible', 'off'); + set(handles.h_state_inc, 'Visible', 'off'); + set(handles.h_joint_update, 'Visible', 'off'); + set(handles.h_joint_inc, 'Visible', 'off'); + set(handles.h_ppi_prior, 'Visible', 'off'); + set(handles.h_ppi_post, 'Visible', 'off'); + set(handles.h_joint_ppi, 'Visible', 'off'); + + prior_obs = handles.ens_members(1, :); + prior_state = handles.ens_members(2, :); + h_observation = get(handles.ui_edit_observation); + h_obs_error_sd = get(handles.ui_edit_obs_error_sd); + observation = str2double(h_observation.String); + obs_error_sd = str2double(h_obs_error_sd.String); + + %If ensemble is not empty + if (size(prior_obs,2) > 0) + + % Observation space increment method is controlled by obs PPI distribution choide + if(strcmp(handles.obs_dist_type, 'Normal')) + % The observation space updates are standard EAKF + [obs_increments, ~] = ... + obs_increment_eakf(prior_obs, observation, obs_error_sd^2); + else + % Observation space updates by unbounded RHF + [obs_increments, ~] = ... + obs_increment_rhf(prior_obs, observation, obs_error_sd^2, false); + end + + % Add on increments to get posterior + post_obs = prior_obs + obs_increments; + + % Updating the state variables + [post_state, prior_obs_ppi, post_obs_ppi, prior_state_ppi, post_state_ppi] = ... + ppi_update(prior_obs, prior_state, post_obs, handles.state_dist_type, ... + handles.obs_dist_type); + + %Set the y-coordinate of the ensembles, to be halfway between 0 and + %the bottom of the graph; + y(1:handles.ens_size) = -handles.y_max/10; + + handles.h_update_ens = plot(post_obs, y, '*', 'MarkerSize', 16, 'Color', atts.blue); + + % Plot the PPI priors and posteriors + axes(handles.h_ppi_joint); + handles.h_ppi_prior = ... + plot(prior_obs_ppi, prior_state_ppi, '*', 'Markersize', 16, 'Color', atts.green); + handles.h_ppi_post = ... + plot(post_obs_ppi, post_state_ppi, '*', 'Markersize', 16, 'Color', atts.blue); + % Plot the connecting segments + for i = 1:handles.ens_size + cx = [prior_obs_ppi(i), post_obs_ppi(i)]; + cy = [prior_state_ppi(i), post_state_ppi(i)]; + handles.h_joint_ppi(i) = plot(cx, cy, 'c'); + end + + % Plot the increments in the state marginal plot + axes(handles.h_obMarginal); + + % Need to sort ensemble to get nice ordering for increments + [~, sort_obs_ind] = sort(prior_obs); + for i = 1:handles.ens_size + y(i) = i / (handles.ens_size + 1); + handles.h_marg_update(i) = ... + plot(post_obs(sort_obs_ind(i)), y(i), '*', 'MarkerSize', 16, 'Color', atts.blue); + % Also plot a segment in blue + handles.h_marg_inc(i) = ... + plot([prior_obs(sort_obs_ind(i)), post_obs(sort_obs_ind(i))], ... + [y(i), y(i)], 'c'); + end + + axes(handles.h_unobMarginal); + + % Now need to sort the state variable ensemble to get nice ordering + [~, sort_ind] = sort(prior_state); + for i = 1:handles.ens_size + handles.h_marg_state(i) = ... + plot(y(i), post_state(sort_ind(i)), '*', 'MarkerSize', 16, 'Color', atts.blue); + % Also plot a segment in blue + handles.h_state_inc(i) = plot([y(i), y(i)], ... + [prior_state(sort_ind(i)), post_state(sort_ind(i))], 'c'); + end + + % Plot the updated joint distribution points + axes(handles.h_joint); + for i = 1:handles.ens_size + handles.h_joint_update(i) = plot(post_obs(i), post_state(i), ... + '*', 'MarkerSize', 16, 'Color', atts.blue); + handles.h_joint_inc(i) = plot([prior_obs(i), post_obs(i)], ... + [prior_state(i), post_state(i)], 'c'); + end + + % Return the focus to the window with pushbuttons + axes(handles.h_obs_likelihood); + + end + end + +%% ------------------------------------------------------------------------- + + function edit_observation_Callback(~, ~) + + % Enable things that an error might have turned off + set(handles.ui_edit_obs_error_sd, 'Enable', 'on'); + set(handles.ui_button_create_ensemble, 'Enable', 'on'); + + % Only enable the update ensemble pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_update_ensemble, 'Enable', 'on'); + end + + % Get the value of the observation + if( isfinite( str2double( get(handles.ui_edit_observation, 'String')))) + observation = str2double(get(handles.ui_edit_observation, 'String')); + + if (observation > handles.xmax) + set(handles.ui_edit_observation, 'String', ['< ', num2str(handles.xmax)]); + input_error('observation'); + return; + elseif (observation < handles.xmin) + set(handles.ui_edit_observation, 'String', ['> ', num2str(handles.xmin)]); + input_error('observation'); + return; + end + + %Set background color to normal and error text off + set(handles.ui_edit_observation, 'BackgroundColor', 'white'); + set(handles.ui_text_error, 'Visible', 'Off'); + set(handles.h_ens_member, 'Visible', 'On'); + set(handles.h_best_fit, 'Visible', 'On'); + set(handles.h_joint_update, 'Visible', 'On'); + set(handles.h_joint_inc, 'Visible', 'On'); + set(handles.h_ppi_prior, 'Visible', 'On'); + set(handles.h_ppi_post, 'Visible', 'On'); + set(handles.h_joint_ppi, 'Visible', 'On'); + else + set(handles.ui_edit_observation, 'String', '?'); + input_error('observation'); + return + end + + % Get the value of the observation error sd + h_obs_error_sd = get(handles.ui_edit_obs_error_sd); + obs_error_sd = str2double(h_obs_error_sd.String); + + % Plot the updated distribution + set(handles.h_marg_obs_plot, 'Visible', 'off'); + handles.h_marg_obs_plot = plot_gaussian(observation, obs_error_sd, 1); + set(handles.h_marg_obs_plot, 'Color', atts.red, 'Linestyle', '--', 'Linewidth', 2); + + % Update the observation asterisk + set(handles.h_obs_ast, 'Visible', 'off'); + handles.h_obs_ast = plot(observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0); + set(handles.h_obs_ast,'Color',atts.red) + + % Plot the updated obs distribution on the marginal subplot + axes(handles.h_obMarginal); + + % Plot the updated observation in the marginal + set(handles.h_obs_marg, 'Visible', 'off'); + handles.h_obs_marg = plot(observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0); + set(handles.h_obs_marg,'Color',atts.red) + + % Replot the update ensemble members so the correlate to new + % observation + update_ensemble_Callback(); + + axes(handles.h_obs_likelihood); + + end + +%% -------------------------------------------------------------------- + + function edit_obs_error_sd_Callback(~, ~) + + % Enable things that an error might have turned off + set(handles.ui_edit_observation, 'Enable', 'on') + set(handles.ui_button_create_ensemble, 'Enable', 'on') + + % Only enable the update ensemble pushbutton if an ensemble has been created + if(handles.ens_size > 0) + set(handles.ui_button_update_ensemble, 'Enable', 'on'); + end + + % Get the value of the observation error standard deviation + if(isfinite(str2double(get(handles.ui_edit_obs_error_sd, 'String'))) && ... + str2double(get(handles.ui_edit_obs_error_sd, 'String')) > 0) + obs_error_sd = str2double(get(handles.ui_edit_obs_error_sd, 'String')); + + %Set background color to normal and error text off + set(handles.ui_edit_obs_error_sd, 'BackgroundColor', 'white'); + set(handles.ui_text_error, 'Visible', 'Off'); + set(handles.h_ens_member, 'Visible', 'On'); + set(handles.h_best_fit, 'Visible', 'On'); + set(handles.h_joint_update, 'Visible', 'On'); + set(handles.h_joint_inc, 'Visible', 'On'); + set(handles.h_ppi_prior, 'Visible', 'On'); + set(handles.h_ppi_post, 'Visible', 'On'); + set(handles.h_joint_ppi, 'Visible', 'On'); + + else + set(handles.ui_edit_obs_error_sd, 'String', '?'); + input_error('standard deviation'); + return + end + + % Get the value of the observation + h_observation = get(handles.ui_edit_observation); + observation = str2double(h_observation.String); + handles.y_max = norm_pdf(observation, observation, obs_error_sd); + + %Give 0.2 cushion to y_max + handles.y_max = handles.y_max + 0.2; + + %Update the axis based on the new y_max + axis([handles.xmin handles.xmax -handles.y_max/5 handles.y_max]); + + set(gca,'YTickMode','auto') + ticks = get(gca,'YTick'); + inds = (ticks >= 0); % Only show ticks for values greater than 0 + newticks = ticks(inds); + set(gca,'YTick',newticks) + + % Replot the update ensemble members so the correlate to new obs_sd + update_ensemble_Callback(); + + % Plot the updated distribution on the menu plot + + set(handles.h_marg_obs_plot, 'Visible', 'off'); + handles.h_marg_obs_plot = plot_gaussian(observation, obs_error_sd, 1); + set(handles.h_marg_obs_plot, 'Color', atts.red, 'Linestyle', '--', 'Linewidth', 2); + + % Update the observation asterisk + + set(handles.h_obs_ast, 'Visible', 'off'); + handles.h_obs_ast = plot(observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0); + set(handles.h_obs_ast, 'Color', atts.red) + + % Plot the updated observation in the marginal + axes(handles.h_obMarginal); + + set(handles.h_obs_marg, 'Visible', 'off'); + handles.h_obs_marg = plot(observation, 0, 'r*', 'MarkerSize', 16,'LineWidth',2.0); + set(handles.h_obs_marg, 'Color', atts.red) + + % Reset focus to the menu gui window + axes(handles.h_obs_likelihood); + + end + +%% ---------------------------------------------------------------------------- + function initialize() + % Insert the ensemble structure into this + handles.ens_size = 0; + handles.ens_members = []; + handles.h_update_ens = []; + handles.h_ens_member = []; + handles.h_best_fit = []; + handles.h_marg_obs_plot = []; + handles.h_obs_ast = []; + handles.h_obs_marg = []; + handles.h_gui_marg = []; + handles.h_unobs = []; + handles.h_marg = []; + handles.h_marg_update = []; + handles.h_marg_inc = []; + handles.h_marg_state = []; + handles.h_state_inc = []; + handles.h_joint_update = []; + handles.h_ppi_prior = []; + handles.h_ppi_post = []; + handles.h_joint_ppi = []; + handles.h_joint_inc = []; + handles.h_correl = []; + handles.first_correl = true; + handles.y_max = 1; + end + +%% ---------------------------------------------------------------------------- + + function Assimilation_selection(~, eventdata) + % Function is called whenever a radio button has been selected, it sets + % the global filter variable + + % eventdata refers to the data in the GUI when a radio button in the + % group is changed + + % Set the filter_type string to newest radiobutton Value + + handles.state_dist_type = get(eventdata.NewValue,'String'); + end +%% ---------------------------------------------------------------------------- + + function obs_ppi_selection(~, eventdata) + % Function is called whenever a radio button has been selected, it sets + % the global filter variable + + % eventdata refers to the data in the GUI when a radio button in the + % group is changed + + % Set the filter_type string to newest radiobutton Value + + handles.obs_dist_type = get(eventdata.NewValue,'String'); + end + +%% ---------------------------------------------------------------------------- + + function input_error(mystring) + switch(lower(mystring)) + case 'observation' + set(handles.ui_edit_observation, 'BackgroundColor', atts.red); + set(handles.h_ens_member, 'Visible', 'Off'); + set(handles.h_best_fit, 'Visible', 'Off'); + set(handles.h_joint_update, 'Visible', 'Off'); + set(handles.h_joint_inc, 'Visible', 'Off'); + set(handles.h_ppi_prior, 'Visible', 'Off'); + set(handles.h_ppi_post, 'Visible', 'Off'); + set(handles.h_joint_ppi, 'Visible', 'Off'); + msg_string = ['Observation must be a number between ', ... + num2str(handles.xmin), ' and ', num2str(handles.xmax)]; + set(handles.ui_text_error, 'String' , msg_string); + set(handles.ui_text_error, 'Visible', 'On'); + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_obs_error_sd, 'Enable', 'off'); + set(handles.ui_button_create_ensemble, 'Enable', 'off'); + set(handles.ui_button_update_ensemble, 'Enable', 'off'); + + otherwise + + set(handles.ui_edit_obs_error_sd, 'BackgroundColor', atts.red); + set(handles.h_ens_member, 'Visible', 'Off'); + set(handles.h_best_fit, 'Visible', 'Off'); + set(handles.h_joint_update, 'Visible', 'Off'); + set(handles.h_joint_inc, 'Visible', 'Off'); + set(handles.h_ppi_prior, 'Visible', 'Off'); + set(handles.h_ppi_post, 'Visible', 'Off'); + set(handles.h_joint_ppi, 'Visible', 'Off'); + set(handles.ui_text_error, 'String' , 'Observation Error SD must be a number greater than 0'); + set(handles.ui_text_error, 'Visible', 'On'); + + % Disable other input to guarantee only one error at a time! + set(handles.ui_edit_observation, 'Enable', 'off') + set(handles.ui_button_create_ensemble, 'Enable', 'off') + set(handles.ui_button_update_ensemble, 'Enable', 'off') + end % of switch + end % function input_error + +%% ---------------------------------------------------------------------------- + + function [x] = gpx(f) + + x = f * handles.figureNormalization / handles.figureWidth; + + end + +%% ---------------------------------------------------------------------------- + + function [y] = gpy(f) + + y = f * handles.figureNormalization / handles.figureHeight; + + end + +end + + +% +% $URL$ +% $Revision$ +% $Date$ diff --git a/guide/DART_LAB/presentation/DART_LAB_Section01.pdf b/guide/DART_LAB/presentation/DART_LAB_Section01.pdf index cc309eff84..e4b697365c 100644 Binary files a/guide/DART_LAB/presentation/DART_LAB_Section01.pdf and b/guide/DART_LAB/presentation/DART_LAB_Section01.pdf differ diff --git a/guide/DART_LAB/presentation/DART_LAB_Section02.pdf b/guide/DART_LAB/presentation/DART_LAB_Section02.pdf index 9a492d5f2e..5c5531bac9 100644 Binary files a/guide/DART_LAB/presentation/DART_LAB_Section02.pdf and b/guide/DART_LAB/presentation/DART_LAB_Section02.pdf differ diff --git a/guide/DART_LAB/presentation/DART_LAB_Section03.pdf b/guide/DART_LAB/presentation/DART_LAB_Section03.pdf index 0f9410ad54..9794f6466c 100644 Binary files a/guide/DART_LAB/presentation/DART_LAB_Section03.pdf and b/guide/DART_LAB/presentation/DART_LAB_Section03.pdf differ diff --git a/guide/DART_LAB/presentation/DART_LAB_Section04.pdf b/guide/DART_LAB/presentation/DART_LAB_Section04.pdf index 7efce32869..ebad45d1cc 100644 Binary files a/guide/DART_LAB/presentation/DART_LAB_Section04.pdf and b/guide/DART_LAB/presentation/DART_LAB_Section04.pdf differ diff --git a/guide/DART_LAB/presentation/DART_LAB_Section05.pdf b/guide/DART_LAB/presentation/DART_LAB_Section05.pdf index 28b8f5df63..4b6c56f5d5 100644 Binary files a/guide/DART_LAB/presentation/DART_LAB_Section05.pdf and b/guide/DART_LAB/presentation/DART_LAB_Section05.pdf differ