% ------------------------------------------------------------
% Exemplo: 
% 1. Aproximar a superficie de um Vulcão 
% 2. Dividir o modelo em 2 submodelos: montanha + buraco
% ------------------------------------------------------------
clear all;

fprintf('\n (1) Aproximar a superfície VULCÃO \n');

projectname='SplitNN_vulcan';

if (0)
[X,Y]=meshgrid(-1:0.1/3:1);
D=(X.^2+Y.^2);
Z1=exp(-D*4).^2;
Z2=-exp(-D*8).^2;
Z=Z1+Z2;
mZ=max(Z,[],'all');
Z=Z./mZ;
Z1=Z1./mZ;
Z2=Z2./mZ;

[nx,ny]=size(X);
x=[reshape(X,1,nx*ny);reshape(Y,1,nx*ny)];
d1=reshape(Z1,1,nx*ny);
d2=reshape(Z2,1,nx*ny);
d=d1+d2;

% Neural Network (NN)
m=0.5;
niter=100000;

nNE=[20 10]; % Número de neurónios das camadas escondidas
%nNE=50; % Número de neurónios das camadas escondidas
[W,b,E,FN]=my_bptrain(x,d,niter,m,nNE);

y=my_net(x,W,b,FN);

%%
[W,b,E,FN]=my_bptrain(x,d,niter,m,W,b,FN);
% 
fprintf('\n MSE= %d',E(end));
%%
save dados_split1;
%%
else
    load dados_split1;
end

% -----------------------
% IMPRIMIR FIGURAS (jpg)
% -----------------------
figure(1);
surf(X,Y,Z);
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$f(x_1,x_2)$','Interpreter','Latex');

figure(2);
surf(X,Y,Z1);
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$f_1(x_1,x_2)$','Interpreter','Latex');

figure(3);
surf(X,Y,Z2);
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$f_2(x_1,x_2)$','Interpreter','Latex');

figure(4);
surf(X,Y,reshape(y,nx,ny));
zlabel('$y=NN(x_1,x_2)$','Interpreter','Latex');

figure(5);
surf(X,Y,Z-reshape(y,nx,ny));
zlabel('$Error$','Interpreter','Latex')
grid on;
for i=1:5
    figure(i);
    h(1)=xlabel('$x_1$','Interpreter','Latex');
    h(2)=ylabel('$x_2$','Interpreter','Latex');
    set(h,'FontSize',12); %,'FontWeight','bold','FontName','Times New Roman');
end

figure(6);
h=semilogy((1:niter)',E,'LineWidth',2);
h(1)=xlabel('$Iteration$','Interpreter','Latex');
h(2)=ylabel('$MSE$','Interpreter','Latex');
set(h,'FontSize',12); %,'FontWeight','bold','FontName','Times New Roman');
%%
str=['.\figures\',projectname,'_'];
for i=1:6
    strfig=[str,'fig',num2str(i),'.png'];
    print(i, strfig, '-dpng','-r1200');
end

% ======================
% Algoritmo SPLIT NN
% ======================
fprintf('\n (2) Dividir a superfície VULCÃO em 2\n');

% Factor de repartição (inicialmente é atribuido 1 ao modelo existente)
nc=length(W); % n. de camadas (incluindo a de entrada e de saída)
alfW=cell(nc,1);
alfb=cell(nc,1);
for i=1:nc
   alfW{i}=ones(size(W{i}))*0.5; % Vector (linha) da repartição dos pesos dos neurónios da camada i
   alfb{i}=ones(size(b{i}))*0.5; % Vector (linha) da repartição dos pesos dos neurónios da camada i
end

% Backpropagation no ajuste dos parametros SPLIT - alf
%%
niter=50000
m=0.1
[alfW,alfb,E]=splitNNbp(x,d,d1,niter,m,W,b,FN,alfW,alfb);

[ys,as]=my_net_split(x,W,b,FN,alfW,alfb);

%%
figure(7);
surf(X,Y,reshape(ys{1},nx,ny));
zlabel('$NN_1(x_1,x_2)$','Interpreter','Latex');
    
figure(8);
surf(X,Y,reshape(ys{2},nx,ny));
zlabel('$NN_2(x_1,x_2)$','Interpreter','Latex');

figure(9);
surf(X,Y,reshape(ys{1}+ys{2},nx,ny));
zlabel('$NN_1(x_1,x_2)+NN_2(x_1,x_2)$','Interpreter','Latex');

for i=7:9
    figure(i);
    h(1)=xlabel('$x_1$','Interpreter','Latex');
    h(2)=ylabel('$x_2$','Interpreter','Latex');
    set(h,'FontSize',12); %,'FontWeight','bold','FontName','Times New Roman');
end

MSN=mean((y-(ys{1}+ys{2})).^2)
MSN_d1=mean((d1-ys{1}).^2)

figure(10);
h=semilogy((1:niter)',E,'LineWidth',2);
h(1)=xlabel('$Iteration$','Interpreter','Latex');
h(2)=ylabel('$MSE$','Interpreter','Latex');
set(h,'FontSize',12); %,'FontWeight','bold','FontName','Times New Roman');

figure(11);
str1='${{{{\left| {\tilde \alpha } \right|}^T}\left| {\tilde w} \right|} \over {\left\| {\tilde w} \right\|}}$';
str2='$\alpha _{n+1}$';
for i=1:nc
    subplot(2,nc,i);
    sW=sum(abs(W{i}),2);
    bar(sqrt(sum(abs(alfW{i}.*W{i}),2)./sW));
    str3=['W (layer',num2str(i),')'];
    title(['$',str3,'$'],'Interpreter','Latex');
    ylabel(str1,'Interpreter','Latex'); 
    xlabel('$neuron \; number$','Interpreter','Latex');
    subplot(2,nc,i+nc);
    sb=sum(abs(b{i}),2);
    bar(sqrt(sum(abs(alfW{i}),2)./sb));
    str3=['b (layer',num2str(i),')'];
    title(['$',str3,'$'],'Interpreter','Latex');
    ylabel(str2,'Interpreter','Latex');
    xlabel('$neuron \; number$','Interpreter','Latex');
end


figure(12);
subplot(2,3,1);
surf(X,Y,Z1);
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$y_1$','Interpreter','Latex');

subplot(2,3,2);
surf(X,Y,Z2);
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$y_2$','Interpreter','Latex');

subplot(2,3,3);
surf(X,Y,Z);
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$y$','Interpreter','Latex');

subplot(2,3,4);
surf(X,Y,reshape(ys{1},nx,ny));
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$NN_1()$','Interpreter','Latex');

subplot(2,3,5);
surf(X,Y,reshape(ys{2},nx,ny));
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$NN_2()$','Interpreter','Latex');

subplot(2,3,6);
surf(X,Y,reshape(ys{1}+ys{2},nx,ny));
xlabel('$x_1$','Interpreter','Latex');
ylabel('$x_2$','Interpreter','Latex');
zlabel('$NN_1()+NN2()$','Interpreter','Latex');

str=['.\figures\',projectname,'_'];
for i=7:12
    strfig=[str,'fig',num2str(i),'.png'];
    print(i, strfig, '-dpng','-r1200');
end

save dados_split1_final;

% ------------------------------------------------------------
% ------------------------------------------------------------
function [y,yc]=my_net(x,W,b,Fn)
% Função de transferência da Rede Neuronal Multi-camada
% W{i} - pesos que ligam a (i-1)^esima à i^esima camada
% Fn - funções de activação de cada camada da rede (defeito "tansig" para
% as camadas escondidas e "linear" para a última camada).
% x- entrada da rede (nx x n), em que n o números de pontos de entrada)
% y - saída da rede (m x n), com m a dimensão da saída
% Paulo Salgado @ UTAD2020

nc=length(W); % n. de camadas
[nx,n]=size(x); % n. de pontos de treino

% INICIALIZAÇÃO DE PARAMETROS POR DEFEITO
if nargin<4 % Tipo de função de activação por defeito em cada camada
    Fn=cell(nc,1); % Regista o tipo de função de activação
    for i=1:nc-1
        Fn{i}='tansig'; % as camadas escondida usarão a função 'tansig'
    end
    Fn{nc}='purelin'; % a última camada usará função linear
end

yc=cell(nc,1); % Célula com as saídas de cada camada interior da NN
yc{1}=feval(Fn{1},W{1}*x+b{1}); % saida da camada 1
for i=2:nc
    yc{i}=feval(Fn{i},W{i}*yc{i-1}+b{i}); % saida da camada i
end
y=yc{nc}; % Saída da Rede (= saída da última camada)

end

% ------------------------------------------------------------
function [W,b,E,Fn]=my_bptrain(x,d,niter,m,W,b,Fn)
% Algoritmo BackPropagation (BP) no treino da rede neuroral (NN)
% Nota: se a entrada W for um vector de números, tal indica o número de
% neuronios das camadas escondidas. De outro modo será uma célula com os
% pesos da rede.
%
% Paulo Salgado @ UTAD2020

[nx,n]=size(x); % nx - dimensão do vector de entrada; n - n. de dados
[ny,n]=size(d);

nc=length(W); % n. de camadas

% INICIALIZAÇÃO DOS PESOS DA REDE (aletáorio)
% Caso W não seja uma celula será um vector com o número de neuronios da 
% camada intermédia
if ~iscell(W) 
    nc=nc+1; % n. de camadas
    D=W; % D = n. de neurónio da camada (com excepção da última camada que é igual à dimensão do vector de saída
    W=cell(nc,1);
    b=cell(nc,1);
    W{1}=randn(D(1),nx).*0.1;
    b{1}=randn(D(1),1).*0.1;
    for i=2:nc-1
        W{i}=randn(D(i),D(i-1)).*0.1;
        b{i}=randn(D(i),1).*0.1;
    end
    W{nc}=randn(ny,D(nc-1)).*0.1;
    b{nc}=randn(ny,1).*0.1;
end

flag=nargin<7;
if ~flag
    if isempty(Fn)
        flag=true;
    end
end
if flag % Tipo de função de activação por defeito em cada camada
    Fn=cell(nc,1); % Regista tipo de função de activação
    for i=1:nc-1
        Fn{i}='tansig';
    end
    Fn{nc}='purelin';
end

best.W=W;
best.b=b;
best.error=Inf;

% ------------------------------
% Processo iterativo
% ------------------------------
delta=cell(nc,1);

E=zeros(niter,1); % Erro da rede em cada iteração de aprendizagem
    
for ni=1:niter

    [y,yc]=my_net(x,W,b,Fn);

    % Calcular o Gradiente da função custo C face à saída dos neuronio da
    % última camada

    C=0.5*(y-d).^2/n;
    % dC_dy (Gadiente) derivada da função escalar C relativamente ao vector
    % de saída da rede, y
    
    E(ni)=sum(C,2);  % Erro na iteracção "ni"

    if E(ni)< best.error
        best.W=W;
        best.b=b;
        best.error=E(ni);
    end
    % Retropropagar o Erro
    if strcmp(Fn{nc},'purelin')        
        % Método dos minimos quadrados
        A=[yc{nc-1,1};b{nc}.*ones(1,n)];
        wb=d*A'*inv(A*A'+0.01*eye(size(A,1)));
        W{nc}=wb(1:end-1);
        b{nc}=wb(end);
        yc{nc}=feval(Fn{nc},W{nc}*yc{nc-1}+b{nc});
        y=yc{nc};
        C=0.5*(y-d).^2/n;
        Grad_C=(y-d)/n;
        % ...........
        delta{nc}=Grad_C;
        
        nct=nc-1;
    else
        Grad_C=(y-d)/n;
        delta{nc}=((1-yc{nc}).*(1+yc{nc})).*Grad_C;
        nct=nc;
    end

    for i=nc-1:-1:1
        delta{i}=((1-yc{i}).*(1+yc{i})).*(W{i+1}'*delta{i+1});
    end
    
    % -----------------
    % Ajustar os pesos
    % -----------------
    % Camada de entrada
    for j=1:size(yc{1},1)
        DW=sum(delta{1}(j,:).*x,2)';
        W{1}(j,:)=W{1}(j,:)-m*0.5.*DW;
    end
    b{1}=b{1}-m*sum(delta{1},2);

    % restantes camadas
    for i=2:nct
        miu=m*(0.5+1.5*(i-1)/(nc-1));
        for j=1:size(yc{i},1)
            DW=sum(delta{i}(j,:).*yc{i-1},2)';
            W{i}(j,:)=W{i}(j,:)-miu.*DW;
        end
        b{i}=b{i}-miu*sum(delta{i},2);
    end
end

W=best.W;
b=best.b;

end

% ------------------------------------------------------------
% FUNÇÕES SPLIT
% ------------------------------------------------------------
function [y,yc,a]=my_net_split(x,W,b,Fn,alfW,alfb)
% Função de transferência da Rede Neuronal Multi-camada SPLIT
% W{i} - pesos que ligam a (i-1)^esima à i^esima camada
% b{i} - pesos offset das ligações entre a (i-1)^esima à i^esima camada
% alf{i} - matriz de repartição. Tamanho [nn_i x m] , onde "nn_i" é o 
% numero de neurónios da camada i e "m" o n. de termos da serie de Fourier.
% Fn - funções de activação de cada camada da rede (defeito "tansig" para
% as camadas escondidas e "linear" para a última camada).
% x- entrada da rede (nx x n), em que n o números de pontos de entrada)
% y - saída da rede (m x n), com m a dimensão da saída
% Paulo Salgado @ UTAD2020

nc=length(W); % n. de camadas
[nx,n]=size(x); % n. de pontos de treino

% INICIALIZAÇÃO DE PARAMETROS POR DEFEITO
flag=false;
if nargin<4 % Tipo de função de activação por defeito em cada camada
    flag=true;
else
   if isempty(Fn)
       flag=true;
   end
end

if flag
    Fn=cell(nc,1); % Regista o tipo de função de activação
    for i=1:nc-1
        Fn{i}='tansig'; % as camadas escondida usarão a função 'tansig' decomposta em serie de Fourier
    end
    Fn{nc}='purelin'; % a última camada usará função linear Split
end

yc=cell(nc,2); % As saídas Split de cada camada interior da NN
a=cell(nc,2); % Entrada das funções de activação de cada camada

wx=W{1}*x+b{1};
wx1=(W{1}.*alfW{1})*x+b{1}.*alfb{1};
wx2=wx-wx1;

I=wx>4;
if any(I,'all')
    wx1(I)=wx1(I)./wx(I).*4;
    wx2(I)=wx2(I)./wx(I).*4;
    wx(I)=4;
end

I=wx<-4;
if any(I,'all')
    wx1(I)=-wx1(I)./wx(I).*4;
    wx2(I)=-wx2(I)./wx(I).*4;
    wx(I)=-4;
end

% I=wx1<-4 | wx1>4;
% if any(I,'all')
%     I
% end
% I=wx2<-4 | wx2>4;
% if any(I,'all')
%     I
% end
a(1,:)={wx1,wx2};
yc(1,:)=feval(Fn{1},a(1,:)); % saida da camada 1
yy=feval(Fn{1},wx); % saida da camada 1
%mean(sum((yc{1,1}+yc{1,2}-yy).^2))

for i=2:nc
    wx=W{i}*yy+b{i};
    wx1=(W{i}.*alfW{i})*yc{i-1,1}+b{i}.*alfb{i};
    wx2=wx-wx1;

    I=wx>4;
    if any(I)
        wx1(I)=wx1(I)./wx(I).*4;
        wx2(I)=wx2(I)./wx(I).*4;
        wx(I)=4;
    end

    I=wx<-4;
    if any(I)
        wx1(I)=-wx1(I)./wx(I).*4;
        wx2(I)=-wx2(I)./wx(I).*4;
        wx(I)=-4;
    end
    a(i,:)={wx1,wx2};
    yc(i,:)=feval(Fn{i},a(i,:)); % saidas da camada i
    yy=feval(Fn{i},wx); % saida (total) da camada i
  
%    yc{i,2}=yy-yc{i,1};
end
y=yc(nc,:); % Saída da Rede (= saída da última camada)

end

% --------------------------------------
function [alfW,alfb,E]=splitNNbp(x,d,d1,niter,m,W,b,Fn,alfW,alfb)
% Algoritmo de SPLIT da rede neuronal.
% Os pesos SPLIT são treinados pelo algoritmo de BackPropagation (BP)
%
% Paulo Salgado @ UTAD2020

[nx,n]=size(x); % nx - dimensão do vector de entrada; n - n. de dados
[ny,n]=size(d);

nc=length(W); % n. de camadas

best.alfW=alfW;
best.alfb=alfb;
best.error=Inf;

% ------------------------------
% Processo iterativo
% ------------------------------.
delta=cell(nc,1);

E=zeros(niter,1); % Erro da rede em cada iteração de aprendizagem

ne=0;
flag=strcmp(Fn{nc},'purelin');
if flag
    nct=nc-1;
else
    nct=nc;
end

for ni=1:niter

    [y,yc,a]=my_net_split(x,W,b,Fn,alfW,alfb);

    % Calcular o Gradiente da função custo C face à saída dos neuronio da
    % última camada

    C=(y{1}-d1).^2/n;
    % dC_dy (Gadiente) derivada da função escalar C relativamente ao vector
    % de saída da rede, y
    
    E(ni)=sum(C,2);  % Erro na iteracção "ni"

    if E(ni)<best.error
        best.alfW=alfW;
        best.alfb=alfb;
        best.error=E(ni);
        ne=0;
    else
        ne=ne+1;
        if ne>10000
        for i=1:nc
            alfW{i}=rand(size(alfW{i})).*2-1;
            alfb{i}=rand(size(alfb{i})).*2-1;
        end
        ne=0;
        end
    end

% ====================
    % Retropropagar o Erro
    if flag
        
        % Método dos minimos quadrados
        A=[W{nc}'.*yc{nc-1,1};b{nc}.*ones(1,n)];
        wb=d1*A'*inv(A*A'+0.001*eye(size(A,1)));
        
        alfW{nc}=wb(:,1:end-1);
        alfb{nc}=wb(:,end);

        wx=W{nc}*(yc{nc-1,1}+yc{nc-1,2})+b{nc};
        wx1=(W{nc}.*alfW{nc})*yc{nc-1,1}+b{nc}.*alfb{nc};
        wx2=wx-wx1;
        yc(nc,:)={wx1,wx2}; % saidas da camada i
        a(nc,:)={wx1,wx2}; % saidas da camada i
        
        C=(wx1-d1).^2/n;
        Grad_C=2.*(wx1-d1)/n;
        
        delta{nc}=Grad_C;
    else
        %delta{nc}=((1-a{nc}).*(1+a{nc})).*Grad_C;
        Grad_C=2*(y{1}-d1)/n;
        
        dF=deriv_tansig_submodels(a{nc,1},a{nc,2});
        delta{nc}=dF.*Grad_C;
    end

    for i=nc-1:-1:1
        dF=deriv_tansig_submodels(a{i,1},a{i,2});
        delta{i}=dF.*((W{i+1}.*alfW{i+1})'*delta{i+1});
    end
    
    % -----------------
    % Ajustar os pesos
    % -----------------

    % camadas escondidas e de saída    
    for i=nct:-1:2
        miu=m*(0.5+0.5*(i-1)/(nc-1));
        for j=1:size(yc{i,1},1)
            %DW=sum(delta{i}(j,:).*a{i-1},2)';
            %W{i}(j,:)=W{i}(j,:)-miu.*DW;
            DalfW=sum(delta{i}(j,:).*(W{i}(j,:)'.*yc{i-1,1}),2)';
            alfW{i}(j,:)=alfW{i}(j,:)-miu.*DalfW;
        end
        alfb{i}=alfb{i}-miu*4*(sum(delta{i},2).*b{i});
    end
    
    % Camada de entrada
    for j=1:size(yc{1,1},1)
        %DW=sum(delta{1}(j,:).*x,2)';
        %W{1}(j,:)=W{1}(j,:)-m*0.5.*DW;
        DalfW=sum(delta{1}(j,:).*(W{1}(j,:)'.*x),2)';
        alfW{1}(j,:)=alfW{1}(j,:)-m*0.5.*DalfW;
    end
    alfb{1}=alfb{1}-m*4*(sum(delta{1},2).*b{1});
        
end

alfW=best.alfW;
alfb=best.alfb;
end

% --------------------------------------
function out=tansig(a)
% Função de activação Sigmoidal bipolar "tansig" ou "tansig SPLIT"
persistent coef;
if isempty(coef) % Coeficientes de Fourir da função TanSig
   coef=[1.19611055624426e+000;...
    252.971280799505e-003;...
    72.1236944285778e-003;...
    20.8874712293142e-003;...
    6.11703185161781e-003;...
    1.78229278188272e-003;...
    483.935721928844e-006;...
    202.802859460432e-006];
end

if ~iscell(a)
    out=2./(1+exp(-2*a))-1;
else
    % Função de activação Sigmoidal bipolar "tansig" decomposta em serie de
    % Fourier

    m=length(coef); % n. de harmonicos
   
    w0=pi/8;
    
    a1=a{1};
    a2=a{2};

    a=a1+a2;

    y1=coef(1).*(sin(w0.*a1).*cos(w0.*a2));
    y2=coef(1).*(cos(w0.*a1).*sin(w0.*a2));
    for i=2:m
        wi=(2*i-1).*w0;
        y1=y1+coef(i).*(sin(wi.*a1).*cos(wi.*a2));
        y2=y2+coef(i).*(cos(wi.*a1).*sin(wi.*a2));
    end
    out={y1 y2};
end

end

% --------------------------------------
function dF=deriv_tansig_submodels(a1,a2)
% Derivada do submodelo Fourier da função TanSig
persistent coef;
if isempty(coef) % Coeficientes de Fourir da função TanSig
   coef=[1.19611055624426e+000;...
    252.971280799505e-003;...
    72.1236944285778e-003;...
    20.8874712293142e-003;...
    6.11703185161781e-003;...
    1.78229278188272e-003;...
    483.935721928844e-006;...
    202.802859460432e-006];
end

m=length(coef); % n. de harmonicos
   
a=a1+a2;
I=a>4;
if any(I)
    a1(I)=a1(I)./a(I)*4;
    a2(I)=a2(I)./a(I)*4;
end
I=a<-4;
if any(I)
    a1(I)=-a1(I)./a(I)*4;
    a2(I)=-a2(I)./a(I)*4;
end

da=a1-a2;

w0=pi/8;
dF=0;
for i=1:m
    wi=(2*i-1).*w0;
    dF=dF+coef(i).*wi.*cos(wi*da);   
end

end
% --------------------------------------
function s=purelin(a)
% Função de activação linear (tipo y=x)
    s=a;
end
% --------------------------------------