%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%Rede Neural com uma camada escondidada%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exemplo - Problema do Ou-Exclusivo
% X =[0 0;0 1;1 0;1 1];
% Y = [0;1;1;0]
% [Aval,Bval]=rede_neural(X,Y,X,Y,2,2000)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Aval,Bval]=rede_neural(X,Yd,Xval,Ydval,h,nepocasmax)

[N,ne]= size(X);
Nval = size(Xval,1);
X = [X,ones(N,1)];
Xval=[Xval,ones(Nval,1)];
ns= size(Yd,2);

A = rands(h,ne+1);
B = rands(ns,h+1);
Yr = calc_saida(X,A,B,h,N);
[dJdA,dJdB]=calc_grad(X,Yd,A,B,N);
grad = [dJdA(:);dJdB(:)];
erro = Yr-Yd;
EQM = 1/N*sum(sum(erro.*erro));

Yrval = calc_saida(Xval,A,B,h,Nval);
erroval = Yrval-Ydval;
EQMval = 1/N*sum(sum(erroval.*erroval));

vetEQM =[];
vetEQM = [vetEQM,EQM];
nep = 0;
%alfa = 0.1;
while norm(grad)>1e-5 & nep<nepocasmax
    nep = nep+1;
    alfa=calc_alfa(X,Yd,A,B,dJdA,dJdB,N)
    A = A - alfa*dJdA;
    B = B - alfa*dJdB;
    Yr = calc_saida(X,A,B,h,N);
    [dJdA,dJdB]=calc_grad(X,Yd,A,B,N);
    grad = [dJdA(:);dJdB(:)];
    erro = Yr-Yd;
    EQM = 1/N*sum(sum(erro.*erro));
    vetEQM = [vetEQM,EQM]; 
    Yrval = calc_saida(Xval,A,B,h,Nval);
    erroval = Yrval-Ydval;
    EQMval_new = 1/N*sum(sum(erroval.*erroval));
    if EQMval_new <EQMval
        EQMval = EQMval_new;
        Aval = A;
        Bval = B;
    end
    
end
plot(vetEQM)
vetEQM(end)
end

function Y = calc_saida(X,A,B,h,N)
% Zin = zeros(N,h);
% for n=1:N
%     for j=1:h
%         Zin(n,j)=A(j,:)*X(n,:)';
%     end
% end
Zin = X*A';
Z = 1./(1+exp(-Zin));
Zb = [Z,ones(N,1)];
Yin = Zb*B';
Y = 1./(1+exp(-Yin));
end    

function [dJdA,dJdB]=calc_grad(X,Yd,A,B,N)
Zin = X*A';
Z = 1./(1+exp(-Zin));
Zb = [Z,ones(N,1)];
Yin = Zb*B';
Y = 1./(1+exp(-Yin));
erro = Y-Yd;
fl = (1-Z).*Z;
gl = (1-Y).*Y;

dJdB = 1/N*(erro.*gl)'*Zb;
dJdZ = (erro.*gl)*B(:,1:end-1);
dJdA = 1/N*(dJdZ.*fl)'*X;


% dJddB = zeros(ns,h+1);
% for k=1:ns
%     for i=1:h+1
%         for n=1:N
%             dJdB(k,i)=  dJdB(k,i)+erro(n,k)*gl(n,k)*Z(n,i);
%         end
%     end
% end
% dJdZ = zeros(N,h);
% for n=1:N
%     for i=1:h
%         for k=1:ns
%             dJdZ(n,i) = dJdZ(n,i)+_erro(n,k)*gl(n,k)*B(k,i);
%         end
%     end
% end

end

function alfa_m=calc_alfa(X,Yd,A,B,dJdA,dJdB,N)

d = [-dJdA(:);-dJdB(:)];
alfa_l =0;
alfa_u =rand(1,1);
Aaux = A - alfa_u*dJdA;
Baux = B - alfa_u*dJdB;
[dJdAaux,dJdBaux]=calc_grad(X,Yd,Aaux,Baux,N);
grad = [dJdAaux(:);dJdBaux(:)];
hl = grad'*d;
while hl<0
    alfa_l = alfa_u;
    alfa_u = 2*alfa_u;
    Aaux = A - alfa_u*dJdA;
    Baux = B - alfa_u*dJdB;
    [dJdAaux,dJdBaux]=calc_grad(X,Yd,Aaux,Baux,N);
    grad = [dJdAaux(:);dJdBaux(:)];
    hl = grad'*d;
end
alfa_m = (alfa_l+alfa_u)/2;
Aaux = A - alfa_m*dJdA;
Baux = B - alfa_m*dJdB;
[dJdAaux,dJdBaux]=calc_grad(X,Yd,Aaux,Baux,N);
grad = [dJdAaux(:);dJdBaux(:)];
hl = grad'*d;

k =0;
kmax = ceil(log((alfa_u - alfa_l)/1e-4));

while k<kmax & abs(hl)>1e-4 
      
    if hl>0
        alfa_u = alfa_m;
    elseif hl<0
        alfa_l = alfa_m;
    else
        break;
    end
    alfa_m = (alfa_l+alfa_u)/2;
    Aaux = A - alfa_m*dJdA;
    Baux = B - alfa_m*dJdB;
    [dJdAaux,dJdBaux]=calc_grad(X,Yd,Aaux,Baux,N);
    grad = [dJdAaux(:);dJdBaux(:)];
    hl = grad'*d;
end

end