function [A,B]=rna(X,Yd,h)
[N,ne]=size(X);
ns=size(Yd,2);
A=rands(h,ne+1)/2; % [-1,1]-->[-0.2,0.2]
B=rands(ns,h+1)/2; % [-1,1]-->[-0.2,0.2]
Y=calc_saida(X,A,B,N);
erro = Y-Yd;
MSE = 1/N*sum(sum(erro.*erro));
ep = 0;
epmax = 20000;
E = 1e-3;
alfa =0.95;
vet=[];
vet=[vet;MSE];
while MSE>E & ep<epmax,
    ep=ep+1;
    [dEdAatual,dEdBatual]=calc_grad(X,Yd,A,B,N);

    if ep==1 | ep==numpar
        dEdA=dEdAatual;
        dEdB=dEdBatual;
        dEdAant=dEdAatual;
        dEdBant=dEdBatual;
    else
        dEdA=dEdAatual+beta*dEdAant;
        dEdB=dEdBatual+beta*dEdBaant;
        dEdAant=dEdAatual;
        dEdBant=dEdBatual;
    end
        
    A= A - alfa*dEdA;
    B= B - alfa*dEdB;
    Y=calc_saida(X,A,B,N);
    erro = Y-Yd;
    MSE = 1/N*sum(sum(erro.*erro));
    vet=[vet;MSE];
    
end
plot(vet)
end

function [dEdA,dEdB]=calc_grad(X,Yd,A,B,N)
X=[X,ones(N,1)];
Zin = X*A';
Z = 1./(1+exp(-Zin));

Zb=[Z,ones(N,1)];
Yin = Zb*B';
Y = 1./(1+exp(-Yin));

erro = Y-Yd;
dg = (1-Y).*Y;
df = (1-Z).*Z;

dEdB = 1/N*(erro.*dg)'*Zb;

dEdZ = (erro.*dg)*B(:,1:end-1);

dEdA = 1/N*(dEdZ.*df)'*X;
end

function Y=calc_saida(X,A,B,N)
X=[X,ones(N,1)];
Zin = X*A';
Z = 1./(1+exp(-Zin));

Z=[Z,ones(N,1)];
Yin = Z*B';
Y = 1./(1+exp(-Yin));
end