clear all; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Traccia due tipi di funzione di attivazione: 
% gaussiana (radbas) e triangolare (tribas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
x = linspace(-3,3,250);
subplot(2,2,1)
plot(x,radbas(x))
% Attenzione perche' radbas non coincide esattamente 
% con la funzione di attivazione gaussiana specificata
% nei lucidi
subplot(2,2,2)
x = linspace(-3,3,250);
plot(x,tribas(x))

subplot(2,1,2)
x = linspace(-4,6,250);
mu1 = 0; sigma1= 1;
y1 = exp(-(x-mu1).^2/(sigma1^2));
mu2 = 2; sigma2= 2;
y2 = exp(-(x-mu2).^2/(sigma2^2));
plot(x,y1,x,y2)
mu3 = -2; sigma3= .25;
y3 = exp(-(x-mu3).^2/(sigma3^2));
plot(x,y1,x,y2,x,y3)
axis([-4 6 0 1.1])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo della distanza di Minkowski 
% di ordine r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Paramentro di localizzazione mu_k
mu_k = [3 -2 5 6];
% Vettore di input
x = [1 2 -4 5];
% Ordine
r = 3;
% Metodo diretto
dr = (sum(abs((x-mu_k).^r)))^(1/r)
% Uso del comando norm
dr = norm(x-mu_k, r)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Imposta i parametri del modello generatore
% dei dati
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
randn('seed',12345)
rand('seed',12345)
N=400;       % Numerosità campionaria
ntrain = floor(2*N/3); % Numerosità del training set
ninp = 5;    % No. variabili X di input
devstd = 0.2 % Deviazione standard della componente stocastica

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Genera i dati campionari
% Data generating process: 
% y = 1.2*(x1^3 - 0.5) + sin(3.46*pi*x1) + eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X =rand(N,ninp);                              % Variabili indipendenti
ydet= 1.2*(X(:,1).^3 -.5)+sin(X(:,1)*pi*1.1); % Componente deterministica
eps = devstd*randn(N,1);                       % Componente stocastica
y = ydet + eps;                               % Variabile dipendente y


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preprocessamento dei dati
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verLessThan('matlab', '7.1.0')
    [Xn,minX,maxX,yn,miny,maxy] = premnmx(X',y');
else
    [Xn,psx] = mapminmax(X',-1,1);
    [yn,psy] = mapminmax(y',-1,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Costruisci training and testing set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xtrain = Xn(:,1:ntrain);
ytrain = yn(1:ntrain);
Xtest = Xn(:,ntrain+1:end);
ytest = yn(ntrain+1:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sintassi: net = newrb(X,Y,goal,spread,MN)
% X - Matrice NinpxN di N vettori di input
% Y - Matrice NoutxN di N vettori target
% goal   - Errore medio quadratico da raggiungere
% spread - Parametro di scala sigma
% MaxNeu - Numero massimo di neuroni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
goal = 0;
spread = 1;
MaxNeu = 50;
optnet=newrb(Xtrain,ytrain,goal,spread,MaxNeu,10);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Attenzione newrb provvede anche ad addestrare la rete !
% Il comando "train" qui non serve quindi.
% La stima del modello neurale RBF viene fatta mediante
% un metodo "iterativo":
%
% (1) Data la matrice X delle N osservazioni delle variabili 
% indipendenti, trova le N previsioni di y della rete (yhat)
% e calcola gli N errori commessi dalla rete
%
% (2) Trova il massimo di questi errori (corrispondente alla
% i-esima osservazione Xi) e aggiungi un neurone nello strato
% nascosto avente Xi come parametro di localizzazione
%
% (3) Ristima i pesi delle connessioni fra strato nascosto e
% output
%
% (4) Ripeti finche' l'errore complessivo commesso dalla rete
% scende al di sotto di goal oppure finche' il numero massimo
% di neuroni MaxNeu è stato raggiunto
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Previsione con la rete neurale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yout_train=sim(optnet, Xtrain);
yout_test=sim(optnet, Xtest); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Postprocessamento output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verLessThan('matlab', '7.1.0')
    yhat_train = postmnmx(yout_train,miny,maxy);
    yhat_test  = postmnmx(yout_test ,miny,maxy);
else
    yhat_train = mapminmax('reverse',yout_train,psy);
    yhat_test  = mapminmax('reverse',yout_test ,psy);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rappresentazione grafica della bontà della rete
% di adattarsi ai dati
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
hold on;
X_train = X(1:ntrain,:);
y_train = y(1:ntrain);
ydet_train = ydet(1:ntrain);
[x1sort,order]=sort(X_train(:,1));
plot(x1sort,y_train(order),'b.')
plot(x1sort,yhat_train(order),'-','Color',[0.23 0.62 0.17],'LineWidth',2);
plot(x1sort,ydet_train(order),'r-','LineWidth',1);
box on;
legend('Dati osservati','Previsioni rete','Comp.determinist.',...
       'location','SouthEast');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo della bontà di adattamento 
% ai dati mediante indice U-Theil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UT = 100*sum((y_train-yhat_train').^2)/sum((y_train-mean(y_train)).^2);
title(['U-Theil = ', num2str(UT), '%']);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rappresentazione grafica della 
% capacita' predittiva della rete
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
hold on;
X_test   = X(ntrain+1:end,:);
y_test   = y(ntrain+1:end);
ydet_test = ydet(ntrain+1:end);
[x1sort,order]=sort(X_test(:,1));
plot(x1sort,y_test(order),'b.')
plot(x1sort,yhat_test(order),'-','Color',[0.23 0.62 0.17],'LineWidth',2);
plot(x1sort,ydet_test(order),'r-','LineWidth',1);
box on;
legend('Dati osservati','Previsioni rete','Comp.determinist.',...
       'location','SouthEast');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo della capacità predittiva della rete
% mediante indice U-Theil
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UT = 100*sum((y_test-yhat_test').^2)/sum((y_test-mean(y_test)).^2);
title(['U-Theil = ', num2str(UT), '%']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualizzazione dei pesi della rete
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Strato nascosto
optnet.layers{1}

% Secondo strato
optnet.layers{2}

% Parametri di localizzazione (centri)
mu=optnet.IW{1}
[m,S]=size(mu)

% Pesi fra strato nascosto e output
Omega=optnet.LW{2,1}
size(Omega)



% Marco Sandri
% Università di Brescia
% Dipartimento Metodi Quantitativi
% info@msandri.it