function [ofcc,oy] = gmpe_fcc(ip, varargin)
% function [ofcc,oy] = gmpe_fcc(ip,...)
% Fahlman Cascade Correlation
% Estimation of the Fahmal's cascade correlation network.
% 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)
%   Optional input parameters:
%       'Max levels' maksymalna iloœæ poziomów sieci; default 50
%       'Min progress' minimalny progres R^2; default 0.0
%       'No candidates' iloœæ kanadydatów przy dodawinmiu kolejnego poziomu; default 100 x N
%       'Max fails' maksymalna iloœæ poziomów nie spe³niaj¹cych minimalnego progresu
%       '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) ;
    maxnn = arparameters('Max levels',50,varargin{:}) ;
    minprogres = arparameters('Min progress',0,varargin{:}) ;
    minnn = arparameters('Min levels',minprogres,varargin{:}) ;
    ipsylon_rp = arparameters('Regularization parameter',0,varargin{:}) ;
    NTryes = arparameters('No candidates',ceil(100*sqrt(xM)),varargin{:}) ;
    max_fails = arparameters('Max fails',1,varargin{:}) ;
    nosriter = arparameters('OSR iterations',50,varargin{:}) ;
    regression_method = arparameters('Regression method','LR',varargin{:}) ;
    isprogress = max_fails ;
    p = xN ; 
    cidx = 0 ; % index of levels
    SSE0 = sum(sqr(Y - mean(Y))) / p ;
    oldR2 = 0;
    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 (((isprogress > 0)||(cidx <= minnn)) && (cidx <= maxnn))
        if cidx
            e = wyn_osr.e;
            fprintf('e: %.2e: ',std(e))
            temp_fval= zeros(1,NTryes) ;
            temp_w = cell(1,NTryes) ;
            parfor iii = 1:NTryes
                w0 = (rand(cidx+xM,1)-0.5) ./ sqrt(max(abs(X),[],1))';
                [temp_w{iii},temp_fval(iii)] = fminsearch(@(w)gmpe_fcc_S(X,e,w,ipsylon_rp), w0, options) ;
            end
            [val,iiidx] = min(temp_fval) ;
            w = temp_w{iiidx} ;
            on = tanh(X*w) ;
            fprintf('O: %.2f, %.2f; Val %.4f at %d\n',mean(on),std(on),-val,iiidx) ;
        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 ;
        if (R2-minprogres > oldR2)
            isprogress = max_fails ;
        else
            isprogress = isprogress - 1 ;
        end
        fprintf('%d: SSE = %f(%f, p=%d) R2 = %f(%d); on %f %f\n', cidx, wyn_osr.SSE, SSEp, p, R2, isprogress, max(on), min(on))
        if ((isprogress > 0)||(cidx <= minnn)) 
            ofcc(cidx) = struct('b',wyn_osr, 'w',w, 'X',X, 'R2',R2) ;
        end
        oldR2 = max(oldR2,R2) ;
    end
    if nargout > 1
        oy = sim_gmpe_fcc(ofcc,ip.X) ;
    end
end
