function [mse,times] = genmse(inputs,outputs,scales,model,type,P,reali,scalei,verbose)
%GENMSE Generates normalized MSE curves of a model with homogeneous components
%
%   [MSE,TIMES] = GENMSE(INPUTS,OUTPUTS,SCALES,MODEL,'TYPE',P) returns
%   normalized MSE curves and execution times for model MODEL of type VF
%   (Volterra filter), PARAFAC-FIR, or PARAFAC-IIR up to homogeneous
%   component of order P.
%
%   INPUTS is a cell structure with zero mean and unit variance realizations of
%   the  modelled  system iput
%   OUTPUTS{r}{s} is the modelled system output for realization r and scale
%   factor SCALES(s) of the input.
%
%   [MSE,TIMES] = GENFIT(INPUTS,OUTPUTS,SCALES,MODEL,'TYPE',P,REALI,SCALEI
%   VERBOSE) uses only realizations with indexes REALI and scales SCALES(SCALEI)
%   REALI=[] or SCALEI=[] defaults to all realizations or scales, respectively
%
%   See also VFDIAGS, VCPD
%
%   CBVK Toolbox
%   JHMG, PMSB 2016


% Defaults
nRealiz = length(inputs);
nScales = length(scales);

if nargin < 7
    reali = (1:nRealiz);
    scalei = (1:nScales);
    verbose = false;
else
    if ~isempty(reali)
        nRealiz = length(reali);
    else
        reali = (1:nRealiz);
    end
    if ~isempty(scalei)
        nScales = length(scalei);
    else
        scalei = (1:nScales);
    end
    
    if nargin < 9
        verbose = false;
    end
end


% Pre-computation of model output

sims = cell(nRealiz,1);
times = zeros(nRealiz,P);
tim=cell(nRealiz,1);
yest=cell(nRealiz,1);

try
    if matlabpool('size')==0
      matlabpool
    end
end


if verbose
    disp('Pre-computing the model output...');
end;

parfor r = 1:nRealiz
    
    u0 = inputs{reali(r)};
    switch upper(type)
        case 'VF'       % Volterra filter (in diagonal coordinates representation)
            [yest{r}, tim{r}] = vfdiags(u0',model,length(model{1})*ones(1,length(model)));
            
        case 'SPIVF'    %  VF with sparse-interpolated diagonals
            [yest{r}, tim{r}] = vfdiags(u, model.kernels, model.widths, model.spiopts);
            
        case 'CPD'      % VF-CPD model
            [yest{r}, tim{r}] = vcpd(u0,model);            
            
%         case 'PARAFAC-IIR'   %%%  --> Not yet implemented.
%             [yest, tim] = compiirPARAFAC(u0,model);

        otherwise
            error(['Unrecognized model type "' type '"'])
    end;
    
    % Stores computed output and spent times
    sims{r} = yest{r};
    times(r,:) = tim{r};
    
    if verbose
        disp(['Output for realization ' num2str(r) ' computed.']);
    end;
end;


% Computation of outputs for all scales and fitness curves generation

fit = zeros(nScales,P);

if verbose
    disp('Computing fitness ...');
end;

for r = 1:nRealiz
    
    for n=1:nScales
        
        % Modelled system output
        y = outputs{reali(r)}{scalei(n)};
        
        % Fitness calculation
        yest = sims{r};
        
        for p=1:P
            % yest(p,:) = (scales(scalei(n))^p)*yest(p,:);
            % fit(n,p) = fit(n,p) + compfit(y,sum(yest(1:p,:),1));
            yest(:,p) = (scales(scalei(n))^p)*yest(:,p);
            fit(n,p) = fit(n,p) + compfit(y,sum(yest(:,1:p),2));
        end;
        
        
    end;
    if verbose
        disp(['Run = ' num2str(r) ' computed.']);
    end;
end;

fit = fit / nRealiz;
mse = (1 - fit/100).^2;
