merging in the changes from the trunk (to the DART_LAB functions -

adding Moha's inflation-enhanced GUIs)



git-svn-id: https://svn-dares-dart.cgd.ucar.edu/DART/releases/Lanai@11679 dfa8782c-da17-4c45-ba5c-5625b50a00d6

Former-commit-id: ec11f1a997992e13f193896d42b95878416d8e6d

DART_LAB/matlab/private/compute_new_density.m

@@ -0,0 +1,26 @@
+function density = compute_new_density(dist_2, sigma_p_2, sigma_o_2, lambda_mean, lambda_sd, gamma, lambda)
+
+%% 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$
+
+
+% Compute probability of this lambda being correct
+exponent_prior = - 0.5 * (lambda - lambda_mean)^2 / lambda_sd^2;
+
+% Compute probability that observation would have been observed given this lambda
+theta_2 = (1.0 + gamma * (sqrt(lambda) - 1.0))^2 * sigma_p_2 + sigma_o_2;
+theta = sqrt(theta_2);
+
+exponent_likelihood = dist_2 / ( -2.0 * theta_2);
+
+% Compute the updated probability density for lambda
+% Have 1 / sqrt(2 PI) twice, so product is 1 / (2 PI)
+density = exp(exponent_likelihood + exponent_prior) / (2.0 * pi * lambda_sd * theta);
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$

DART_LAB/matlab/private/plot_polar.m

@@ -10,7 +10,7 @@
 % 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);
+h = polar_dares(y, mean_dist + x_t, string);
 
 end
 

DART_LAB/matlab/private/polar_dares.m

@@ -0,0 +1,230 @@
+function hpol = polar_dares(THETA, RHO, SPEC)
+% This is a modified version of MATLAB's polar function. 
+% The modifications are done in such a way that the visualization of the
+% DART tutorial is faster. 
+%       - The X- and Y- labels are turned off. 
+%       - The function allows for zooming 
+%
+% DART $Id$
+
+    if nargin > 2
+        varargin{1} = THETA;
+        varargin{2} = RHO;
+        varargin{3} = SPEC;
+    else
+        varargin{1} = THETA;
+        varargin{2} = RHO;
+    end
+
+    % Parse possible Axes input
+    [cax, args, nargs] = axescheck(varargin{:});
+
+    if nargs < 1
+        error(message('MATLAB:narginchk:notEnoughInputs'));
+    elseif nargs > 3
+        error(message('MATLAB:narginchk:tooManyInputs'));
+    end
+
+    if nargs < 1 || nargs > 3
+        error(message('MATLAB:polar:InvalidDataInputs'));
+    elseif nargs == 2
+        theta = args{1};
+        rho = args{2};
+        if ischar(rho)
+            line_style = rho;
+            rho = theta;
+            [mr, nr] = size(rho);
+            if mr == 1
+                theta = 1 : nr;
+            else
+                th = (1 : mr)';
+                theta = th(:, ones(1, nr));
+            end
+        else
+            line_style = 'auto';
+        end
+    elseif nargs == 1
+        theta = args{1};
+        line_style = 'auto';
+        rho = theta;
+        [mr, nr] = size(rho);
+        if mr == 1
+            theta = 1 : nr;
+        else
+            th = (1 : mr)';
+            theta = th(:, ones(1, nr));
+        end
+    else % nargs == 3
+        [theta, rho, line_style] = deal(args{1 : 3});
+    end
+    if ischar(theta) || ischar(rho)
+        error(message('MATLAB:polar:InvalidInputType'));
+    end
+    if ~isequal(size(theta), size(rho))
+        error(message('MATLAB:polar:InvalidInputDimensions'));
+    end
+    try
+        theta = full(double(theta));
+        rho = full(double(rho));
+    catch
+        error(message('MATLAB:specgraph:private:specgraph:nonNumericInput'));
+    end
+
+    % get hold state
+    cax = newplot(cax);
+
+    next = lower(get(cax, 'NextPlot'));
+    hold_state = ishold(cax);
+
+    if isa(handle(cax),'matlab.graphics.axis.PolarAxes')
+        error(message('MATLAB:polar:PolarAxes'));
+    end
+
+    % get x-axis text color so grid is in same color
+    % get the axis gridColor
+    axColor = get(cax, 'Color');
+    gridAlpha = get(cax, 'GridAlpha');
+    axGridColor = get(cax,'GridColor').*gridAlpha + axColor.*(1-gridAlpha);
+    tc = axGridColor;
+    ls = get(cax, 'GridLineStyle');
+
+    % Hold on to current Text defaults, reset them to the
+    % Axes' font attributes so tick marks use them.
+    fAngle = get(cax, 'DefaultTextFontAngle');
+    fName = get(cax, 'DefaultTextFontName');
+    fSize = get(cax, 'DefaultTextFontSize');
+    fWeight = get(cax, 'DefaultTextFontWeight');
+    fUnits = get(cax, 'DefaultTextUnits');
+    set(cax, ...
+        'DefaultTextFontAngle', get(cax, 'FontAngle'), ...
+        'DefaultTextFontName', get(cax, 'FontName'), ...
+        'DefaultTextFontSize', get(cax, 'FontSize'), ...
+        'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
+        'DefaultTextUnits', 'data');
+
+    % only do grids if hold is off
+    if ~hold_state
+
+        % make a radial grid
+        hold(cax, 'on');
+        % ensure that Inf values don't enter into the limit calculation.
+        arho = abs(rho(:));
+        maxrho = max(arho(arho ~= Inf));
+        hhh = line([-maxrho, -maxrho, maxrho, maxrho], [-maxrho, maxrho, maxrho, -maxrho], 'Parent', cax);
+        set(cax, 'DataAspectRatio', [1, 1, 1], 'PlotBoxAspectRatioMode', 'auto');
+        v = [get(cax, 'XLim') get(cax, 'YLim')];
+        ticks = sum(get(cax, 'YTick') >= 0);
+        delete(hhh);
+        % check radial limits and ticks
+        rmin = 0;
+        rmax = v(4);
+        rticks = max(ticks - 1, 2);
+        if rticks > 5   % see if we can reduce the number
+            if rem(rticks, 2) == 0
+                rticks = rticks / 2;
+            elseif rem(rticks, 3) == 0
+                rticks = rticks / 3;
+            end
+        end
+
+        % define a circle
+        th = 0 : pi / 50 : 2 * pi;
+        xunit = cos(th);
+        yunit = sin(th);
+        % now really force points on x/y axes to lie on them exactly
+        inds = 1 : (length(th) - 1) / 4 : length(th);
+        xunit(inds(2 : 2 : 4)) = zeros(2, 1);
+        yunit(inds(1 : 2 : 5)) = zeros(3, 1);
+        % plot background if necessary
+        if ~ischar(get(cax, 'Color'))
+            patch('XData', xunit * rmax, 'YData', yunit * rmax, ...
+                'EdgeColor', tc, 'FaceColor', get(cax, 'Color'), ...
+                'HandleVisibility', 'off', 'Parent', cax);
+        end
+
+        % draw radial circles
+        c82 = cos(82 * pi / 180);
+        s82 = sin(82 * pi / 180);
+        rinc = (rmax - rmin) / rticks;
+        for i = (rmin + rinc) : rinc : rmax - 1
+            hhh = line(xunit * i, yunit * i, 'LineStyle', ls, 'Color', tc, 'LineWidth', 1, ...
+                'HandleVisibility', 'off', 'Parent', cax);
+            text((i + rinc / 20) * c82, (i + rinc / 20) * s82, ...
+                ['  ' num2str(i)], 'VerticalAlignment', 'bottom', ...
+                'HandleVisibility', 'off', 'Parent', cax);
+        end
+        set(hhh, 'LineStyle', '-'); % Make outer circle solid
+
+        % plot spokes
+        th = (1 : 6) * 2 * pi / 12;
+        cst = cos(th);
+        snt = sin(th);
+        cs = [-cst; cst];
+        sn = [-snt; snt];
+        line(rmax * cs, rmax * sn, 'LineStyle', ls, 'Color', tc, 'LineWidth', 1, ...
+            'HandleVisibility', 'off', 'Parent', cax);
+
+        % annotate spokes in degrees
+        rt = 1.1 * rmax;
+        for i = 1 : length(th)
+            text(rt * cst(i), rt * snt(i), int2str(i * 30),...
+                'HorizontalAlignment', 'center', ...
+                'HandleVisibility', 'off', 'Parent', cax, 'FontSize', 12);
+            if i == length(th)
+                loc = int2str(0);
+            else
+                loc = int2str(180 + i * 30);
+            end
+            text(-rt * cst(i), -rt * snt(i), loc, 'HorizontalAlignment', 'center', ...
+                'HandleVisibility', 'off', 'Parent', cax, 'FontSize', 12);
+        end
+
+        % set view to 2-D
+        view(cax, 2);
+        % set axis limits
+        axis(cax, rmax * [-1, 1, -1.15, 1.15]);
+    end
+
+    % Reset defaults.
+    set(cax, ...
+        'DefaultTextFontAngle', fAngle , ...
+        'DefaultTextFontName', fName , ...
+        'DefaultTextFontSize', fSize, ...
+        'DefaultTextFontWeight', fWeight, ...
+        'DefaultTextUnits', fUnits );
+
+    % transform data to Cartesian coordinates.
+    xx = rho .* cos(theta);
+    yy = rho .* sin(theta);
+
+    % plot data on top of grid
+    if strcmp(line_style, 'auto')
+        q = plot(xx, yy, 'Parent', cax);
+    else
+        q = plot(xx, yy, line_style, 'Parent', cax);
+    end
+
+    if nargout == 1
+        hpol = q;
+    end
+
+    if ~hold_state
+        set(cax, 'DataAspectRatio', [1, 1, 1]), axis(cax, 'off');
+        set(cax, 'NextPlot', next);
+    end
+
+    % Disable pan and zoom
+    p = hggetbehavior(cax, 'Pan');
+    p.Enable = true;
+    z = hggetbehavior(cax, 'Zoom');
+    z.Enable = true;
+
+    if ~isempty(q) && ~isdeployed
+        makemcode('RegisterHandle', cax, 'IgnoreHandle', q, 'FunctionName', 'polar');
+    end
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$

DART_LAB/matlab/private/update_inflate.m

@@ -0,0 +1,103 @@
+function [new_cov_inflate, new_cov_inflate_sd] = update_inflate(x, sigma_p_2, obs, sigma_o_2, inflate_prior_val, lambda_mean, ...
+                                                    lambda_mean_LB, lambda_mean_UB, gamma, lambda_sd, lambda_sd_LB)
+% Adaptive scheme to update the inflation both in space and time. 
+% Based on the algorithm in the following article:
+% Anderson, J. L., 2009: Spatially and temporally varying adaptive covariance 
+% inflation for ensemble filters. Tellus A, 61, 72-83. doi: 10.1111/j.1600-0870.2008.00361.x 
+%
+% More documentation available at:
+% https://svn-dares-dart.cgd.ucar.edu/DART/releases/Manhattan/assimilation_code/modules/assimilation/adaptive_inflate_mod.html
+
+%% 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$
+
+
+%% FIRST, update the inflation mean:
+
+% Get the "non-inflated" variance of the sample
+% lambda here, is the prior value before the update.
+sigma_p_2 = sigma_p_2 / ( 1 + gamma*(sqrt(inflate_prior_val) - 1) )^2;
+
+% Squared-innovation
+dist_2 = (x - obs)^2;
+
+% d (distance between obs and ens) is drawn from Gaussian with 0 mean and
+% variance: E(d^2) = \lambda^o*\sigma_p^2 + \sigma_o^2
+theta_bar_2  = ( 1 + gamma * (sqrt(lambda_mean) - 1) )^2 * sigma_p_2 + sigma_o_2;
+theta_bar    = sqrt(theta_bar_2);
+u_bar        = 1 / (sqrt(2 * pi) * theta_bar);
+like_exp_bar = - 0.5 * dist_2 / theta_bar_2;
+v_bar        = exp(like_exp_bar);
+
+gamma_terms  = 1 - gamma + gamma*sqrt(lambda_mean);
+dtheta_dinf  = 0.5 * sigma_p_2 * gamma * gamma_terms / (theta_bar * sqrt(lambda_mean));              
+
+% The likelihood: p(d/lambda)
+like_bar     = u_bar * v_bar;
+
+% Derivative of the likelihood; evaluated at the current inflation mean
+like_prime   = (like_bar * dtheta_dinf / theta_bar) * (dist_2 / theta_bar_2 - 1);
+like_ratio   = like_bar / like_prime;
+
+% Solve a quadratic equation
+a = 1;
+b = like_ratio - 2*lambda_mean;
+c = lambda_mean^2 - lambda_sd^2 - like_ratio*lambda_mean ;
+
+o = max( [ abs(a), abs(b), abs(c) ] );
+a = a/o;
+b = b/o;
+c = c/o;
+d = b^2 - 4*a*c;
+
+if b < 0
+    s1 = 0.5 * ( -b + sqrt(d) ) / a;
+else 
+    s1 = 0.5 * ( -b - sqrt(d) ) / a;
+end
+s2 = ( c/a ) / s1;
+
+% Select the updated-mean that is closest to the prior
+if abs(s2 - lambda_mean) < abs(s1 - lambda_mean)
+   new_cov_inflate = s2;
+else
+   new_cov_inflate = s1;
+end
+
+% Make sure the update is not smaller than the lower bound
+if new_cov_inflate < lambda_mean_LB || new_cov_inflate > lambda_mean_UB || isnan(new_cov_inflate)
+    new_cov_inflate = lambda_mean_LB; 
+    new_cov_inflate_sd = lambda_sd;
+    return
+end
+
+
+%% Now, update the inflation variance
+if lambda_sd <= lambda_sd_LB 
+    new_cov_inflate_sd = lambda_sd;
+    return
+else
+    % First compute the new_max value for normalization purposes
+    new_max = compute_new_density(dist_2, sigma_p_2, sigma_o_2, lambda_mean, lambda_sd, gamma, new_cov_inflate);
+
+    % Find value at a point one OLD sd above new mean value
+    new_1_sd = compute_new_density(dist_2, sigma_p_2, sigma_o_2, lambda_mean, lambda_sd, gamma, new_cov_inflate + lambda_sd);
+
+    ratio = new_1_sd / new_max;
+
+    % Can now compute the standard deviation consistent with this as
+    % sigma = sqrt(-x^2 / (2 ln(r))  where r is ratio and x is lambda_sd (distance from mean)
+    new_cov_inflate_sd = sqrt( - 0.5 * lambda_sd^2 / log(ratio) );
+
+    % Prevent an increase in the sd of lambda
+    if new_cov_inflate_sd > lambda_sd, new_cov_inflate_sd = lambda_sd; end
+    if new_cov_inflate_sd < lambda_sd_LB, new_cov_inflate_sd = lambda_sd_LB; end
+end
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$

DART_LAB/matlab/run_lorenz_96.m

@@ -891,7 +891,7 @@ function step_ahead()
                         obs_prior, obs_increments);
 
                     % Compute distance between obs and state for localization
-                    dist = abs(i - j) / 40;
+                    dist = abs(i - j) / MODEL_SIZE;
                     if(dist > 0.5), dist = 1 - dist; end
 
                     % Compute the localization factor