% ------------------------------------------------------------
% 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=200000;

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;
%%
% -----------------------
% IMPRIMIR FIGURAS (jpg)
% -----------------------
figure(1);
surf(X,Y,Z);

figure(2);
surf(X,Y,Z1);

figure(3);
surf(X,Y,Z2);

figure(4);
surf(X,Y,reshape(y,nx,ny));

figure(5);
surf(X,Y,Z-reshape(y,nx,ny));
title('Error');
grid on;

for i=1:5
    figure(i);
    h(1)=xlabel('\it{x}');
    h(2)=ylabel('\it{y}');
    set(h,'FontSize',14,'FontWeight','bold','FontName','Times New Roman');
end

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


str=['.\figures\',projectname,'_'];
for i=1:6
    strfig=[str,'fig',num2str(i),'.jpg'];
    print(i, '-djpeg', strfig);
end

else
    load dados_split1;
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})); % Vector (linha) da repartição dos pesos dos neurónios da camada i
   alfb{i}=ones(size(b{i})); % Vector (linha) da repartição dos pesos dos neurónios da camada i
   alfT{i}=ones(size(b{i})); % Vector (linha) da repartição dos pesos da serie de Fourier da Função de activação da camada i
end

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

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

figure(7);
surf(X,Y,reshape(ys{1},nx,ny));

figure(8);
surf(X,Y,reshape(ys{2},nx,ny));

figure(9);
surf(X,Y,reshape(ys{1}+ys{2},nx,ny));

for i=7:9
    figure(i);
    h(1)=xlabel('\it{x}');
    h(2)=ylabel('\it{y}');
    set(h,'FontSize',14,'FontWeight','bold','FontName','Times New Roman');
end

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

figure(10);
plot((1:niter)',E);
xlabel('Iteration');
ylabel('MSE');

figure(11);
for i=1:nc
    subplot(2,nc,i);
    sW=sum(abs(W{i}),2);
    bar(sqrt(sum(abs(alfW{i}.*W{i}),2)./sW));
    title(['layer',num2str(i)]);
    subplot(2,nc,i+nc);
    bar(sqrt(sum(abs(alfW{i}.*b{i}),2)));
end
legend(num2str(1:nc));

figure(12);
subplot(2,3,1);
surf(X,Y,Z1);
subplot(2,3,2);
surf(X,Y,Z2);
subplot(2,3,3);
surf(X,Y,Z);
subplot(2,3,4);
surf(X,Y,reshape(ys{1},nx,ny));
subplot(2,3,5);
surf(X,Y,reshape(ys{2},nx,ny));
subplot(2,3,6);
surf(X,Y,reshape(ys{1}+ys{2},nx,ny));

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

%save dados_split1_final;
% ------------------------------------------------------------
function [y,a]=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

a=cell(nc,1); % Célula com as saídas de cada camada interior da NN
a{1}=feval(Fn{1},W{1}*x+b{1}); % saida da camada 1
for i=2:nc
    a{i}=feval(Fn{i},W{i}*a{i-1}+b{i}); % saida da camada i
end
y=a{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

if nargin<7 % 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
% ------------------------------
% Processo iterativo
% ------------------------------
delta=cell(nc,1);

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

    [y,a]=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
    Grad_C=(y-d)/n;
    
    E(ni)=sum(C,2);  % Erro na iteracção "ni"

    % Retropropagar o Erro
    if strcmp(Fn{nc},'purelin')
        delta{nc}=Grad_C;    
    else
        delta{nc}=((1-a{nc}).*(1+a{nc})).*Grad_C;
    end

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

    % Camada de entrada
    for j=1:size(a{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:nc
        miu=m*(0.5+1.5*(i-1)/(nc-1));
        for j=1:size(a{i},1)
            DW=sum(delta{i}(j,:).*a{i-1},2)';
            W{i}(j,:)=W{i}(j,:)-miu.*DW;
        end
        b{i}=b{i}-miu*sum(delta{i},2);
    end
end

end

% ------------------------------------------------------------
% FUNÇÕES SPLIT
% ------------------------------------------------------------
function [y,a]=my_net_split(x,W,b,Fn,alfW,alfb,alfT)
% 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

a=cell(nc,2); % Célula com as saídas Split de cada camada interior da NN

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

a(1,:)=feval(Fn{1},{wx1,wx2}); % saida da camada 1
aa=feval(Fn{1},W{1}*x+b{1},alfT{1}); % saida da camada 1
%mean(sum((a{1,1}+a{1,2}-aa).^2))
for i=2:nc
    wx1=(W{i}.*alfW{i})*a{i-1,1}+b{i}.*alfb{i};
    %wx2=W{i}*aa+b{i}-(W{i}.*alfW{i})*a{i-1,1}-b{i}.*alfb{i};
    %aa=a{i-1,1}+a{i-1,2};
    wx2=W{i}*aa+b{i}-wx1;
    a(i,:)=feval(Fn{i},{wx1,wx2},alfT{i}); % saidas da camada i
    aa=feval(Fn{i},W{i}*aa+b{i}); % saida (total) da camada i

%    mean(sum((a{i,1}+a{i,2}-aa).^2))
end
y=a(nc,:); % Saída da Rede (= saída da última camada)

end

% --------------------------------------
function [W,b,alfW,alfb,alfT,E,Fn]=splitNNbp(x,d,d1,niter,m,W,b,Fn,alfW,alfb,alfT)
% 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.alfT=alfT;
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=true;

for ni=1:niter

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

    % 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
    Grad_C=2.*(y{1}-d1)/n;
    
    E(ni)=sum(C,2);  % Erro na iteracção "ni"

    if E(ni)<best.error
        best.alfW=alfW;
        best.alfb=alfb;
        best.alfT=alfT;
        best.error=E(ni);
%         if flag
%             fprintf('\n %d ->',ni);
%             flag=false;
%         end
%         %m=m*1.01;
%     else
%         ne=ne+1;
%         if ~flag
%             fprintf(' %d | ',ni);
%             flag=true;
%         end
%         if ne>niter*0.1
%         for i=1:nc
%             alfW{i}=rand(size(alfW{i})).*2-1;
%             alfb{i}=rand(size(alfb{i})).*2-1;
%         end
%         ne=0;
%         fprintf('*');
%         end
    end
% ====================
    nct=max(1,nc-min(floor(ni/2000),nc)); % numero de camada a treinar
% ====================
    % Retropropagar o Erro
    if strcmp(Fn{nc},'purelin')
        delta{nc}=Grad_C;
    else
        %delta{nc}=((1-a{nc}).*(1+a{nc})).*Grad_C;
        dF=deriv_tansig_submodels(y{1}-y{2},alfT{nc});
        delta{nc}=dF.*Grad_C;
    end

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

    % camadas escondidas e de saída
    for i=nc:-1:max(nct,2)
        miu=m*(0.5+0.5*(i-1)/(nc-1));
        for j=1:size(a{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,:)'.*a{i-1,1}),2)';
            alfW{i}(j,:)=alfW{i}(j,:)-miu.*DalfW;
        end
        alfb{i}=alfb{i}-miu*(sum(delta{i},2).*b{i});
        dFc=derivCoefs_tansig_submodels(a(i,:),alfT{i});
        alfT{i}=alfT{i}-miu*sum(dFc,2);
    end
    
    if nct==1
    % Camada de entrada
    for j=1:size(a{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,:)-miu*0.5.*DalfW;
    end
    alfb{1}=alfb{1}-m*(sum(delta{1},2).*b{1});
    dFc=derivCoefs_tansig_submodels(a(i,:),alfT{i});
    alfT{i}=alfT{i}-miu*sum(dFc,2);
    end

end

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

% --------------------------------------
function out=tansig(a,alfCoef)
% 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;
    r=0;
    s=0;
    if nargin==1
        for i=1:m
            wi=(2*i-1).*w0;
            r=r+coef(i).*sin(wi*a{1}).*cos(wi*a{2});
            s=s+coef(i).*cos(wi*a{1}).*sin(wi*a{2});
        end
    else
        for i=1:m
            wi=(2*i-1).*w0;
            r=r+coef(i).*alfCoef(i).*sin(wi*a{1}).*cos(wi*a{2});
            s=s+coef(i).*(1-alfCoef(i)).*cos(wi*a{1}).*sin(wi*a{2});
        end
    end
    out={r s};
end

end

% --------------------------------------
function dF=deriv_tansig_submodels(da,alfCoef)
% 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
   
w0=pi/8;
dF=0;

if nargin==1
    for i=1:m
        wi=(2*i-1).*w0;
        dF=dF+coef(i).*wi.*cos(wi*da);
    end
else
    for i=1:m
        wi=(2*i-1).*w0;
        dF=dF+coef(i).*alfCoef(i).*wi.*cos(wi*da);   
    end
end

end

% --------------------------------------
function dFc=derivCoefs_tansig_submodels(a,alfCoef)
% Derivada dos coeficientes 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
   
w0=pi/8;
n=size(a{1},1);
dFc=zeros(n,m);

for i=1:m
    wi=(2*i-1).*w0;
    dFc(:,i)=dFc(:,i)+coef(i).*alfCoef(i).*sin(wi*a{1}).*cos(wi*a{2});   
end
    
end

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