function [regressor ,exponentMatrix] = data2xRegressorCentered(obj,input,output, numberOfLocalModel)
% DATA2XRERESSOR is used to build the x-regression matrix depending on the
% number of inputs and outputs and their delays. In the x-Regressor matrix
% the offset is intended to be in the last column, i.e. the whole last
% column contains ones.
%
% [regressor, exponentMatrix] = data2xRegressor(input,output)
%
% Input:
%   obj                 object      regressor object containing all
%                                   relevant information and variables
%
%   input               (N x p)     data matrix containing physical inputs
%
%   output              (N x q)     data matrix containing physical
%                                   outputs (optional)
%
%   numberOfLocalModel  (1 x 1)     idicates for which local model the
%                                   transformation should be done
%
% Output:
%   regressor           (N x nx)    regression matrix containing inputs
%                                   (and outputs if dynamic model)
%
%   exponentMatrix      (nx x p)    exponent matrix with as many coloums as
%                                   inputs
%
%
% 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, 24-April-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles

% 30.09.11 v0.1
% 04.10.11 v0.2
% 05.10.11
% 08.11.11 case 'sparsePolynomial' implemented
% 18.11.11  help updated

if ~exist('output','var')
    % only inputs are used to assemble to regressor
    output = [];
end

% create datamatrix for regressor
if ~isempty(input)
    % test if an input is given
    
    if ~isempty(obj.xInputDelay)
        % use given delay
        inputDelay = obj.xInputDelay;
    else
        % use default delays
        inputDelay = num2cell(zeros(1,size(input,2)));
    end
else
    inputDelay = [];
end

if ~isempty(output)
    % test if an output is given
    
    if ~isempty(obj.xOutputDelay)
        % use given delay
        outputDelay = obj.xOutputDelay;
    else
        % use default delays
        outputDelay = cell(1,size(output,2));
    end
else
    outputDelay = [];
end

% assemble all the delayed physical inputs and outputs to a data matrix
data = obj.delay2DataMatrix(input,output,inputDelay,outputDelay);

% Scale the data into local coordinates
% Initialize scaled data matrix
dataScaled = zeros(size(data));
% copy centers of the z space
center = obj.localModels(numberOfLocalModel).center;
    
for idx = 1:length(obj.relationZtoX)% For-loop over all physical inputs and outputs
    
    if obj.relationZtoX(idx) == 0
        dataScaled(:,idx) = data(:,idx) - mean(data(:,idx));
    elseif obj.relationZtoX(idx) > 0;
        dataScaled(:,idx) = data(:,idx) - center(obj.relationZtoX(idx));
    else
        error('regressor:data2xRegressor','The property realtionZtoX must have non negative entries!')
    end
end

% get N and number of inputs
[N, numInputs] = size(dataScaled);

% check exponent matrix
if isempty(obj.xRegressorExponentMatrix)
    % build exponent matrix
    exponentMatrix = obj.buildExponentMatrix(obj.xRegressorType,obj.xRegressorDegree,numInputs);
    
    % delete higher orders of exponents if necessary
    idx2Delete = any(exponentMatrix > obj.xRegressorMaxPower,2);
    exponentMatrix = exponentMatrix(~idx2Delete,:);
    
    % sort the exponent Matrix, according to lowest order to highest
    % calculate the sum of all exponents
    sumExponents = sum(exponentMatrix,2);
    % get highest power for all expoents of each regressor
    maxExponents = max(exponentMatrix,[],2);
    % initialise vector for the new indice sequence to use
    idxOrder = [];
    for k = unique(sumExponents') % loop over sum of exponents
        % get the indices of the regressor to sort for, start with highest
        idx2Sort = (sumExponents == k);
        for kk = unique(maxExponents(idx2Sort))' % loop over max of exponents
            idxOrder = [idxOrder; find((maxExponents == kk) & idx2Sort)];
        end
    end
    % use the idxOrder vector to shape expoentMatrix in the new Order
    exponentMatrix = exponentMatrix(idxOrder,:);
else
    exponentMatrix = obj.xRegressorExponentMatrix;
end

% fill regressor with data
switch obj.xRegressorType
    case {'polynomial', 'sparsePolynomial'}
        if all(all(exponentMatrix== 0))
            % 0th Order Model
            regressor = ones(N,1);
        elseif all(exponentMatrix(exponentMatrix>0)==1)
            % 1st order model, faster than the loop below
            noOfRegressors = size(exponentMatrix,1)+1;
            regressor = ones(N,noOfRegressors);
            regressor(:,2:end) = dataScaled;
        else
            % n-th order model
            noOfRegressors = size(exponentMatrix,1);
            regressor = ones(N,noOfRegressors);
            for h = 1:size(exponentMatrix,2)
                % find all exponent greater than zero.
                % with regressor = ones(...)
                i = exponentMatrix(:,h)>0;
                % multiply inputs
                regressor(:,i) = regressor(:,i) .*  bsxfun(@power,dataScaled(:,h) *ones(1,sum(i)),exponentMatrix(i,h)');
            end
            % add 0th Order
            regressor = [ones(N,1), regressor];
        end
    otherwise
        error('regressor:data2xRegressor','This type of trail function for the regressor is not implemented!')
end



end

