% 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(100000,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=max(-1,min(1,1-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
m=8;  % n. de termos da Serie de Fourier
A=sin(wo.*x);
for i=2:m
    wi=(2*i-1)*wo;
    A=[A sin(wi.*x)./factorial(2*i-1)];
end
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');

figure(2);
coefF=(a./factorial((1:2:(2*m-1))'))'
ar=A.*coefF;
plot(x(I),ar(I,:));

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

% Derivada
x=(-4:0.01:4)';
n=length(x);
dTF=coefF(1)*wo*cos(wo.*x);
for i=2:m
    wi=(2*i-1)*wo;
    dTF=dTF+coefF(i)*wi*cos(wi*x);
end
dT=1-((exp(x)-exp(-x))./(exp(x)+exp(-x))).^2;

figure(3);
plot(x,dT,'b',x,dTF,'--r');

E_dT=mean((dT-dTF).^2)
save coef_Fourier_TanSig coefF;

return;
% --------------------------------------------------------
% Repartição da entrada, x, entre as componentes x1 e x2.
% --------------------------------------------------------
alfa=1.3:-0.1:0; % Factor de repartição
alfa=[2;1.5;1;0.8;0.6;0.55;0.5;0.45;0.4;0.2;0;-0.5;-1;-2];

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 k=1:n
    %Separação do sinal x em duas parcelas: x1 e x2.
    x1=x.*alfa(k);
    x2=x.*(1-alfa(k));
    
    % 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);
    A2=cos(wo.*x1).*sin(wo*x2);
    for i=2:m
        wi=(2*i-1)*wo;
        A1=[A1 sin(wi.*x1).*cos(wi.*x2)./factorial(2*i-1)];
        A2=[A2 cos(wi.*x1).*sin(wi.*x2)./factorial(2*i-1)];
    end

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

a1=a(1:m); % Parametros do modelo 1
a2=a(m+1:2*m);% 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=1:n
    subplot(2,7,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