// training set 
function tset=trainingSet()
   tset=tlist(["Training set","L","X","d"],0,zeros(1,100),zeros(1,100));
   tset.X=[0:1:10*%pi];
   tset.L=length(tset.X);
   tset.d=sin(tset.X)+0.2*rand(tset.X)-0.1*ones(tset.X);
endfunction  

// network creation
function network=RBFcreate(N,K,M,sigma)
   network=tlist(["RBF NN","N","K","M","X0","Y0","X1","Y1","X2","Y2","C","W","f","sigma"],N,K,M,zeros(N,1),zeros(N,1),zeros(K,1),zeros(K+1,1),zeros(M,1),zeros(M,1),zeros(K,N),zeros(M,K+1),gaussian,sigma);
endfunction

// gaussian activation function
function output=gaussian(x)
   output=exp(-x^2/2);
endfunction

// euclidean distance
function d=dist(x,y)
   d=sqrt((x-y)'*(x-y));
endfunction

// output signal computation
function network=forward(network,X)
 network.X0=X; network.Y0=network.X0;
 for k=1:network.K 
   network.X1(k)=dist(network.C(k)',network.Y0);
 end  
 network.Y1=[-1; feval(network.X1./network.sigma,network.f)]; 
 network.X2=network.W*network.Y1;
 network.Y2=network.X2;
endfunction

// error computation
function err=error_computation(network, tset)
err=0;
[N,L]=size(tset.X);
for i=1:L
     network=forward(network,tset.X(:,i));
     delta=(tset.d(i)-network.Y2).^2;
     err=err+sum(delta);
end
err=err/L
endfunction


// network training
function [network,ap,aE]=train(network,tset,pmax,Emax,eta)
L=tset.L;
network.C=tset.X';
network.W=2*rand(network.W)-ones(network.W);
E=error_computation(network,tset);
p=0;
ap=[];aE=[]; 
while p<pmax & E>Emax
 
  E=error_computation(network,tset);
  for i=1:L
     network=forward(network,tset.X(:,i)); 
     y=network.W*network.Y1;
     delta=tset.d(i)-y;
     network.W=network.W+eta*(delta*network.Y1')
  end;
  E=error_computation(network,tset);
  if (modulo(p,10)==0) then
      disp(p,"Iteration:");
      disp(E,"error=");
      ap=[ap p]; aE=[aE E]; 
  end
  p=p+1; 
 end 
 disp(E,"error=");
endfunction

// network testing
function output=test(network,X)
  network=forward(network,X);
  output=network.Y2;
endfunction  

function test_plot(tset,network)
  x=0:0.1:10*%pi;
  for i=1:length(x)
    network=forward(network,x(i));
    y(i)=network.Y2;
  end
  plot(tset.X,tset.d,'r*',x,y,'b');
endfunction


// regresie neliniara folosind o retea RBF
function tset=regressionRBF(eta,epochs,sigma)
tset=trainingSet();           // construirea setului de antrenare
rbf=RBFcreate(1,length(tset.X),1,sigma); // crearea retelei
disp(tset,"tset=");
[rbfa,ap,aE]=train(rbf,tset,epochs,0.0001,eta); // antrenarea retelei
// testarea retelei
test_plot(tset,rbfa);
pause;    // pauza in executie;  se continua tastand <resume>
clf;
plot(ap,aE);
endfunction

// apel functie
regressionRBF(0.1,1000,1)
