function [coeff pseudoInv] = estimateParametersLocal(obj, xRegressor, output, validityFunctionValue, numberOfCoeff, dataWeighting)
% ESTIMATEPARAMETERSLOCAL estimates the parameters of all local models. In order to perform local
% estimation, the regressor matrix and the output targets are weighted with their associated
% validity function value. The user has the possibility to incorporate prior knowledge in terms of a
% data weighting vector. This allows a specific weighting of the data, e.g., if the user knows about
% uncertainties in some data regions. If the process has several outputs, the data weighting is a
% matrix with q columns. The estimated parameters are stored in coeff.
%
%
% [coeff] = estimateParametersLocal(obj, xRegressor, output, ...
%   validityFunctionValue, numberOfCoeff, dataWeighting);
%
%
% INPUT:
%   xRegressor:             (N x nx)    Matrix containing the complete
%                                       regression matrix (columns of
%                                       used regressors).
%   output:                 (N x q)     Output vector.
%   validityFunctionValue:  (N x 1)     Array of validity of the local
%                                       model.
%   numberOfCoeff:          (1 x 1)     Indicates how many regressor are
%                                       to be chosen.
%   dataWeighting:          (N x 1)     Weight for each data point; to
%                                       include previous knowledge.
%
%
% OUTPUT:
%   coeff:                  (nx x q)    Array containing the estimated
%                                       coefficient of the polynomial for
%                                       each local model.
%
%
% LMNtool - Local Model Network Toolbox
% Benjamin Hartmann, 04-April-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles

% 16.11.11 added ridge regression (TE)
% 03.02.12 Subset Selection with all data points via Single-step local estimation.
%          Two new cases:  'STEWISE_ALL_LS' and 'REDUCED_STEPWISE_ALL_LS', Benjamin Hartmann.
% 08.02.12 Totally updated stepwise regression procedure that allows subset selection with weighted
%          regressors.


% Get some constants
[numberOfSamples numberOfxRegressors] = size(xRegressor);
numberOfOutputs = size(output,2);

%% Reliability check for input variables

% Set dataWeighting to default if not defined
if ~exist('dataWeighting','var') || isempty(dataWeighting)
    dataWeighting = ones(numberOfSamples,1);
end

% all variables must have the same number of rows corresponding to the
% number of samples
if any([size(output,1) ~= numberOfSamples, size(dataWeighting,1) ~= numberOfSamples])
    error('estimateLocalModel:estimateParametersLocal: missmatch number of samples between xRegressor, output, validity or data weighting!')
end

% Set validityFunctionValue to default if not defined
if ~exist('validityFunctionValue','var') || isempty(validityFunctionValue)
    validityFunctionValue = ones(numberOfSamples,1);
end

% check if the number of validity function values match to the number
% of samples
if all(size(validityFunctionValue) ~= [numberOfSamples 1])
    error('<validityFunctionValue> must be a vector with N entries')
end

% if subset selection is used, check if the number of coefficients is possible to handle
if strcmp(obj.estimationProcedure,'STEPWISE_LS') || strcmp(obj.estimationProcedure,'REDUCED_STEPWISE_LS') || strcmp(obj.estimationProcedure,'SIMPLE_SELECTION') || strcmp(obj.estimationProcedure,'FORWARD_SELECTION')
    if ~exist('numberOfCoeff','var') || ~isscalar(numberOfCoeff) || isempty(numberOfCoeff)
        error('estimateLocalModel:estimateParametersLocal','an incorrect number of coefficients (numberOfCoeff) is given')
    elseif numberOfCoeff > numberOfxRegressors % too many coefficients shall be chosen: reduce the number of coefficients to the number of given Xregressors
        %if obj.history.displayMode
        %    fprintf('\n estimateLocalModel:estimateParametersLocal: the number of coefficient to select is higher than or equal to the number of potential regressors \n the complete regressor is used for estimation! \n')
        %end
        numberOfCoeff = numberOfxRegressors;
    elseif numberOfCoeff < 1 % no coefficients are sought, no proper estimation is possible
        error('estimateLocalModel:estimateParametersLocal','for a proper estimation at least one coefficient has to be chosen!')
    end
end

%% Local estimation of the model parameters

% Initialize the coefficient vector and the pseudo inverse
coeff = zeros(numberOfxRegressors, numberOfOutputs);
pseudoInv = [];

% Weight the x-regressor and the output with sqrt(validityfunctionvalue)
r = sqrt(validityFunctionValue.*dataWeighting);             % Weighting vector for data.
rMat = r(:,ones(numberOfxRegressors,1));                    % Generate an N x nx matrix [r r ... r].
xRegressorWeighted = xRegressor.*rMat;                      % Weighted regression matrix.
rMat = r(:,ones(numberOfOutputs,1));                        % Generate an N x q matrix [r r ... r].
outputWeighted = output.*rMat;                              % Weighted output matrix.

% Option for cutting the data points such that only the points corresponding to the investingated
% local model are left.
cutLMPoints = false;
if cutLMPoints
    idxLMPoints = (validityFunctionValue>=0.5);
    while sum(idxLMPoints) < 2*numberOfCoeff
        phiCut = validityFunctionValue(~idxLMPoints);
        if isempty(phiCut)
            % no additional data samples can be added
            if obj.history.displayMode
                fprintf('\n\nValid data samples: %g, coefficient limit : %g\n',sum(idxLMPoints),2*numberOfCoeff)
                fprintf('\nNo additional data samples to add where found!\n')
            end
            break
        end
        [~,idxMax] = max(phiCut);
        idxCut = (validityFunctionValue==phiCut(idxMax));
        idxLMPoints(idxCut) = true;
        
        % Select lines that correspond to the current LM
        xRegressorWeighted = xRegressorWeighted(idxLMPoints,:);
        outputWeighted     = outputWeighted(idxLMPoints,:);
    end
end

% Switch between different estimation methods
switch obj.estimationProcedure
    
    case 'LS'
        % LS-estimation for the coefficients
        if ~obj.LOOCV
            coeff = xRegressorWeighted\outputWeighted; % using the Q-R-factorization for fast computation
        else
            try
                pseudoInv = pinv(xRegressorWeighted);
                coeff = pseudoInv*outputWeighted;
            catch
                warning('pinv was not able to calculate the pseudo-inverse.')
                coeff = xRegressorWeighted\outputWeighted; % using the Q-R-factorization for fast computation
            end
        end
        %         coeff = robustfit(xRegressorWeighted,outputWeighted,'bisquare',4.685,'off');
        
        
    case {'STEPWISE_LS','REDUCED_STEPWISE_LS','FORWARD_SELECTION'}
        
        % initialize some parameters for the stepwise
        pEnter = [];
        pRemove = [];
        initialModel = [];
        
        % Set customized properties for different subset selection approaches
        switch obj.estimationProcedure
            case 'STEPWISE_LS'
                % All regressor can be selected. Only the offset term is strictly included.
                regressorsToSelect = numberOfCoeff;
            case 'REDUCED_STEPWISE_LS'
                % Ensure, that at least the offset and the linear terms are included in each local model; if more regressors are intended use
                % subset selection.
                initialModel = false(1,size(xRegressorWeighted,2)-1); % the offset is not considered
                initialModel(1:obj.numberOfInputs) = true;
                if numberOfCoeff > obj.numberOfInputs
                    regressorsToSelect = numberOfCoeff-obj.numberOfInputs;
                else
                    regressorsToSelect = 0;
                    %warning('hallo!!!')
                end
            case 'FORWARD_SELECTION'
                % Set pEnter and pRemove such that only forward selection is performed.
                pEnter = 1-eps;
                pRemove = 1-eps; % no removing allowed
                regressorsToSelect = numberOfCoeff;
        end
        
        % Select regressors and estimate parameters for all outputs
        for out = 1:numberOfOutputs
            
            % Perform subset selection
            if regressorsToSelect == 1
                
                % Only apply weighted LS estimation
                coeffReduced = xRegressorWeighted(:,1:numberOfCoeff)\outputWeighted(:,out);
                coeff(1:numberOfCoeff,out) = coeffReduced;
                
            else
                
                if strcmp(version('-release'),'2011b') || strcmp(version('-release'),'2012a')
                    [b,~,~,inmodel,stats] = obj.stepwiseWeighted(...
                        xRegressorWeighted(:,2:end),    outputWeighted(:,out), ...
                        'maxiter',                      regressorsToSelect-1, ...
                        'offsetWeighting',              xRegressorWeighted(:,1), ...
                        'inmodel',                      initialModel, ...
                        'keep',                         initialModel, ...
                        'penter',                       pEnter, ...
                        'premove',                      pRemove, ...
                        'display',                      'off');
                elseif strcmp(version('-release'),'2012a')
                    [b,~,~,inmodel,stats] = obj.stepwiseWeighted(...
                        xRegressorWeighted(:,2:end),    outputWeighted(:,out), ...
                        'maxiter',                      regressorsToSelect-1, ...
                        'offsetWeighting',              xRegressorWeighted(:,1), ...
                        'inmodel',                      initialModel, ...
                        'keep',                         initialModel, ...
                        'penter',                       pEnter, ...
                        'premove',                      pRemove, ...
                        'display',                      'off');
                elseif strcmp(version('-release'),'2012b')
                    [b,~,~,inmodel,stats] = obj.stepwiseWeighted(...
                        xRegressorWeighted(:,2:end),    outputWeighted(:,out), ...
                        'maxiter',                      regressorsToSelect-1, ...
                        'offsetWeighting',              xRegressorWeighted(:,1), ...
                        'inmodel',                      initialModel, ...
                        'keep',                         initialModel, ...
                        'penter',                       pEnter, ...
                        'premove',                      pRemove, ...
                        'display',                      'off');
                elseif strcmp(version('-release'),'2011a')
                    [b,~,~,inmodel,stats] = obj.stepwiseWeighted2011a(...
                        xRegressorWeighted(:,2:end),    outputWeighted(:,out), ...
                        'inmodel',                      initialModel, ...
                        'penter',                       pEnter, ...
                        'premove',                      pRemove, ...
                        'display',                      'off', ...
                        'maxiter',                      regressorsToSelect-1, ...
                        'keep',                         initialModel, ...
                        'scale',                        'off',...
                        'offsetWeighting',              xRegressorWeighted(:,1));
                else
                    error('Other versions are not intended for this class!')
                end
                
                % The vectors 'theat' and 'inmodel' has to be extended, because the offset is treated seperately during subset selection
                coeffReduced = [stats.intercept; b(inmodel)];
                subsetSelectionIndexCoeff = [true, inmodel];
                
                % Update coefficient vector with the selected regressors from stepwisefit
                coeff(subsetSelectionIndexCoeff,out) = coeffReduced;
                
            end
            
            % Information statement
            %if obj.history.displayMode && sum(sum(coeff~=0))==1
            %    fprintf('\n No varibale selection performed. The local model only consists of a constant offset term.\n')
            %end
            
        end
        
    case 'SIMPLE_SELECTION'
        % simply select the first xRegressors depending on the required number of coefficients
        coeffReduced = xRegressorWeighted(:,1:numberOfCoeff)\outputWeighted;
        coeff(1:numberOfCoeff,:) = coeffReduced;
        
    case 'OLS'
        % Apply Orthogonal Least Squares as subset selection algorithm
        for out = 1:numberOfOutputs
            
            % NOTICE: The column order of coeffReduced must correspond to the order in
            % subsetSelectionIndexCoeff!
            [coeffReduced subsetSelectionIndexCoeff] = obj.orthogonalLeastSquares(xRegressorWeighted, outputWeighted(:,out), numberOfCoeff, obj.minErrorImprovement, obj.thresholdLinearDependency);
            
            % Update coefficient vector with the selected regressors form OLS
            coeff(subsetSelectionIndexCoeff,out) = coeffReduced;
        end
        
    case 'RIDGE'
        % Perform ridge regression
        coeff = ridge(outputWeighted,xRegressorWeighted(:,2:end),1e-8,0);
        
    case 'WLS_BACKWARD_ELIMINATION'
        % Perform backward elimination with confidence level alpha. xRegressor must be scaled to unit variance and zero mean!
        coeff = obj.WLSBackwardElimination(xRegressor,output,validityFunctionValue);
        
    otherwise
        error('estimateLocalModel:estimateParametersLocal','unknown estimation procedure: use "LS", "STEPWISE_LS", "RIDGE"')
end

end