sbeach_toolbox/fsbeach_modelskill.m

function [bss, rmse, dVc, dVm, ddV] = fsbeach_modelskill(measured, computed, initial)
%edited by Joshua Simmons 04/04/2016 to revise definitions of r2 etc
%XB_SKILL  Computes a variety skill scores
%
%   Computes a variety skill scores: R^2, Sci, Relative bias, Brier Skill
%   Score. Special feature: within the XBeach testbed, the results are
%   stored to be able to show the development of the different skill scores
%   in time.
%
%   Syntax:
%   [r2 sci relbias bss] = xb_skill(measured, computed, varargin)
%
%   Input:
%   measured  = Measured data where the first column contains independent
%               values and the second column contains dependent values
%   computed  = Computed data where the first column contains independent
%               values and the second column contains dependent values
%   initial   = Initial data where the first column contains independent 
%               values and the second column contains dependent values
%   varargin  = var:    Name of the variable that is supplied
%
%   Output:
%   bss       = Brier Skill Score
%   rb20      = Regression through origin R2 value 
%
%   Example
%   [r2 sci relbias bss] = xb_skill(measured, computed)
%   [r2 sci relbias bss] = xb_skill(measured, computed, 'var', 'zb')
%
%   See also xb_plot_skill

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2011 Deltares
%       Bas Hoonhout
%
%       bas.hoonhout@deltares.nl	
%
%       Rotterdamseweg 185
%       2629HD Delft
%
%   This library is free software: you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation, either
%   version 2.1 of the License, or (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://OpenEarth.nl">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and 
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute 
% your own tools.

%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 13 Apr 2011
% Created with Matlab version: 7.9.0.529 (R2009b)

% $Id: xb_skill.m 4725 2011-06-27 12:58:48Z hoonhout $
% $Date: 2011-06-27 14:58:48 +0200 (ma, 27 jun 2011) $
% $Author: hoonhout $
% $Revision: 4725 $
% $HeadURL: https://repos.deltares.nl/repos/OpenEarthTools/trunk/matlab/applications/xbeach/xb_analysis/xb_skill.m $
% $Keywords: $

%% read options

if ~exist('initial', 'var'); initial = []; end;

%% xbeach bugfix    
% bug - sometime xbeach has suprious points on the end (a flat end to
% the profile) this is fake so remove it.
if length(computed)>=10
    if nanmean(computed(1:10,2)) < nanmean(computed(end-9:end,2))
        %xbeach direction
        spr_len = length(diff(computed(:,2))) - find(diff(computed(:,2))~=0,1,'last');
        computed(end-spr_len+1:end,2) = NaN;
    else
        spr_len = find(diff(computed(:,2))~=0,1,'first');
        computed(1:spr_len-1,2) = NaN;
    end
    if spr_len > 2
        disp(['Warning check computed points, ' num2str(spr_len) ' points found to be spurious'] )
    end    
end

%% remove nans

measured = measured(~any(isnan(measured), 2),:);
computed = computed(~any(isnan(computed), 2),:);

if ~isempty(initial)
    initial  = initial (~any(isnan(initial ), 2),:);
end

%% compute skills

[bss rmse dVc dVm ddV] = deal(nan);

if size(measured,1) > 1 && size(computed,1) > 1
    
    x       = unique([computed(:,1) ; measured(:,1)]);
    x       = x(~isnan(x));

    zmt     = interp1(measured(:,1), measured(:,2), x);    
    zct     = interp1(computed(:,1), computed(:,2), x);

    % determine active zone
    if ~isempty(initial)
        zit = interp1(initial(:,1),  initial(:,2),  x);
        
        %average over 10 steps for small spikes
        dzmi = smooth(abs(zmt-zit),10);
        dzmi((isnan(zmt)|isnan(zit))) = 0;
        dzcm = smooth(abs(zct-zmt),10);
        dzcm((isnan(zct)|isnan(zmt))) = 0;
        
        dz  = max(dzmi,dzcm);
        %has to be more change than a standard measurement error (10cm say) 
        zi1 = find(dz>0.1,1,'first');
        zi2 = find(dz>0.1,1,'last');

        x   = x  (zi1:zi2);
        zmt = zmt(zi1:zi2);
        zct = zct(zi1:zi2);
        zit = zit(zi1:zi2);
    end
    
    %% clear all nan values in input (to avoid problems with interpolation)
    %all equal at this point and assign values for BSS calc
    x0 = x;
    xc = x;
    xm = x;
    z0 = zit;
    zc = zct;
    zm = zmt;
    
    x0(isnan(z0)) = [];
    z0(isnan(z0)) = [];
    xc(isnan(zc)) = [];
    zc(isnan(zc)) = [];
    xm(isnan(zm)) = [];
    zm(isnan(zm)) = [];

    %% determine new grid covered by all profiles
    x_new = unique([x0' xc' xm']);

    %% interpolate profiles onto new grid
    [x0,id1] = unique(x0);
    [xc,id2] = unique(xc);
    [xm,id3] = unique(xm);
    z0 = z0(id1);
    zc = zc(id2);
    zm = zm(id3);

    z0_new = interp1(x0, z0, x_new);
    zc_new = interp1(xc, zc, x_new);
    zm_new = interp1(xm, zm, x_new);

    % reduce vectors if any nans are left
    id = ~(isnan(z0_new) | isnan(zc_new) | isnan(zm_new));
    [x_new z0_new zc_new zm_new] = deal(x_new(id), z0_new(id), zc_new(id), zm_new(id));

    %% Calculate volume change
    try
        Vc = volumefind(x_new,zc_new,0);
        Vm = volumefind(x_new,zm_new,0);
        V0 = volumefind(x_new,z0_new,0);
        
        dVm = Vm-V0; %measured volume change
        dVc = Vc-V0; %computed volume change
        ddV = dVc-dVm; %predicted volume change - measured
        nddV = ddV/abs(V0-Vm);
    catch
        dVm = NaN;
        dVc = NaN;
        ddV = NaN;
        nddV = NaN;
    end
    %% derive weight
    % weight is defined as half of diff left and half of diff right of each
    % point; first and last point are considered to have only one side.
    weight = diff([x_new(1) diff(x_new)/2 + x_new(1:end-1) x_new(end)]);
    weight = weight / sum(weight);

    %% calculate BSS only
    % brier skill score
    if ~isempty(initial)
        mse_p = sum( (zm_new - zc_new).^2 .* weight );
        mse_0 = sum( (zm_new - z0_new).^2 .* weight );
        bss = 1. - (mse_p/mse_0);
    else
        bss = 1-(std(zc-zm))^2/var(zm);
    end
    
    %% Calculate others
    
    rmse     = sqrt(mean((zc_new-zm_new).^2));
    rmsm    = sqrt(mean(zm_new.^2));
    sci     = rmse/max(rmsm,abs(mean(zm_new)));
    relbias = mean(zc_new-zm_new)/max(rmsm,abs(mean(zm_new)));
    
end