%% RBF- vs. NRBF-Netz

clc;
clear;
close all;

%% 
% Einfluss Standardabweichung (sigmas: 0.1, 0.05, 0.2)
% Reaktivierung: Unterschiedliche Standardabweichungen verwenden

%% Definition der Testfunktion
pointsPerAxis   = 41;
uAxis           = linspace(0,1,pointsPerAxis)';
U               = makeBenchInputND(1,uAxis);
yUndisturbed    = benchFunktionND(U,'hyperbel');
noiseLevel      = 0.0;
yTrain          = yUndisturbed + range(yUndisturbed)*noiseLevel*randn(size(yUndisturbed));
extraLinks      = 1;
extraRechts     = 1;

%% Definiere Zentren und Standarabweichungen der radialen Basisfunktionen
% (Gauß-Funktionen) händisch
c               = linspace(0,1,4);
sigma           = 0.1*ones(1,size(c,2));
sigma     = [0.1 0.3 0.1 0.1]

%% Definition der Gauß-Funktion
gaussian        = @(u,center,s) exp(-1/2*((u-center)/s).^2);

%% Visualisierung der radialen Basisfunktionen
fh_basis        = figure;
ah_basis        = axes;
hold(ah_basis,'on');
box(ah_basis,'on');
grid(ah_basis,'on');
uPlot           = linspace(min(U)-extraLinks,max(U)+extraRechts,101)';
gaussianValues  = NaN(size(uPlot,1),size(c,2));

for ii=1:size(c,2)
    gaussianValues(:,ii) = gaussian(uPlot,c(ii),sigma(ii));
%     ph_basis(ii) = plot(ah_basis,uPlot,gaussian(uPlot,c(ii),sigma(ii)));
end
ph_basis        = plot(ah_basis,uPlot,gaussianValues);
title(ah_basis,'Radiale Basisfunktionen');

% Visualisierung der normierten radialen Basisfunktionen
fh_nb           = figure;
ah_nb           = axes;
hold(ah_nb,'on');
box(ah_nb,'on');
grid(ah_nb,'on');
ph_nb           = NaN(1,size(c,2));
nrbf_values     = NaN(size(gaussianValues));

for ii=1:size(c,2)
    nrbf_values(:,ii) = gaussianValues(:,ii)./sum(gaussianValues,2);
end
ph_nrbf         = plot(ah_nb,uPlot,nrbf_values);
title(ah_nb,'Normierte Basisfunktionen');

% Einstellungen für die Visualisierung
set([ah_basis,ah_nb],'YLim',[-0.1 1.1]);

%% Trainig der Gewichte
% RBF Training
% Definition des X-Regressors
Xreg_rbf    = NaN(size(U,1),size(c,2));
for ii=1:size(c,2)
    Xreg_rbf(:,ii) = gaussian(U,c(ii),sigma(ii));
end
w_rbf       = Xreg_rbf\yTrain;

% NRBF Training
% Definition des X-Regressors
sumRBFvalues    = sum(Xreg_rbf,2);
Xreg_nrbf       = Xreg_rbf;
for ii=1:size(c,2)
    Xreg_nrbf(:,ii) = Xreg_rbf(:,ii)./sumRBFvalues;
end
w_nrbf          = Xreg_nrbf\yTrain;

%% Visualisierung RBF- und NRBF-Modell
fh_modelle      = figure;
ah_modelle      = axes;
hold(ah_modelle,'on');
box(ah_modelle,'on');
grid(ah_modelle,'on');

% Erzeuge X-Regressionsmatrizen mit uPlot
X_rbf           = [gaussianValues];
X_nrbf          = [nrbf_values];

% Berechnung der Modellausgänge
y_rbf           = X_rbf * w_rbf;
y_nrbf          = X_nrbf * w_nrbf;

% Berechnung des tatsächlichen Prozessausgangs
uProzess        = linspace(0,1,101)';
yProzess        = benchFunktionND(uProzess,'hyperbel');

ph_prozess      = plot(ah_modelle,uProzess,yProzess,'LineWidth',5,'Color',[0.7 0.7 0.7]);
ph_rbf          = plot(ah_modelle,uPlot,y_rbf,'k-','LineWidth',2);
ph_nrbf         = plot(ah_modelle,uPlot,y_nrbf,'r:','LineWidth',2);

% Extra auskommentiert --> Studenten sollen Modellklasse zuordnen.
% legend(ah_modelle,{'Prozess','RBF-Modell','NRBF-Modell'});

set(ah_modelle,'YLim',[-0.05 2.05],'XLim',[min(U)-extraLinks,max(U)+extraRechts]);