%% Lokale vs. Globale Schaetzung

%% Different noise level

clear; clear all; close all; clc;

%% Define Process and 'measured' training data
pointsPerAxis   = 41;
uAxis           = linspace(0,1,pointsPerAxis)';
U               = makeBenchInputND(1,uAxis);
yUndisturbed    = benchFunktionND(U,'hyperbel');
noiseLevel      = [0    0.05    0.05];
smoothness      = [0.1  0.1     1];

% 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;
Uvali           = makeBenchInputND(1,uAxisVali');
yUndisturbedVali = benchFunktionND(Uvali,'hyperbel');

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

%% Plot options
scrsize         = get(0,'ScreenSize');
width           = 600;
height          = 500;
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]};

ah = cell(2,2,size(noiseLevel,2));
% 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
    lmnLocal{ii,1}             = lolimot;
    
    % Pass training data
    lmnLocal{ii,1}.input       = U;
    lmnLocal{ii,1}.output      = yMeasured;
    
    % Pass validation data
    lmnLocal{ii,1}.validationInput = Uvali;
    lmnLocal{ii,1}.validationOutput = yMeasuredVali;
    
    % Some options for the training
    lmnLocal{ii,1}.smoothness  = smoothness(ii);
    lmnLocal{ii,1}.maxNumberOfLM = 6;
    lmnLocal{ii,1}.minError    = 0.00;
    lmnLocal{ii,1}.maxValidationDeterioration = inf;
    lmnLocal{ii,1}.maxPenaltyDeterioration = inf;
    
    % Start the LOLIMOT training algorithm (local)
    lmnLocal{ii,1}             = lmnLocal{ii,1}.train;
    
    % Start the LOLIMOT training algorithm (global)
    lmnGlobal{ii,1}             = lmnLocal{ii,1}.estimateParametersGlobal;
    
    % Create new figure for the current noise level
    figure('Position',pos);
    
    % Create all axes for the current figure
    ah{1,1,ii} = subplot(2,2,1);
    ah{1,2,ii} = subplot(2,2,2);
    ah{2,1,ii} = subplot(2,2,3);
    ah{2,2,ii} = subplot(2,2,4);
    hold(ah{1,1,ii},'on');
    hold(ah{1,2,ii},'on');
    hold(ah{2,1,ii},'on');
    hold(ah{2,2,ii},'on');
    
    figure('Position',pos-[width 0 0 0]);
    ah2 = subplot(2,1,1);
    ah3 = subplot(2,1,2);
    options.plotAxis = ah3;
    hold(ah2,'on');
    
    % Plot validity function and measured data
    plot(ah2,U,yUndisturbed,'LineWidth',linewidth,...
        'Color',[0 0 0]);
    plot(ah2,U,yMeasured,'LineStyle','none','LineWidth',linewidth,...
        'Color',[1 0 0],'Marker','x','MarkerSize',5);
    lmnLocal{ii,1}.plotPartition(1,[],options);
    
    % Plot undisturbed process with the resulting model
    % Global estimation
    plot(ah{1,1,ii},U,yUndisturbed,...
        'LineWidth',linewidth,'Color',[0.3 0.3 0.3]);
    plot(ah{1,1,ii},lmnGlobal{ii,1}.input,lmnGlobal{ii,1}.outputModel,...
        'LineWidth',linewidth,'Color',[1 0 0]);
    visualizeSeperateLocalModels(lmnGlobal{ii,1},[],ah{1,2,ii});
    
    % Local estimation
    plot(ah{2,1,ii},U,yUndisturbed,...
        'LineWidth',linewidth,'Color',[0.3 0.3 0.3]);
    plot(ah{2,1,ii},lmnLocal{ii,1}.input,lmnLocal{ii,1}.outputModel,...
        'LineWidth',linewidth,'Color',[1 0 0]);
    visualizeSeperateLocalModels(lmnLocal{ii,1},[],ah{2,2,ii});
    
    ht1 = title(ah{1,1,ii},'Process and Model (Global estimation)');
    ht2 = title(ah{1,2,ii},'Local Models (Global estimation)');
    ht3 = title(ah{2,1,ii},'Process and Model (Local estimation)');
    ht4 = title(ah{2,2,ii},'Local Models (Local estimation)');
    ht5 = title(ah2,'Process and measured data');
    ht6 = title(ah3,'Validity Functions');
    
    set([ah{1,1,ii},ah{1,2,ii},ah{2,1,ii},ah{2,2,ii},ah2,ah3,ht1,ht2,ht3,ht4,ht5,ht6],...
        'FontSize',fontsize,'FontName',fontname);
    
    pause
    
end




