[Dart-dev] [8125] DART/trunk/DART_LAB/matlab:
nancy at ucar.edu
nancy at ucar.edu
Wed Jun 24 15:20:35 MDT 2015
Revision: 8125
Author: thoar
Date: 2015-06-24 15:20:35 -0600 (Wed, 24 Jun 2015)
Log Message:
-----------
Moving all the 'helper' functions to be private.
This will make the functions the user is supposed to use
more visible and the directory less daunting.
Added Paths:
-----------
DART/trunk/DART_LAB/matlab/private/
DART/trunk/DART_LAB/matlab/private/advance_oned.m
DART/trunk/DART_LAB/matlab/private/comp_cov_factor.m
DART/trunk/DART_LAB/matlab/private/g_prod_plot.m
DART/trunk/DART_LAB/matlab/private/get_ens_rank.m
DART/trunk/DART_LAB/matlab/private/get_state_increments.m
DART/trunk/DART_LAB/matlab/private/kurt.m
DART/trunk/DART_LAB/matlab/private/lorenz_63_adv_1step.m
DART/trunk/DART_LAB/matlab/private/lorenz_63_static_init_model.m
DART/trunk/DART_LAB/matlab/private/lorenz_96_adv_1step.m
DART/trunk/DART_LAB/matlab/private/lorenz_96_static_init_model.m
DART/trunk/DART_LAB/matlab/private/norm_inv.m
DART/trunk/DART_LAB/matlab/private/norm_pdf.m
DART/trunk/DART_LAB/matlab/private/obs_increment_eakf.m
DART/trunk/DART_LAB/matlab/private/obs_increment_enkf.m
DART/trunk/DART_LAB/matlab/private/obs_increment_rhf.m
DART/trunk/DART_LAB/matlab/private/plot_gaussian.m
DART/trunk/DART_LAB/matlab/private/plot_polar.m
DART/trunk/DART_LAB/matlab/private/product_of_gaussians.m
-------------- next part --------------
Copied: DART/trunk/DART_LAB/matlab/private/advance_oned.m (from rev 8120, DART/trunk/DART_LAB/matlab/advance_oned.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/advance_oned.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/advance_oned.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,28 @@
+function x_new = advance_oned(x, alpha, model_bias)
+%% advance_oned(x, alpha, model_bias)
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+x_new = x + comp_dt(x, alpha, model_bias);
+end
+
+%---------------------------------------------------
+
+% Internal function comp_dt
+function dx = comp_dt(x, alpha, model_bias)
+
+% Compute the time tendency; alpha controls nonlinearity
+% model_bias controls a shift in the model dynamics
+dx = (x + model_bias) + alpha .* x .* abs(x);
+
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/comp_cov_factor.m (from rev 8120, DART/trunk/DART_LAB/matlab/comp_cov_factor.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/comp_cov_factor.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/comp_cov_factor.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,26 @@
+function cov_factor = comp_cov_factor(z_in, c)
+%% comp_cov_factor Gaspari Cohn cutoff, z_in is the distance while c is the cutoff
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+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
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/g_prod_plot.m (from rev 8120, DART/trunk/DART_LAB/matlab/g_prod_plot.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/g_prod_plot.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/g_prod_plot.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,108 @@
+function [prior_mean, prior_sd, obs_mean, obs_err_sd, is_err] = g_prod_plot(h)
+%% g_prod_plot Updates the plot of the prior and observation gaussians
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% Successful return as default
+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
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/get_ens_rank.m (from rev 8120, DART/trunk/DART_LAB/matlab/get_ens_rank.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/get_ens_rank.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/get_ens_rank.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,22 @@
+function [rank] = get_ens_rank(ens, x)
+%% get_ens_rank
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+s_ens = sort(ens);
+rank = find(s_ens < squeeze(x),1,'last') + 1;
+if(isempty(rank))
+ rank = 1;
+end
+
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/get_state_increments.m (from rev 8120, DART/trunk/DART_LAB/matlab/get_state_increments.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/get_state_increments.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/get_state_increments.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,21 @@
+function [state_incs] = get_state_increments(state_ens, obs_ens, obs_incs)
+%% get_state_increments Computes state increments given observation increments and
+% the state and obs prior ensembles
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% Compute state variance and covariance
+covar = cov(state_ens, obs_ens);
+
+
+state_incs = obs_incs * covar(1, 2) / covar(2, 2);
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/kurt.m (from rev 8120, DART/trunk/DART_LAB/matlab/kurt.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/kurt.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/kurt.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,30 @@
+function k = kurt(vals)
+% computes the kurtosis of the given input array
+%
+% based on the second formula on this web page:
+% http://www.ats.ucla.edu/stat/mult_pkg/faq/general/kurtosis.htm
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% array of diffs from mean
+del = vals - mean(vals);
+
+% compute the square and 4th power of the diffs from mean
+m2 = mean(del .^ 2);
+m4 = mean(del .^ 4);
+
+% compute the kurtosis value. this is the version
+% of the kurtosis formula that is not nonbiased and
+% does not subtract 3.0 from the result.
+k = (m4 ./ (m2 .^ 2));
+
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
Copied: DART/trunk/DART_LAB/matlab/private/lorenz_63_adv_1step.m (from rev 8120, DART/trunk/DART_LAB/matlab/lorenz_63_adv_1step.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/lorenz_63_adv_1step.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/lorenz_63_adv_1step.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,51 @@
+function[x_new, time_new] = lorenz_63_adv_1step(x, time)
+%% lorenz_63)adv_1step advances the lorenz convective 3 variable model
+% for a single two step runge-kutta time step
+%
+% x is the 3-vector state, time is the 2-vector days and seconds time
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+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
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/lorenz_63_static_init_model.m (from rev 8120, DART/trunk/DART_LAB/matlab/lorenz_63_static_init_model.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/lorenz_63_static_init_model.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/lorenz_63_static_init_model.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,40 @@
+%% lorenz_63_static_init_model Initializes class data for L63, sets up global storage
+% and reads in control data from input file
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% 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;
+TIME_STEP_SECONDS = 0;
+
+
+% Lorenz-63 fixed model parameters
+global MODEL_SIZE
+
+MODEL_SIZE = 3;
+
+global STATE_LOC
+
+STATE_LOC = (0:2) / 3;
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/lorenz_96_adv_1step.m (from rev 8120, DART/trunk/DART_LAB/matlab/lorenz_96_adv_1step.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/lorenz_96_adv_1step.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/lorenz_96_adv_1step.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,70 @@
+function[x_new, time_new] = lorenz_96_adv_1step(x, time)
+%% lorenz_96_adv_1step 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
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+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
+
+dt = zeros(1,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
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/lorenz_96_static_init_model.m (from rev 8120, DART/trunk/DART_LAB/matlab/lorenz_96_static_init_model.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/lorenz_96_static_init_model.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/lorenz_96_static_init_model.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,36 @@
+%% lorenz_96_static_init_model Initializes class data for L96, sets up global storage
+% and reads in control data from input file
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% 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;
+TIME_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;
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/norm_inv.m (from rev 8120, DART/trunk/DART_LAB/matlab/norm_inv.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/norm_inv.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/norm_inv.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,83 @@
+function x = norm_inv(p)
+% Computes inverse CDF for the normal distribution,
+% evaluated at either the scalar value or array P.
+%
+% based on the code in DART assim_tools_mod.f90, which is
+% in turn based on http://home.online.no/~pjacklam/notes/invnorm
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% constants where precision is vital
+a1 = -39.69683028665376;
+a2 = 220.9460984245205;
+a3 = -275.9285104469687;
+a4 = 138.357751867269;
+a5 = -30.66479806614716;
+a6 = 2.506628277459239;
+b1 = -54.4760987982241;
+b2 = 161.5858368580409;
+b3 = -155.6989798598866;
+b4 = 66.80131188771972;
+b5 = -13.28068155288572;
+c1 = -0.007784894002430293;
+c2 = -0.3223964580411365;
+c3 = -2.400758277161838;
+c4 = -2.549732539343734;
+c5 = 4.374664141464968;
+c6 = 2.938163982698783;
+d1 = 0.007784695709041462;
+d2 = 0.3224671290700398;
+d3 = 2.445134137142996;
+d4 = 3.754408661907416;
+
+% Split into an inner and two outer regions which have separate fits
+p_low = 0.02425;
+p_high = 1 - p_low;
+
+% if P is scalar, simplify
+if (numel(p) == 1)
+ if (p < p_low)
+ q = sqrt(-2.0 * log(p));
+ x = (((((c1*q + c2)*q + c3)*q + c4)*q + c5)*q + c6) / ...
+ ((((d1*q + d2)*q + d3)*q + d4)*q + 1.0);
+ elseif (p >= p_low) && (p <= p_high)
+ q = p - 0.5;
+ r = q*q;
+ x = (((((a1*r + a2)*r + a3)*r + a4)*r + a5)*r + a6)*q / ...
+ (((((b1*r + b2)*r + b3)*r + b4)*r + b5)*r + 1.0);
+
+ else % (p > p_high)
+ q = sqrt(-2.0 * log(1.0 - p));
+ x = -(((((c1*q + c2)*q + c3)*q + c4)*q + c5)*q + c6) / ...
+ ((((d1*q + d2)*q + d3)*q + d4)*q + 1.0);
+ end
+
+else % if P is an array, so do it element by element
+
+ inds = (p < p_low);
+ q(inds) = sqrt(-2.0 * log(p(inds)));
+ x(inds) = (((((c1*q(inds) + c2).*q(inds) + c3).*q(inds) + c4).*q(inds) + c5).*q(inds) + c6) ./ ...
+ ((((d1*q(inds) + d2).*q(inds) + d3).*q(inds) + d4).*q(inds) + 1.0);
+
+ inds = (p >= p_low) & (p <= p_high);
+ q(inds) = p(inds) - 0.5;
+ r(inds) = q(inds).*q(inds);
+ x(inds) = (((((a1*r(inds) + a2).*r(inds) + a3).*r(inds) + a4).*r(inds) + a5).*r(inds) + a6).*q(inds) ./ ...
+ (((((b1*r(inds) + b2).*r(inds) + b3).*r(inds) + b4).*r(inds) + b5).*r(inds) + 1.0);
+
+ inds = (p > p_high);
+ q(inds) = sqrt(-2.0 * log(1.0 - p(inds)));
+ x(inds) = -(((((c1*q(inds) + c2).*q(inds) + c3).*q(inds) + c4).*q(inds) + c5).*q(inds) + c6) ./ ...
+ ((((d1*q(inds) + d2).*q(inds) + d3).*q(inds) + d4).*q(inds) + 1.0);
+end
+
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
Copied: DART/trunk/DART_LAB/matlab/private/norm_pdf.m (from rev 8120, DART/trunk/DART_LAB/matlab/norm_pdf.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/norm_pdf.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/norm_pdf.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,32 @@
+function [y] = norm_pdf(x, mu, sigma)
+% computes a gaussian (normal) PDF
+% for the points of X with a given mean (mu) and standard deviation (sigma)
+%
+% normal plot, y given x:
+% y = (1 / (sigma * sqrt(2*pi))) * e ^ ((-1/2 * ((x-mu) / sigma)^2)
+% or
+% g(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{ -\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2 }.
+%
+% see: https://en.wikipedia.org/wiki/Probability_density_function
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+
+e = exp(1);
+
+basen = (1.0 / (sigma * sqrt(2*pi)));
+expon = -0.5 * (((x-mu) / sigma).^2 );
+
+y = basen * (e .^ expon);
+
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/obs_increment_eakf.m (from rev 8120, DART/trunk/DART_LAB/matlab/obs_increment_eakf.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/obs_increment_eakf.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/obs_increment_eakf.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,56 @@
+function [obs_increments, err] = obs_increment_eakf(ensemble, observation, obs_error_var)
+%% obs_increment_eakf Computes increments for an ensemble adjustment filter
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% 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;
+
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/obs_increment_enkf.m (from rev 8120, DART/trunk/DART_LAB/matlab/obs_increment_enkf.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/obs_increment_enkf.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/obs_increment_enkf.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,74 @@
+function [obs_increments, err] = obs_increment_enkf(ensemble, observation, obs_error_var)
+%% obs_increment_enkf Computes increments for an ensemble Kalman filter with perturbed obs mean correction.
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% 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;
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/obs_increment_rhf.m (from rev 8120, DART/trunk/DART_LAB/matlab/obs_increment_rhf.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/obs_increment_rhf.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/obs_increment_rhf.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,185 @@
+function [obs_increments, err] = obs_increment_rhf(ensemble, observation, obs_error_var)
+%% obs_increment_rhf Computes increments for a rank histogram filter
+% Need to discuss the available options eventually
+% For now this implements the default options
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% 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;
+
+% allocate space for some of the density calculation bits
+like_dense = zeros(1,ens_size); % likelihood density
+mass = zeros(1,ens_size);
+height = zeros(1,ens_size);
+
+% Sort the ensemble members and keep the indices
+[x, e_ind] = sort(ensemble);
+
+% 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 width so
+% The height of the prior is 1 / ((n+1) width); Multiplying by width leaves 1/(n+1)
+
+% In prior, have 1/(n+1) mass in each bin, multiply by mean likelihood
+% density to get approximate mass in updated bin
+
+for i = 2:ens_size
+ mass(i) = like_dense(i) / (ens_size + 1);
+ % Height of prior in this bin is mass/width; Only needed for trapezoidal
+ % If two ensemble members are the same, set height to -1 as flag
+ if(x(i) == x(i-1))
+ 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 = zeros(1,ens_size+1);
+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;
+
+new_ens = zeros(1,ens_size);
+obs_increments = zeros(1,ens_size);
+% Find each new ensemble member's location
+for i = 1:ens_size
+ % Each update ensemble member has 1/(n+1) mass before it
+ 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 = norm_inv(np);
+
+% Add in the mean and normalize by sd
+x = mean + x * sd;
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/plot_gaussian.m (from rev 8120, DART/trunk/DART_LAB/matlab/plot_gaussian.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/plot_gaussian.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/plot_gaussian.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,26 @@
+function[plot_handle] = plot_gaussian(mymean, sd, weight)
+%% plot_gaussian Plot gaussian over 5 standard deviations
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+x_min = mymean - 5*sd;
+x_max = mymean + 5*sd;
+x_range = x_max - x_min;
+
+% Number of points is 1001
+num_points = 1001;
+interval = x_range / num_points;
+x = x_min:interval:x_max;
+y = weight * norm_pdf(x, mymean, sd);
+
+plot_handle = plot(x, y);
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/plot_polar.m (from rev 8120, DART/trunk/DART_LAB/matlab/plot_polar.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/plot_polar.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/plot_polar.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,21 @@
+function h = plot_polar(y, x, mean_dist, string, model_size)
+%% plot_polar
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% Y includes a wraparound point, x does not
+x_t(model_size + 1) = x(1);
+x_t(1:model_size) = x;
+h = polar(y, mean_dist + x_t, string);
+
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Copied: DART/trunk/DART_LAB/matlab/private/product_of_gaussians.m (from rev 8120, DART/trunk/DART_LAB/matlab/product_of_gaussians.m)
===================================================================
--- DART/trunk/DART_LAB/matlab/private/product_of_gaussians.m (rev 0)
+++ DART/trunk/DART_LAB/matlab/private/product_of_gaussians.m 2015-06-24 21:20:35 UTC (rev 8125)
@@ -0,0 +1,30 @@
+function [post_mean, post_sd, weight] = ...
+ product_of_gaussians(prior_mean, prior_sd, obs, obs_err_sd)
+%% product_of_gaussians Computes mean, variance and weight of the product of two unit gaussians given the mean and standard deviation of each.
+
+%% DART software - Copyright 2004 - 2013 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% DART $Id$
+
+% Get the prior and observational error variance
+prior_var = prior_sd^2;
+obs_err_var = obs_err_sd^2;
+
+% Compute the posterior variance
+post_var = 1. / (1. / prior_var + 1. / obs_err_var);
+post_sd = sqrt(post_var);
+
+% Compute the posterior mean
+post_mean = post_var * (prior_mean / prior_var + obs / obs_err_var);
+
+% Compute the associated weight
+weight = (1. / (sqrt(2. * pi) * sqrt(prior_var + obs_err_var))) *...
+ exp(-0.5 * (obs - prior_mean).^2 ./ (prior_var + obs_err_var));
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
More information about the Dart-dev
mailing list