function [outputModel,phi] = simulateParallel(obj,xRegressor,zRegressor)
% SIMULATEPARALLEL  Simulate the output of a local linear neuro-fuzzy model
%
% [outputModel,phi] = simulateParallel(obj,xRegressor,zRegressor)
%
% Output:
%   outputModel:    (N x q)    Matrix of model outputs
%   phi:            (N x M)    Matrix of phi (validity function values)
%
% Input:
%   xRegressor
%   zRegressor

% LMNtool - Local Model Network Toolbox
% Tobias Ebert, 24-April-2012
% 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));

% preallocation for function outputs
outputModel = zeros(numberOfSamples,numberOfOutputs);
phi = zeros(numberOfSamples,sum(obj.leafModels));

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

for k = 1:numberOfSamples % Through all samples
    
    % if offset is last entry
    % 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
    
    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
    phi1row = obj.calculatePhi(zRegressor(k,:),obj.leafModels);
    
    % save phi(k,:) to variable for function output
    phi(k,:) = phi1row;
    
    % Calculate model output for sample k
    outputModel(k,:) = obj.calcYhat(xRegressor(k,:),phi1row,localParameter);
    
end

end

