function [ofcc,oy] = gmpe_fcc_1(ip, ifcc, varargin)
% function [ofcc,oy] = gmpe_fcc_1(ip, ifcc,...)
% Fahlman Cascade Correlation
% Modification of the Fahmal's cascade correlation network
%  based on existing network
%  using single candidate neurons. 
% Non clasic computation of correlation coeficient S 
% Counting of the FCC by maximise the correlation of the output of the neuron
% with the residual output error of the network.
% The Pearson product-moment correlation coefficient is applied.
% Version 20160226 - regularization parameter was added
%         20160914 - adding one-stage-regresion and
%                    changes of function parameters
% Inputs:
%   ip - input parameters
%       ip.data - data used for GMPE estimation (dataset)
%       ip.event_key - unique values to recognise the same event (cell of strings)
%       ip.outputs_set - value of ground motion (cell of strings)
%       ip.inputs_set - GMPE inputs (excluding bias = 1) (cell of strings)
%   ifcc - initial FCC
%   Optional input parameters:
%       'Regularization parameter'; default 0 
%       'Regression method'; default 'LR' - linear regresion
%       'OSR iterations' number of iterations of one-stage-regression 
% Outputs
%   ofcc - FCC parameters
%       array of structure coresponding each level of the FCC:
%       b - outputs weights and one step regresion parameters,
%       w - weights of neuronus,
%       p - inputs of one step regresion ,
%       R2 - adjusted coefficient of determination
%   oy - estymated output of the GMPE
%
% (c) Jan Wiszniowki 2016.09.25
    
    X = double(ip.data(:,ip.inputs_set));
    Y = double(ip.data(:,ip.outputs_set));
    [~,ia,ic] = unique(ip.data,ip.event_key);
    nev = length(ia);
    [xN,xM] = size(X) ;
    ipsylon_rp = arparameters('Regularization parameter',0,varargin{:}) ;
    nosriter = arparameters('OSR iterations',50,varargin{:}) ;
    regression_method = arparameters('Regression method','LR',varargin{:}) ;
    p = xN ; 
    cidx = 0 ; % index of levels
    SSE0 = sum(sqr(Y - mean(Y))) / p ;
    wyn_osr = [];
%    fprintf('Start %f (%f, p=%d), max.levels=%d, min.progress=%f\n', SSE0, oldR2, p, maxnn, minprogres)
%    ofcc = struct('b',{}, 'w',{}, 'p',{}, 'R2',{}) ;
    ofcc = struct('b',{}, 'w',{}, 'X',{}, 'R2',{});
    options = optimset('MaxIter', 100000000, 'TolFun', 1e-5, 'TolX',1e-5,'Display','off') ;

    while cidx < length(ifcc)
        if cidx
            e = wyn_osr.e;
%             fprintf('cidx: %.2e: ',cidx)
%             size(ifcc(cidx+1).w)
            w0 = ifcc(cidx+1).w;
%            w0'
            w = fminsearch(@(w)gmpe_fcc_S(X,e,w,ipsylon_rp), w0, options);
%            w'
            on = tanh(X*w) ;
        else
            w = [] ;
            on = ones(xN,1) ;
        end
        X = [X,on] ;
        if strcmp(regression_method,'OSR')
            wyn_osr = one_stage_regression(X,Y,ic,nev,nosriter) ;
        elseif strcmp(regression_method,'LR')
            wyn_osr = gmpe_lr(X, Y) ;
        else
            error('Unknown regression method')
        end
        cidx = cidx+1 ;
        p = p-xM-cidx ;
        SSEp = wyn_osr.SSE/p ;
        R2 = 1 - SSEp/SSE0 ;
        ofcc(cidx) = struct('b',wyn_osr, 'w',w, 'X',X, 'R2',R2) ;
    end
    if nargout > 1
        oy = sim_gmpe_fcc(ofcc,ip.X) ;
    end
end
