[Dart-dev] [4011] DART/trunk: GUI presentation of fundamental assimilation concepts

nancy at ucar.edu nancy at ucar.edu
Tue Aug 25 17:04:35 MDT 2009


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090825/07a8a1fd/attachment-0001.html 
-------------- next part --------------
Added: DART/trunk/DART_LAB/matlab/advance_oned.m
===================================================================
--- DART/trunk/DART_LAB/matlab/advance_oned.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/advance_oned.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,25 @@
+function x_new = advance_oned(x, alpha)
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+x_new = x + comp_dt(x, alpha);
+end
+
+%---------------------------------------------------
+
+% Internal function comp_dt
+function dx = comp_dt(x, alpha)
+
+% Compute the time tendency; alpha controls nonlinearity
+dx = x + alpha .* x .* abs(x);
+
+end


Property changes on: DART/trunk/DART_LAB/matlab/advance_oned.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/comp_cov_factor.m
===================================================================
--- DART/trunk/DART_LAB/matlab/comp_cov_factor.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/comp_cov_factor.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,26 @@
+function cov_factor = comp_cov_factor(z_in, c)
+% Gaspari Cohn cutoff, z_in is the distance while c is the cutoff
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+z = abs(z_in);
+
+if( z >= c*2.0)
+   cov_factor = 0;
+elseif( z <= c )
+   r = z / c;
+   cov_factor = ((( -0.25*r +0.5)*r +0.625)*r -5.0/3.0)*r^2 + 1.0;
+else
+   r = z / c;
+   cov_factor = ((((r/12 -0.5)*r +0.625)*r +5.0/3.0)*r -5.0)*r + 4.0 - 2.0 / (3.0 * r);
+end
+


Property changes on: DART/trunk/DART_LAB/matlab/comp_cov_factor.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/g_prod_plot.m
===================================================================
--- DART/trunk/DART_LAB/matlab/g_prod_plot.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/g_prod_plot.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,107 @@
+function [prior_mean, prior_sd, obs_mean, obs_err_sd, is_err] = g_prod_plot(h)
+% Updates the plot of the prior and observation gaussians
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Successful return as defaul  t
+is_err = false;
+
+% Default failed returns for other quantities
+prior_mean = 0; prior_sd   = -1;
+obs_mean   = 0; obs_err_sd = -1;
+
+h_prior_mean = get(h.edit1);
+prior_mean = str2double(h_prior_mean.String);
+% The mean must be a number
+if(isnan(prior_mean)) 
+   error_banner(h, 'Prior Mean must be a number');
+   is_err = true;
+   return
+end
+
+h_prior_sd = get(h.edit2);
+prior_sd = str2double(h_prior_sd.String);
+
+% Prior sd must be a number
+if(isnan(prior_sd))
+   error_banner(h, 'Prior SD must be a postive number');
+   is_err = true;
+   return
+end
+
+% Prior sd must also be positive 
+if(prior_sd <= 0)
+   error_banner(h, 'Prior SD must be positive')
+   is_err = true;
+   return
+end
+   
+hold off
+prior_handle = plot_gaussian(prior_mean, prior_sd, 1);
+set(prior_handle, 'Color', [0 0.73 0], 'LineWidth', 2);
+hold on
+
+h_obs_mean = get(h.edit3);
+obs_mean = str2double(h_obs_mean.String);
+
+% Obs value must be a number
+if(isnan(obs_mean)) 
+   error_banner(h, 'Obs value must be a number');
+   is_err = true;
+   return
+end
+
+h_obs_err_sd = get(h.edit4);
+obs_err_sd = str2double(h_obs_err_sd.String);
+
+% Obs error sd must be a positive number
+if(isnan(obs_err_sd))
+   error_banner(h, 'Obs Error SD must be a positive number');
+   is_err = true;
+   return
+end
+
+if(obs_err_sd <= 0)
+   error_banner(h, 'Obs Error SD must be positive');
+   is_err = true;
+   return
+end
+
+obs_handle = plot_gaussian(obs_mean, obs_err_sd, 1);
+set(obs_handle, 'Color', 'r', 'LineStyle', '--', 'LineWidth', 2);
+
+% Put on a legend
+legend('Prior', 'Obs. Likelihood');
+
+end
+
+%---------------------------------------------------------
+
+% Internal function to write error banner
+function error_banner(h, message)
+
+   hold off
+   x= [1 2];
+   plot(x, 'Visible', 'off')
+   h_fig = get(h.figure1);
+   x_limits = get(h_fig.CurrentAxes, 'Xlim');
+   y_limits = get(h_fig.CurrentAxes, 'Ylim');
+   text(x_limits(1) * 7/8 + x_limits(2) / 8, mean(y_limits), ...
+      message, 'fontsize', 16, 'Color', 'r');
+
+   % While we're at it, clear out the values for posterior, too
+   set(h.text7, 'String', 'Posterior Mean = ');
+   set(h.text8, 'String', 'Posterior SD = ');
+   set(h.text9, 'String', 'Weight = ');
+   return; 
+
+end


Property changes on: DART/trunk/DART_LAB/matlab/g_prod_plot.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/gaussian_product.fig
===================================================================
(Binary files differ)


Property changes on: DART/trunk/DART_LAB/matlab/gaussian_product.fig
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + application/octet-stream

Added: DART/trunk/DART_LAB/matlab/gaussian_product.m
===================================================================
--- DART/trunk/DART_LAB/matlab/gaussian_product.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/gaussian_product.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,214 @@
+function varargout = gaussian_product(varargin)
+% GAUSSIAN_PRODUCT demonstrates the product of two gaussian distributions.
+%
+%      This is fundamental to ensemble data assimilation. Change the
+%      parameters of the gaussian for the Prior (green) and the Observation
+%      (red) and click on 'Plot Posterior'. 
+%
+%      The product of the two gaussians is a gaussian - in this case, 
+%      the 'Posterior'. If the parameters of the two gaussians are known, 
+%      the parameters of the resulting gaussian can be calculated.
+%
+% See also: oned_model, oned_ensemble, twod_ensemble, run_lorenz_63, 
+%           run_lorenz_96
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Last Modified by GUIDE v2.5 21-Mar-2009 22:11:31
+
+% Begin initialization code - DO NOT EDIT
+gui_Singleton = 1;
+gui_State = struct('gui_Name',       mfilename, ...
+                   'gui_Singleton',  gui_Singleton, ...
+                   'gui_OpeningFcn', @gaussian_product_OpeningFcn, ...
+                   'gui_OutputFcn',  @gaussian_product_OutputFcn, ...
+                   'gui_LayoutFcn',  [] , ...
+                   'gui_Callback',   []);
+if nargin && ischar(varargin{1})
+    gui_State.gui_Callback = str2func(varargin{1});
+end
+
+if nargout
+    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
+else
+    gui_mainfcn(gui_State, varargin{:});
+end
+% End initialization code - DO NOT EDIT
+
+
+% --- Executes just before gaussian_product is made visible.
+function gaussian_product_OpeningFcn(hObject, eventdata, handles, varargin)
+% This function has no output args, see OutputFcn.
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+% varargin   command line arguments to gaussian_product (see VARARGIN)
+
+help gaussian_product
+
+% Choose default command line output for gaussian_product
+handles.output = hObject;
+
+% Update handles structure
+guidata(hObject, handles);
+
+% Plot the initial prior and observation likelihood pdf's
+h = guihandles;
+g_prod_plot(h);
+
+% UIWAIT makes gaussian_product wait for user response (see UIRESUME)
+% uiwait(handles.figure1);
+
+
+% --- Outputs from this function are returned to the command line.
+function varargout = gaussian_product_OutputFcn(hObject, eventdata, handles) 
+% varargout  cell array for returning output args (see VARARGOUT);
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Get default command line output from handles structure
+varargout{1} = handles.output;
+
+
+
+function edit1_Callback(hObject, eventdata, handles)
+% hObject    handle to edit1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit1 as text
+%        str2double(get(hObject,'String')) returns contents of edit1 as a double
+
+g_prod_plot(handles);
+
+
+% --- Executes during object creation, after setting all properties.
+function edit1_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+
+function edit2_Callback(hObject, eventdata, handles)
+% hObject    handle to edit2 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit2 as text
+%        str2double(get(hObject,'String')) returns contents of edit2 as a double
+
+g_prod_plot(handles);
+
+
+% --- Executes during object creation, after setting all properties.
+function edit2_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit2 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+
+
+
+function edit3_Callback(hObject, eventdata, handles)
+% hObject    handle to edit3 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit3 as text
+%        str2double(get(hObject,'String')) returns contents of edit3 as a double
+
+g_prod_plot(handles);
+
+
+% --- Executes during object creation, after setting all properties.
+function edit3_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit3 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+
+
+
+function edit4_Callback(hObject, eventdata, handles)
+% hObject    handle to edit4 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit4 as text
+%        str2double(get(hObject,'String')) returns contents of edit4 as a double
+
+g_prod_plot(handles);
+
+
+% --- Executes during object creation, after setting all properties.
+function edit4_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit4 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+
+% --- Executes on button press in pushbutton1.
+function pushbutton1_Callback(hObject, eventdata, handles)
+% hObject    handle to pushbutton1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Need to replot prior and obs, then compute posterior and plot
+[prior_mean, prior_sd, obs_mean, obs_err_sd, is_err] = g_prod_plot(handles);
+
+% If there is an error, don't try to do posterior computation
+% But do zero out the posterior text values
+if(is_err) return; 
+   set(handles.text7, 'String', strcat('Posterior Mean = '));
+   set(handles.text8, 'String', strcat('Posterior SD = '));
+   set(handles.text9, 'String', strcat('Weight = '));
+end
+
+% Compute the posterior mean, sd and weight
+[post_mean, post_sd, weight] = ...
+    product_of_gaussians(prior_mean, prior_sd, obs_mean, obs_err_sd);
+post_handle = plot_gaussian(post_mean, post_sd, 1);
+set(post_handle, 'Color', 'b', 'LineWidth', 2);
+
+% Print values
+set(handles.text7, 'String', ['Posterior Mean = ', num2str(post_mean)]);
+set(handles.text8, 'String', ['Posterior SD = ', num2str(post_sd)]);
+
+% Also plot the weighted posterior as dashed
+post_handle = plot_gaussian(post_mean, post_sd, weight);
+set(post_handle, 'Color', 'b', 'Linestyle', '--');
+set(handles.text9, 'String', ['Weight = ', num2str(weight)]);
+
+legend('Prior', 'Obs. Likelihood', 'Posterior', 'Weighted Posterior');
+
+


Property changes on: DART/trunk/DART_LAB/matlab/gaussian_product.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/get_ens_rank.m
===================================================================
--- DART/trunk/DART_LAB/matlab/get_ens_rank.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/get_ens_rank.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,25 @@
+function [rank] = get_ens_rank(ens, x)
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+s_ens = sort(ens);
+
+for i = 1 : size(ens, 2)
+   if(x < s_ens(1, i))
+      rank = i;
+      return;
+   end
+end
+
+rank = size(ens, 2) + 1;
+
+end


Property changes on: DART/trunk/DART_LAB/matlab/get_ens_rank.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/get_state_increments.m
===================================================================
--- DART/trunk/DART_LAB/matlab/get_state_increments.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/get_state_increments.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,21 @@
+function [state_incs] = get_state_increments(state_ens, obs_ens, obs_incs)
+% Computes state increments given observation increments and
+% the state and obs prior ensembles
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Compute state variance and covariance
+covar = cov(state_ens, obs_ens);
+
+
+state_incs = obs_incs * covar(1, 2) / covar(2, 2);
+


Property changes on: DART/trunk/DART_LAB/matlab/get_state_increments.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/lorenz_63_adv_1step.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_63_adv_1step.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/lorenz_63_adv_1step.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,51 @@
+function[x_new, time_new] = lorenz_63_adv_1step(x, time)
+% Does a single time step advance for lorenz convective 3 variable model
+% using two step runge-kutta time step
+%
+% x is the 3-vector state, time is the 2-vector days and seconds time
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2007, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+global DELTAT
+
+% Set the fraction for the rk-2 time step
+fract = 1;
+
+% Compute the first intermediate step
+dx = comp_dt(x);
+x1 = x + fract * DELTAT * dx;
+
+% Compute the second intermediate step
+dx = comp_dt(x1);
+x2 = x1 + fract * DELTAT * dx;
+
+% New value for x is average of original value and second intermediate
+x_new = (x + x2) / 2;
+
+% Update the time ; Non-dimensional single unit for development
+time_new = time + 1;
+
+end
+
+
+function[dt] = comp_dt(x)
+
+global SIGMA
+global R
+global B
+
+dt(1) = SIGMA * (x(2) - x(1));
+dt(2) = -1 * x(1)*x(3) + R*x(1) - x(2);
+dt(3) = x(1)*x(2) - B*x(3);
+
+end
+


Property changes on: DART/trunk/DART_LAB/matlab/lorenz_63_adv_1step.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,39 @@
+% Initializes class data for L63, sets up global storage
+% and reads in control data from input file
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Lorenz-63 model parameters
+global SIGMA
+global R
+global B
+global DELTAT
+global TIME_STEP_DAYS
+global TIME_STEP_SECONDS
+
+% Set default values for the model parameters
+SIGMA = 10;
+R = 28;
+B = 8/3;
+DELTAT = 0.01;
+TIME_STEP_DAYS = 0;
+TIMES_STEP_SECONDS = 0;
+
+
+% Lorenz-63 fixed model parameters
+global MODEL_SIZE
+
+MODEL_SIZE = 3;
+
+global STATE_LOC
+
+STATE_LOC = (0:2) / 3;


Property changes on: DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,66 @@
+function[x_new, time_new] = lorenz_96_adv_1step(x, time)
+% Does a single time step advance for lorenz_96 40-variable model using four step runge-kutta time step
+%
+% x is the 40-vector state, time is the 2-vector days and seconds time
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+global DELTA_T
+
+% Compute first intermediate step
+dx = comp_dt(x);
+x1 = DELTA_T * dx;
+inter = x + x1 / 2;
+
+% Compute second intermediate step
+dx = comp_dt(inter);
+x2 = DELTA_T * dx;
+inter = x + x2 / 2;
+
+% Compute third intermediate step
+dx = comp_dt(inter);
+x3 = DELTA_T * dx;
+inter = x + x3;
+
+% Compute fourth intermediate step
+dx = comp_dt(inter);
+x4 = DELTA_T * dx;
+
+% Compute new value for x
+x_new = x + x1/6 + x2/3 + x3/3 + x4/6;
+
+% Increment time step
+time_new = time + 1;
+
+
+end
+
+%------------------------------------------------------------------------------
+
+function[dt] = comp_dt(x)
+
+global FORCING
+global MODEL_SIZE
+
+for j = 1:MODEL_SIZE
+   jp1 = j + 1;
+   if(jp1 > MODEL_SIZE) jp1 = 1; end
+   jm2 = j - 2;
+   if(jm2 < 1) jm2 = MODEL_SIZE + jm2; end
+   jm1 = j - 1;
+   if(jm1 < 1) jm1 = MODEL_SIZE; end
+
+   dt(j) = (x(jp1) - x(jm2)) * x(jm1) - x(j) + FORCING;
+end
+
+end
+


Property changes on: DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m
===================================================================
--- DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,35 @@
+% Initializes class data for L96, sets up global storage
+% and reads in control data from input file
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Lorenz-96 model parameters
+global FORCING
+global DELTA_T
+global TIME_STEP_DAYS
+global TIME_STEP_SECONDS
+
+% Set default values for the model parameters
+FORCING = 8;
+DELTA_T = 0.05;
+TIME_STEP_DAYS = 0;
+TIMES_STEP_SECONDS = 0;
+
+
+% Lorenz-96 fixed model parameters
+global MODEL_SIZE
+
+MODEL_SIZE = 40;
+
+global STATE_LOC
+
+STATE_LOC = (0:MODEL_SIZE - 1) / MODEL_SIZE;


Property changes on: DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/obs_increment_eakf.m
===================================================================
--- DART/trunk/DART_LAB/matlab/obs_increment_eakf.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/obs_increment_eakf.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,58 @@
+function [obs_increments, err] =  obs_increment_eakf(ensemble, observation, obs_error_var)
+% Computes increments for an ensemble adjustment filter
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Set error return to default successful
+err = 0;
+
+% Compute prior ensemble mean and variance
+prior_mean = mean(ensemble);
+prior_var = var(ensemble);
+
+% If both prior and observation error variance are return error
+if(prior_var <= 0 && obs_error_var <= 0)
+   err = 1;
+   return;
+end
+
+% Compute the posterior mean and variance
+% If prior variance is 0, posterior mean is prior_mean and variance is 0
+if(prior_var == 0)
+   post_mean = prior_mean;
+   post_var = 0;
+elseif(obs_error_var == 0)
+% If obs_error_var is 0, posterior mean is observation and variance is 0
+   post_mean = observation;
+   post_var = 0;
+else
+% Use product of gaussians 
+   % Compute the posterior variance
+   post_var = 1 / (1 / prior_var + 1 / obs_error_var);
+
+   % Compute posterior mean
+   post_mean = post_var * (prior_mean / prior_var + observation / obs_error_var);
+end
+
+
+% Shift the prior ensemble to have the posterior mean
+updated_ensemble = ensemble - prior_mean + post_mean;
+
+% Contract the ensemble to have the posterior_variance
+var_ratio = post_var / prior_var;
+updated_ensemble = sqrt(var_ratio) * (updated_ensemble - post_mean) + post_mean;
+
+% Compute the increments
+obs_increments = updated_ensemble - ensemble;
+
+
+


Property changes on: DART/trunk/DART_LAB/matlab/obs_increment_eakf.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/obs_increment_enkf.m
===================================================================
--- DART/trunk/DART_LAB/matlab/obs_increment_enkf.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/obs_increment_enkf.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,74 @@
+function [obs_increments, err] =  obs_increment_eakf(ensemble, observation, obs_error_var)
+% Computes increments for an ensemble Kalman filter with perturbed obs mean correction.
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Set error return to default successful
+err = 0;
+
+% Compute prior ensemble mean and variance
+prior_mean = mean(ensemble);
+prior_var = var(ensemble);
+
+% If both prior and observation error variance are zero return error
+if(prior_var <= 0 && obs_error_var <= 0)
+   err = 1;
+   return;
+end
+
+% Compute the posterior mean and variance
+% If prior variance is 0, posterior mean is prior_mean and variance is 0
+if(prior_var == 0)
+   post_mean = prior_mean;
+   post_var = 0;
+elseif(obs_error_var == 0)
+% If obs_error_var is 0, posterior mean is observation and variance is 0
+   post_mean = observation;
+   post_var = 0;
+else
+% Use product of gaussians 
+   % Compute the posterior variance
+   post_var = 1 / (1 / prior_var + 1 / obs_error_var);
+
+   % Compute posterior mean
+   post_mean = post_var * (prior_mean / prior_var + observation / obs_error_var);
+end
+
+% Generate the perturbed observations by adding 
+% draw from Normal(0, obs_error_sd)
+temp_obs = observation + sqrt(obs_error_var) * randn(size(ensemble));
+
+% Adjust so that perturbed observations have mean = to observation
+% This is technically an enhancement of earliest EnKFs
+temp_obs = temp_obs - mean(temp_obs) + observation;
+
+% Compute new ensemble members by taking product of prior ensemble
+% members and perturbed obs pairs
+updated_ens = post_var * (ensemble / prior_var + temp_obs / obs_error_var);
+
+% Increments are difference between updated and original ensemble
+obs_increments = updated_ens - ensemble;
+
+% Following are enhancements that change characteristics of EnKF
+% Coding is not finalized so do NOT just uncomment
+% Could also adjust the mean and variance of the final sample;
+% This can greatly improve and change the behavior
+% Shift the prior ensemble to have the posterior mean
+%%%updated_ensemble = ensemble - prior_mean + post_mean;
+
+% Contract the ensemble to have the posterior_variance
+%%%var_ratio = post_var / prior_var;
+%%%updated_ensemble = sqrt(var_ratio) * (updated_ensemble - post_mean) + post_mean;
+
+% Compute the increments
+%%%obs_increments = updated_ensemble - ensemble;
+


Property changes on: DART/trunk/DART_LAB/matlab/obs_increment_enkf.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/obs_increment_rhf.m
===================================================================
--- DART/trunk/DART_LAB/matlab/obs_increment_rhf.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/obs_increment_rhf.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,177 @@
+function [obs_increments, err] =  obs_increment_eakf(ensemble, observation, obs_error_var)
+% Computes increments for a rank histogram filter
+% Need to discuss the available options eventually
+% For now this implements the default options
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Set error return to default successful
+err = 0;
+
+% Get the ensemble size
+ens_size = size(ensemble, 2);
+prior_sd = std(ensemble);
+prior_var = prior_sd^2;
+
+% Sort the ensemble members and keep the indices
+[x, e_ind] = sort(ensemble);
+
+% Compute the likelihood of each member given the observation
+like = exp(-1 * (x - observation).^2 / (2 * obs_error_var));
+
+% Compute the mean likelihood density in each interior bin
+for i = 2:ens_size
+   like_dense(i) = (like(i - 1) + like(i)) / 2;
+end
+
+% For unit normal, find distance from mean to where cdf is 1/(n+1)
+dist_for_unit_sd =  -1 * weighted_norm_inv(1, 0, 1, 1/(ens_size + 1));
+
+% Variance of tails is just sample prior variance
+% Mean is adjusted so that 1/(ens_size + 1) is outside
+left_mean = x(1) + dist_for_unit_sd * prior_sd;
+left_var = prior_var;
+left_sd = prior_sd;
+
+% Same for the right tail
+right_mean = x(ens_size) - dist_for_unit_sd * prior_sd;
+right_var = prior_var;
+right_sd = prior_sd;
+
+% Eventually want to support options, for now 
+gaussian_likelihood_tails = false;
+
+if(gaussian_likelihood_tails)
+   % Need to fill this in
+else
+   % Block to do flat tails for likelihood follows
+   % This removes assumptions about likelihood and cuts cost
+   new_var_left = left_var;
+   new_sd_left = left_sd;
+   new_mean_left = left_mean;
+   prod_weight_left = like(1);
+   mass(1) = like(1) / (ens_size + 1);
+
+   % Same for right tail
+   new_var_right = right_var;
+   new_sd_right = right_sd;
+   new_mean_right = right_mean;
+   prod_weight_right = like(ens_size);
+   mass(ens_size + 1) = like(ens_size) / (ens_size + 1);
+   % End block for flat tail likelihood
+end
+
+% The mass in each interior box is the height times the width
+% The height of the likelihood is like_dense
+% For the prior, mass is 1 / ((n+1) width), and mass = height x widt so
+% The height of the prior is 1 / ((n+1) width); Multiplying by width leaves 1/(n+1)
+
+% In prior, have 1/(n+1) mass in each bin, multiply by mean likelihood
+% density to get approximate mass in updated bin
+for i = 2:ens_size
+   mass(i) = like_dense(i) / (ens_size + 1);
+   % Height of prior in this bin is mass/width; Only needed for trapezoidal
+   % If two ensemble members are the same, set height to -1 as flag
+   if(x(i) == x(i-1)) 
+      height(i) = -1;
+   else
+      height(i) = 1 / ((ens_size + 1) * (x(i) - x(i-1)));
+   end
+end
+
+% Now normalize the mass in the different bins to get a pdf
+mass_sum = sum(mass);
+nmass = mass / mass_sum;
+
+% 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;
+% This is the same as right_amp=(ens_size + 1)*nmass(ens_size+1)
+right_amp = prod_weight_right / mass_sum;
+
+% Find cumulative mass at each box boundary and middle boundary
+cumul_mass(1) = 0;
+for i = 1:ens_size + 1
+   cumul_mass(i+1) = cumul_mass(i) + nmass(i);
+end
+
+% Begin internal box search at bottom of lowest box
+lowest_box = 1;
+
+% Find each new ensemble member's location
+for i = 1:ens_size
+   % Each update ensemble member has 1/(n+1) mass before it
+   umass = i / (ens_size + 1);
+
+   % 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);
+   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
+      new_ens(i) = weighted_norm_inv(right_amp, new_mean_right, ...
+         new_sd_right, 1 - umass);
+      % Coming in from the right, use symmetry after pretending it's on left
+      new_ens(i) = new_mean_right + (new_mean_right - new_ens(i));
+   else
+      % In one of the inner boxes
+      for j = lowest_box:ens_size - 1
+         % Find the box that this mass is in
+         if(umass >= cumul_mass(j+1) && umass <= cumul_mass(j+2))
+
+            % Only rectangular is implemented for now
+            rectangular_quadrature = true;
+
+            if(rectangular_quadrature)
+               % Rectangular quadrature block first
+               % Linearly interpolate in mass
+               new_ens(i) = x(j) + ((umass - cumul_mass(j+1)) / ...
+                  (cumul_mass(j+2) - cumul_mass(j+1))) * (x(j+1) - x(j));
+            else
+               % Trapezoidal interpolation block goes here
+            end
+
+            % Don't need to search lower boxes again
+            lowest_box = j;
+            break
+         end
+      end
+   end
+end
+
+% Convert to increments for unsorted
+for i = 1:ens_size
+   obs_increments(e_ind(i)) = new_ens(i) - x(i);
+end
+
+
+%-----------------------------------------------
+
+
+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 = norminv(np, 0, 1);
+
+% Add in the mean and normalize by sd
+x = mean + x * sd;
+


Property changes on: DART/trunk/DART_LAB/matlab/obs_increment_rhf.m
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + text/x-matlab
Name: svn:keywords
   + Date Rev Author URL Id
Name: svn:eol-style
   + native

Added: DART/trunk/DART_LAB/matlab/oned_ensemble.fig
===================================================================
(Binary files differ)


Property changes on: DART/trunk/DART_LAB/matlab/oned_ensemble.fig
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + application/octet-stream

Added: DART/trunk/DART_LAB/matlab/oned_ensemble.m
===================================================================
--- DART/trunk/DART_LAB/matlab/oned_ensemble.m	                        (rev 0)
+++ DART/trunk/DART_LAB/matlab/oned_ensemble.m	2009-08-25 23:04:35 UTC (rev 4011)
@@ -0,0 +1,473 @@
+function varargout = oned_ensemble(varargin)
+% ONED_ENSEMBLE explore the details of ensemble data assimilation for a scalar.
+%
+%      Click on the 'Create New Ensemble' button to activate the interactive 
+%      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.)
+%
+%      After you have an ensemble and an observation, choose an assimilation 
+%      algorithm and click 'Update Ensemble'. The algorithm is applied and the 
+%      Posterior (blue) is plotted below the Prior (green).  Choose 'EAKF' 
+%      and click 'Update' ...  multiple times. Do the same for 'EnKF' and
+%      'RHF'. Hmnnnn ....
+%
+%      change the Observation Error SD, lay down an ensemble pretty far away
+%      from the observation - have fun with it.
+%
+% See also: gaussian_product, oned_model, twod_ensemble, run_lorenz_63, 
+%           run_lorenz_96
+
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Last Modified by GUIDE v2.5 25-Mar-2009 14:34:04
+
+% Begin initialization code - DO NOT EDIT
+gui_Singleton = 1;
+gui_State = struct('gui_Name',       mfilename, ...
+                   'gui_Singleton',  gui_Singleton, ...
+                   'gui_OpeningFcn', @oned_ensemble_OpeningFcn, ...
+                   'gui_OutputFcn',  @oned_ensemble_OutputFcn, ...
+                   'gui_LayoutFcn',  [] , ...
+                   'gui_Callback',   []);
+if nargin && ischar(varargin{1})
+    gui_State.gui_Callback = str2func(varargin{1});
+end
+
+if nargout
+    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
+else
+    gui_mainfcn(gui_State, varargin{:});
+end
+% End initialization code - DO NOT EDIT
+
+
+% --- Executes just before oned_ensemble is made visible.
+function oned_ensemble_OpeningFcn(hObject, eventdata, handles, varargin)
+% This function has no output args, see OutputFcn.
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+% varargin   command line arguments to oned_ensemble (see VARARGIN)
+
+help oned_ensemble
+
+% Choose default command line output for oned_ensemble
+handles.output = hObject;
+
+% Insert the ensemble structure into this
+handles.ens_size = 0;
+handles.ens_members = 0;
+handles.h_obs_plot = 0;
+handles.h_update_ens = 0;
+handles.h_ens_member = 0;
+handles.h_obs_ast = 0;
+
+% Update handles structure
+guidata(hObject, handles);
+
+% Go ahead and plot the initial observational error distribution
+h_observation = get(handles.edit1);
+h_obs_error_sd = get(handles.edit2);
+observation = str2double(h_observation.String);
+obs_error_sd = str2double(h_obs_error_sd.String);
+handles.h_obs_plot = plot_gaussian(observation, obs_error_sd, 1);
+set(handles.h_obs_plot, 'Color', 'r', 'Linestyle', '--', 'Linewidth', 2);
+hold on
+
+% Plot an asterisk 
+handles.h_obs_ast = plot(observation, 0, 'r*', 'MarkerSize', 16);
+
+% Set a basic plotting domain range that includes mean +/- 3 obs SDs
+lower = observation - 3*obs_error_sd;
+upper = observation + 3*obs_error_sd;
+axis([lower upper -0.2 1]);
+
+
+set(gca, 'YTick', [0 0.2 0.4 0.6 0.8]);
+set(gca, 'YTickLabel', [0 0.2 0.4 0.6 0.8]);
+
+hold on
+plot([lower upper], [0 0], 'k', 'Linewidth', 2);
+
+% Update handles structure
+guidata(hObject, handles);
+
+
+% UIWAIT makes oned_ensemble wait for user response (see UIRESUME)
+% uiwait(handles.figure1);
+
+
+%----------------------------------------------------------------------
+
+
+% --- Outputs from this function are returned to the command line.
+function varargout = oned_ensemble_OutputFcn(hObject, eventdata, handles) 
+% varargout  cell array for returning output args (see VARARGOUT);
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Get default command line output from handles structure
+varargout{1} = handles.output;
+
+
+%----------------------------------------------------------------------
+
+
+% --- Executes on button press in pushbutton1.

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list