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 .
% --------------------------------------------------------------------
% This tool is part of OpenEarthTools.
% 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
% 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