function [outputModel,MSFValue,validityFunctionValue] = simulateParallel(obj,xRegressor,zRegressor)
% SIMULATEPARALLEL simulates the output of a local linear neuro-fuzzy
% model. This is need to be done for dynamic systems, because each time
% step is depending on its antecessor. The outputs are calculated parallel.
%
% [outputModel,MSFValue,validityFunctionValue] = simulateParallel(obj,xRegressor,zRegressor)
%
%
% OUTPUTS:
%
% outputModel:          (N x q)     Matrix of model outputs
%
% MSFValue:             (N x M)     Matrix of (non-normalized) membership
%                                   function values
%                       
% validityFunctionValue (N x M)     Matrix of normalized validity function
%                                   values
%
%
% INPUTS:
%
% obj:                  (object)    local model object containing all
%                                   relevant net and data set information
%                                   and variables
%
% xRegressor:           (N x nx)    regression matrix of the consequents
%
% zRegressor:           (N x nz)    regression matrix of the premises
%
%
% SYMBOLS AND ABBREVIATIONS:
%
% LM:  Local model
%
% p:   Number of inputs (physical inputs)
% q:   Number of outputs
% N:   Number of data samples
% M:   Number of LMs
% nx:  Number of regressors (x)
% nz:  Number of regressors (z)
%
%
% LMNtool - Local Model Network Toolbox
% Tobias Ebert, 18-October-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles


% number of data samples
numberOfSamples = size(xRegressor,1);

% number of outputs
numberOfOutputs = obj.numberOfOutputs;

% determine number of regressor inputs
xNumberOfInputRegressorsSum = sum(cellfun(@length,obj.xInputDelay));
zNumberOfInputRegressorsSum = sum(cellfun(@length,obj.zInputDelay));

% initialize the outputs
outputModel = zeros(numberOfSamples,numberOfOutputs);
MSFValue = cell(1,length(obj.leafModels));
validityFunctionValue = cell(1,length(obj.leafModels));

% get the parameters of all leaf models
%localParameter = cellfun(@(cobj) cobj.parameter,obj.localModels(obj.leafModels),'UniformOutput',false);
%localParameter = arrayfun(@(cobj) cobj.parameter,obj.localModels(obj.leafModels),'UniformOutput',false);
localParameter = {obj.localModels(obj.leafModels).parameter};

for k = 1:numberOfSamples % through all samples
    
    % if offset is last entry in x-regressor
    %idx = xNumberOfInputRegressorsSum + 1; % for regression matrix x
    
    % if offset is first regressor, the xOutputDelays start at index
    % xNumberOfInputRegressorsSum+2! (delayed inputs + Offset + first idx after that)
    idx = xNumberOfInputRegressorsSum + 2; % for regression matrix x
    
    % fill matrix with output regressors
    for out = 1:numberOfOutputs % through all outputs
        for outReg = 1:length(obj.xOutputDelay{out}) % through all output regressors
            kDelay = k-obj.xOutputDelay{out}(outReg);
            if kDelay > 0
                xRegressor(k,idx) = outputModel(kDelay,out);
            end
            idx = idx + 1;
        end
    end
    
    % as there is no offset in z-reressor, the idx with +1 is corret
    idx = zNumberOfInputRegressorsSum + 1; % for regression matrix z
    
    % fill matrix with output regressors
    for out = 1:numberOfOutputs % through all outputs
        for outReg = 1:length(obj.zOutputDelay{out}) % through all output regressors
            kDelay = k-obj.zOutputDelay{out}(outReg);
            if kDelay > 0
                zRegressor(k,idx) = outputModel(kDelay,out);
            end
            idx = idx + 1;
        end
    end
    
    % calculate validity function values, only one row
    % MSFValue1row = cellfun(@(cobj) cobj.calculateMSF(zRegressor(k,:)),obj.localModels(obj.leafModels),'UniformOutput',false);
    MSFValue1row = arrayfun(@(cobj) cobj.calculateMSF(zRegressor(k,:)),obj.localModels(obj.leafModels),'UniformOutput',false);
    validityFunctionValue1row = obj.calculateVFV(MSFValue1row);
    
    % update each leaf model for the calculated sample k
    for kk = 1:sum(obj.leafModels)
        MSFValue{kk}(end+1,:) = MSFValue1row{kk};
        validityFunctionValue{kk}(end+1,1) = validityFunctionValue1row{kk};
    end
    
    % calculate model output for sample k
    outputModel(k,:) = obj.calcYhat(xRegressor(k,:),validityFunctionValue1row,localParameter);
    
end

end

