% Teste da decomposição (split) de um neurónio
% Nota: A decomposição faz-se pela divisão da função dos sinais de entrada
% x, em dois sub-sinais, x1 e x2, com a correspondente separação da função
% de activação.

Ts=16; % Periodo (?) da função tangente hiperbólica
wo=2*pi/Ts;
L=2;

W=randn(5,1); % Peso do neurónio
xi=randn(1000,5); % Entradas dos neurónios
x=xi*W; % Entrada da função de activação
x=max(-4,x);
x=min(4,x);

% Função de activação Tangente Hiperbólica
T=(exp(x)-exp(-x))./(exp(x)+exp(-x));
T2=(exp(L*x)-exp(-L*x))./(exp(L*x)+exp(-L*x));

figure(1);
[~,I]=sort(x);
plot(x(I),T(I),'b',x(I),tansig(x(I)),'r--');

% Decomposição da função Tangente Hiperbólica em serie de Fourier
A=[sin(wo.*x) sin(3*wo.*x)./6 sin(5*wo.*x)./120, ...
   sin(7*wo.*x)./5040 sin(9*wo.*x)./362880];
a=inv(A'*A)*A'*T;
Tf=A*a;

figure(1);
[~,I]=sort(x);
plot(x(I),T(I),'b',x(I),Tf(I),'r');

E=mean((T-Tf).^2)

% --------------------------------------------------------
% Repartição da entrada, x, entre as componentes x1 e x2.
% --------------------------------------------------------
alfa=1:-0.05:0; % Factor de repartição

n=length(alfa); % N. de diferentes valores de repartição

% Limitar dominio da função de activação ao intervalo [-4,4]

Ta=zeros(1000,n); % Regista os valores da partição 1
T1=zeros(1000,n); % Regista os valores da partição 1
T2=zeros(1000,n); % Regista os valores da partição 2
A=[];
for i=1:n
    %Separação do sinal x em duas parcelas: x1 e x2.
    x1=x.*alfa(i);
    x2=x.*(1-alfa(i));
    
    % Decomposição da função de activação (tangente hiperbólica) em duas
    % sub-funções (expressa em serie de Fourier)
    
    A1=[sin(wo.*x1).*cos(wo.*x2),...
       sin(3*wo.*x1).*cos(3*wo.*x2)./6,...
       sin(5*wo.*x1).*cos(5*wo.*x2)./120,...
       sin(7*wo.*x1).*cos(7*wo.*x2)./5040,...
       sin(9*wo.*x1).*cos(9*wo.*x2)./362880];
    A2=[cos(wo.*x1).*sin(wo.*x2),...
       cos(3*wo.*x1).*sin(3*wo.*x2)./6,...
       cos(5*wo.*x1).*sin(5*wo.*x2)./120,...
       cos(7*wo.*x1).*sin(7*wo.*x2)./5040,...
       cos(9*wo.*x1).*sin(9*wo.*x2)./362880];
   A=[A;[A1 A2]];
end

% Pesos das series de Fourier dos 2 sub-funções da função de activação
a=inv(A'*A+0.00*eye(10))*A'*repmat(T,n,1); % Parametros da decomposição de Fourier

a1=a(1:5); % Parametros do modelo 1
a2=a(6:10);% Parametros do modelo 2
[a1 a2]

Ta=A*a;
Ta=reshape(Ta,1000,n);

% Erro da aproximação
E=mean((T-Ta).^2)
figure(2);
plot(E);
title('Approximation error of TanSin function');

% Resposta dos sub-modelos
T1=A(:,1:5)*a1; % Resposta do sub-modelo 1
T2=A(:,6:10)*a2; % Resposta do sub-modelo 1
T1=reshape(T1,1000,n);
T2=reshape(T2,1000,n);

figure(3)
[~,I]=sort(x);
plot(x(I),T(I),'b',x(I),Ta(I,:),'r');
title('Approximation of TanSin function');

figure(4);
clf;
k=1;
for i=3:2:n
    subplot(2,5,k);
    plot(x(I),T(I),'k',x(I),T1(I,i),'b',x(I),T2(I,i),'m',x(I),T1(I,i)+T2(I,i),'r');
    title(['alfa=',num2str(alfa(i))]);
    k=k+1;
end

% Alfa entre 0.4 e 0.6 com mais detalhe

alfa=0.6:-0.05:0.4; % Factor de repartição

n=length(alfa); % N. de diferentes valores de repartição

% Limitar dominio da função de activação ao intervalo [-4,4]

Ta=zeros(1000,n); % Regista os valores da partição 1
T1=zeros(1000,n); % Regista os valores da partição 1
T2=zeros(1000,n); % Regista os valores da partição 2
A=[];
for i=1:n
    %Separação do sinal x em duas parcelas: x1 e x2.
    x1=x.*alfa(i);
    x2=x.*(1-alfa(i));
    
    % Decomposição da função de activação (tangente hiperbólica) em duas
    % sub-funções (expressa em serie de Fourier)
    
    A1=[sin(wo.*x1).*cos(wo.*x2),...
       sin(3*wo.*x1).*cos(3*wo.*x2)./6,...
       sin(5*wo.*x1).*cos(5*wo.*x2)./120,...
       sin(7*wo.*x1).*cos(7*wo.*x2)./5040,...
       sin(9*wo.*x1).*cos(9*wo.*x2)./362880];
    A2=[cos(wo.*x1).*sin(wo.*x2),...
       cos(3*wo.*x1).*sin(3*wo.*x2)./6,...
       cos(5*wo.*x1).*sin(5*wo.*x2)./120,...
       cos(7*wo.*x1).*sin(7*wo.*x2)./5040,...
       cos(9*wo.*x1).*sin(9*wo.*x2)./362880];
   A=[A;[A1 A2]];
end
