%% Different noise level

% clear; clear all; close all; clc;
function differentNoiseLevel_globalEstimation(pointsPerAxis)

%% Define Process and 'measured' training data
% pointsPerAxis   = 51; % 51, 301
uAxis           = linspace(0,1,pointsPerAxis)';
% uAxis           = rand(pointsPerAxis,1);
U               = uAxis;
yUndisturbed    = benchFunktionND(U,'hyperbelGauss');
noiseLevel      = [0 0.05 0.1];
noiseLevel      = [0.05 0.1 0.2];

% Define validation data
deltaU          = (uAxis(2)-uAxis(1))/2;    % 'Zwischengitterstellen' + andere Rauschauspraegung
% deltaU          = uAxis(2)-uAxis(1);        % Gleiche Stellen, andere Rauschauspraegung
uAxisVali       = deltaU:deltaU*2:max(uAxis)-deltaU;
% uAxisVali       = rand(pointsPerAxis,1);
Uvali           = uAxisVali';
yUndisturbedVali = benchFunktionND(Uvali,'hyperbelGauss');

% Predefine cell array, that will contain the local model networks for
% different noise levels
lmnGlobal             = cell(size(noiseLevel,2),1);

%% Plot options
scrsize         = get(0,'ScreenSize');
width           = 500;
height          = 400;
pos             = [(scrsize(3)-width)/2 (scrsize(4)-height)/2 width height];
linewidth       = 2;
fontname        = 'Times';
fontsize        = 14;
colors          = {[0 0 1],[0 0.8 0],[1 0 0]};

figure('Position',pos);
ahLeft = axes;
figure('Position',pos+[width 0 0 0]);
ahRight = axes;
hold(ahLeft,'on');
hold(ahRight,'on');
grid(ahLeft,'on');
grid(ahRight,'on');
% options.plotAxis = ah;

%% Training und Visualisierung
for ii=1:size(noiseLevel,2)
    
    % Add noise to training data
    rng(24,'twister');
    yMeasured       = yUndisturbed + noiseLevel(ii)*randn(size(yUndisturbed,1),1);
    
    % Add noise to validation data
    yMeasuredVali   = yUndisturbedVali + 0*noiseLevel(ii)*randn(size(yUndisturbedVali,1),1);
    
    % Create LOLIMOT model
    lmnGlobal{ii,1}             = lolimot;
    lmnGlobal{ii,1}.xRegressorDegree = 1;
    lmnGlobal{ii,1}.lossFunctionGlobal = 'RMSE';
    
    % Pass training data
    lmnGlobal{ii,1}.input       = U;
    lmnGlobal{ii,1}.output      = yMeasured;
    
    % Pass validation data
    lmnGlobal{ii,1}.validationInput = Uvali;
    lmnGlobal{ii,1}.validationOutput = yMeasuredVali;
    
    % Pass test data
    lmnGlobal{ii,1}.testInput = Uvali;
    lmnGlobal{ii,1}.testOutput = yMeasuredVali;
    
    % Some options for the training
    lmnGlobal{ii,1}.smoothness  = 1.0;
    lmnGlobal{ii,1}.maxNumberOfLM = 15;
    lmnGlobal{ii,1}.minError    = 0.00;
    lmnGlobal{ii,1}.maxValidationDeterioration = inf;
    lmnGlobal{ii,1}.maxPenaltyDeterioration = inf;
    
    % Start the LOLIMOT training algorithm
    lmnGlobal{ii,1}             = lmnGlobal{ii,1}.train;
    
    GEglobalLossFun_train = NaN(1,size(lmnGlobal{ii,1}.history.globalLossFunction,2));
    GEglobalLossFun_vali = NaN(1,size(lmnGlobal{ii,1}.history.globalLossFunction,2));
    for jj=1:size(lmnGlobal{ii,1}.history.globalLossFunction,2)
        % Set model complexity
        lmnGlobal{ii,1}.leafModels = lmnGlobal{ii,1}.history.leafModelIter{jj};
        
        % Estimate parameters with current net configuration globally
        lmnGlobal{ii,1} = lmnGlobal{ii,1}.estimateParametersGlobal;
        lmnTmp = lmnGlobal{ii,1};
        
        % Calculate model output on training data
        yHat_train = lmnTmp.calculateModelOutput(lmnTmp.unscaledInput);
        
        % Calculate model output on validation data
        yHat_vali = lmnTmp.calculateModelOutput(lmnTmp.validationInput);
        
        % Calculate loss functions
        GEglobalLossFun_train(1,jj) = lmnTmp.calcGlobalLossFunction(...
            lmnTmp.unscaledOutput,yHat_train);
        GEglobalLossFun_vali(1,jj) = lmnTmp.calcGlobalLossFunction(...
            lmnTmp.validationOutput,yHat_vali);
    end
    
    % Global loss function (training data)
    plot(ahLeft,GEglobalLossFun_train,...
        'LineWidth',linewidth,'Color',colors{ii});
    
    % Global loss function (validation data)
    plot(ahRight,GEglobalLossFun_vali,...
        'LineWidth',linewidth,'Color',colors{ii});

end

hlLeft = legend(ahLeft,{['\sigma_n = ',num2str(noiseLevel(1))],...
    ['\sigma_n = ',num2str(noiseLevel(2))],['\sigma_n = ',num2str(noiseLevel(3))]});
hlRight = legend(ahRight,{['\sigma_n = ',num2str(noiseLevel(1))],...
    ['\sigma_n = ',num2str(noiseLevel(2))],['\sigma_n = ',num2str(noiseLevel(3))]});
hxLeft = xlabel(ahLeft,'Model Complexity');
hyLeft = ylabel(ahLeft,'Loss Function Value');
htLeft = title(ahLeft,'Training Data (global)');
hxRight = xlabel(ahRight,'Model Complexity');
hyRight = ylabel(ahRight,'Loss Function Value');
htRight = title(ahRight,'Validation Data (Noise Free - global)');

set([hlLeft,hlRight,ahLeft,ahRight,hxLeft,hyLeft,htLeft,hxRight,hyRight,htRight],...
    'FontSize',fontsize,'FontName',fontname);
set([hlLeft,hlRight],'Location','best');

set([ahLeft,ahRight],'YScale','log');
% set([ahLeft,ahRight],'YLim',[0.003 1]);


