function [CV,wyn_all,Err] = gmpe_fcc_CV(ip,varargin)
% function [CV,GMPE,Err] = gmpe_fcc_CV(ip,...)
% Ground Motion Preduction Equation
% Fahlman Cascade Correlation
% Calculates Cross-Validation (CV) or Nonlinear Cross-Validation (NCV)
% Instead of removing one sample as in classical CV,
% all samples of each event are removed.
% Inputs:
% ip - input parameters
%   ip.data - data used for regression
%   ip.event_key - unique values to recognise the same event
%   ip.outputs_set - value of ground motion
%   ip.inputs_set - GMPE inputs (excluding bias = 1)
% Input options:
%   'nonlinear' - first the model is estimated from all data
%                 and it is start point to CV
%                 It is called nonlinear cross validation (NCV) (Moody, 1994) 
%   'FCC'       - instead of computing inital FCC e.g. for the NCV
%                 this parameter are taken           
%   other options see gmpe_fcc
% Outputs:
%   CV    - value for cross validation prediction risk
%           for all units od the FCC
%   GMPE  - the GMPE FCC model
%   Err   - residuals for all steps of the CV
%
% (c) Jan Wiszniowki 2016.09.25

    p = ip;
    p.data = sortrows(p.data,p.event_key);
    N = length(p.data) ;
    [~,ia,ic] = unique(p.data,p.event_key);
    p.ic = ic;
    p.no_events = length(ia);
    wyn_all = arparameters('FCC',[],varargin{:});
    if isempty(wyn_all)
        fprintf('Compute FCC\n')
        wyn_all = gmpe_fcc(p,varargin{:});
    end
    Nlevels = length(wyn_all);
    SSE = zeros(Nlevels,length(p.outputs_set));
    if nargout > 2
        Err = NaN(length(ic),Nlevels);
    end
    for i = 1:p.no_events
        idxes = (ic == i);
        fprintf('================compute %4d / %4d events; %d/%d records ================\n',i,p.no_events, sum(idxes),N)
        p1 = p;
        p1.data(idxes,:)=[];
        if aroptions('nonlinear',varargin{:})
            wyn1 = gmpe_fcc_1(p1,wyn_all,varargin{:});
        else
            wyn1 = gmpe_fcc(p1,'Max levels',Nlevels,'Min levels',Nlevels,varargin{:});
        end
        x1 = double(p.data(idxes,p.inputs_set));
        y1 = double(p.data(idxes,p.outputs_set));
        for level = 1:Nlevels
            [oy] = sim_gmpe_fcc(wyn1(1:level),x1);
            if nargout > 2
                Err(idxes,level) = (oy - y1);
            end
            SSE(level,:) = SSE(level,:) + sum((oy - y1).^2,1);
        end
%         SSE'
    end
    CV = SSE ./ N ;
end