function [obj thetaAll] = estimateParametersGlobal(obj, estimationProcedure)
% ESTIMATEPARAMETERSGLOBAL estimates the parameters of all local models globally with one single LS
% estimation, if 'LS' is used as estimationProcedure. Furthermore, the parameter estimation can be
% regularized using ridge regression, if 'RIDGE' is used as estimationProcedure. Then, the
% regularization parameter lambda is optimized w.r.t. the Leave-One-Out-score or PRESS-value,
% respectively. The parameter estimates are stored in the object.
%
%
% [obj thetaAll] = estimateParametersGlobal(obj, estimationProcedure)
%
%
% INPUT:
%   estimationProcedure:  (string)  Method for estimation procedure. Either 'LS'
%                                   for global estimation or 'RIDGE' global estimation
%                                   combined with ridge regression.
%
%
% OUTPUT:
%   obj:                  (object)  Object containing newly estimated LM parameters.
%   thetaAll:             (1 x M)   Cell array containing the estimated output
%                                   weights for each local model.
%
%
%
% LMNtool - Local Model Network Toolbox
% Benjamin Hartmann & Torsten Fischer, 11-December-2012
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles



%% Reliability check for input variables

% If no estimationProcedure is explicitly given, use the same as stored in the object
if nargin < 2
    estimationProcedure = obj.estimationProcedure;
end

% If no explicit input arguments are given
try
    xRegressor = obj.xRegressor;
    output     = obj.output;
    if any(strcmp(superclasses(obj),'gaussianOrthoGlobalModel'))
        validityFunctionValue = cell2mat(obj.calculateVFV(obj.MSFValue(obj.leafModels)));
    else
        validityFunctionValue = obj.phi(obj.leafModels);
    end
catch
    error('estimateLocalModel:estimateParametersGlobal','No matrices available for estimation! Check X, y and phi.')
end

% Get some constants
[numberOfSamples numberOfxRegressors] = size(xRegressor);
numberOfOutputs = size(output,2);
M = sum(obj.leafModels);

% Set dataWeighting to default, if not defined already.
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(validityFunctionValue,1) ~= numberOfSamples, size(dataWeighting,1) ~= numberOfSamples])
    error('estimateLocalModel:estimateParametersGlobal','Missmatching number of samples between xRegressor, output, validity or data weighting!')
end
% Check size of validity function matrix.
if all(size(validityFunctionValue) ~= [numberOfSamples M])
    error('estimateLocalModel:estimateParametersGlobal', '<validityFunctionValue> must be a matrix whose number of rows and columns must match with the number of samples and the number of local models.')
end
% Check for complex entries.
if ~isreal(xRegressor) || ~isreal(output) || ~isreal(validityFunctionValue)
    error('estimateLocalModel:estimateParametersLocal','<xRegressor>, <output> and <validityFunctionValue> must be real-valued.');
end
% Check if data is too spares and give warning
if obj.history.displayMode && (numberOfSamples < M*numberOfxRegressors)
    warning('estimateModel:estimateParametersGlobal','Data is actually too sparse for a reliable estimate of the parameters!');
end


%% Global estimation of the local model parameters

% Initialize
thetaAll = cell(1,M);

% Perform data weighting, if any element of dataWeighting has not value 1.
if any(dataWeighting~=1)
    validityFunctionValue = bsxfun(@times,dataWeighting,validityFunctionValue);
end

% Switch between different estimation methods
switch estimationProcedure
    
    case 'LS'  % Ordinary global least-squares estimation.
        
        % Assembly of the weighted regression matrix for global estimation.
        XAll = zeros(numberOfSamples,numberOfxRegressors*sum(obj.leafModels));
        for lm = 1:M
            XAll(:,(lm-1)*numberOfxRegressors + (1:numberOfxRegressors)) = bsxfun(@times,validityFunctionValue(:,lm),xRegressor);  % Global regression matrix
        end
        
        % Global LS estimation for all LM at once
        theta = pinv(XAll)*output; % For numerical robustnes: pseudoinverse using SVD.
        
        % Re-arrange the parameters in the cell thetaAll
        for lm = 1:M
            thetaAll{lm} = theta((lm-1)*numberOfxRegressors + (1:numberOfxRegressors),:);
        end
        
    case 'RIDGE'  % Estimation using SVD. Lambda is optimized w.r.t. the PRESS value / LOO-cross validation.
        
        % Assembly of the weighted regression matrix for global estimation.
        XAll = zeros(numberOfSamples,numberOfxRegressors*M);
        for lm = 1:M
            XAll(:,(lm-1)*numberOfxRegressors + (1:numberOfxRegressors)) = bsxfun(@times,validityFunctionValue(:,lm),xRegressor);  % Global regression matrix
        end
        
        % Produce the singular value decomposition of XAll.
        [U,S,V]=svd(XAll,'econ');
        p = size(S,1);
        
        % Optimiere lambda bezüglich Leave-One-Out-Fehler = PRESS
        myOptions = optimset;
        myOptions = optimset(myOptions,'Display','iter','Diagnostics','off');
        
        % First, find a good initial value for lambda in the bounds [1e-10,1e10]
        [lambdaIni PRESSIni] = fminbnd(@(lambda) obj.findMinPRESSLambda(lambda,U,S,output) ,1e-10,1e10,myOptions)
        
        % Second, perform an unconstrained optimization with the previously found lambdaIni.
        [lambdaOpt PRESSOpt] = fminunc(@(lambda) obj.findMinPRESSLambda(lambda,U,S,output) ,lambdaIni, myOptions)
        
        % Calculate the LM parameters
        theta = V*S*((S.^2 + lambdaOpt*eye(p))\eye(p))*U'*output;
        
        % Re-arrange the parameters in the cell thetaAll
        for lm = 1:M
            thetaAll{lm} = theta((lm-1)*numberOfxRegressors + (1:numberOfxRegressors),:);
        end    
        
                
    otherwise
        % intended estimation prodedure is not implemented?
        error('estimateModel:estimateParametersGlobal','Unknown estimation procedure: use "LS" or "RIDGE"')
end

% Store parameters in the object
leaves = find(obj.leafModels);
for k = 1:M
    obj.localModels(leaves(k)).parameter = thetaAll{k};
end
% Delete all parameters in knots, because they're not valid anymore.
knots = find(~obj.leafModels);
for k = 1:sum(~obj.leafModels)
    obj.localModels(knots(k)).parameter = [];
end

% Calculate model output
obj.outputModel = calculateModelOutput(obj);

% Calculate the local loss function values for all leaf models
localLossFunctionValue = calcLocalLossFunction(obj, obj.unscaledOutput, obj.outputModel , validityFunctionValue, obj.dataWeighting, obj.outputWeighting);
for k = 1:M
    % Allocate the corrected local loss function value to each local model
    obj.localModels(leaves(k)).localLossFunctionValue = localLossFunctionValue(k);
end

% Update the history object
obj.history = writeHistory(obj.history,obj);


