////////////////////////////////////////////////////////////////////////////////
//
// CARGA en Scilab: getf rn.sci;
//

////////////////////////////////////////////////////////////////////////////////
// IMPLEMENTACION

// crea una capa de proceso vectorial
function c = capa(tipo,Ni,No)
  if Ni<1 | No< 1 then error(tipo+": entradas no calzan con salidas"); end
  select tipo
  case 'lineal'
    c = tlist( ['capa' 'tipo' 'Ni' 'No' 'Rw'],...
               tipo, Ni, No, 1:Ni*No )
  case 'afin'
    c = tlist( ['capa' 'tipo' 'Ni' 'No' 'Rw'],...
               tipo, Ni, No, 1:(Ni+1)*No )
  case 'logit'
    if Ni<>No then error("logit: entradas no calzan con salidas"); end
    c = tlist( ['capa' 'tipo' 'Ni' 'No' 'Rw'],...
               tipo, Ni, No, 0 )
  case 'tanh'
    if Ni<>No then error("tanh: entradas no calzan con salidas"); end
    c = tlist( ['capa' 'tipo' 'Ni' 'No' 'Rw'],...
               tipo, Ni, No, 0 )
  case 'errQM'  // error cuadratico medio
    if Ni<>No then error("errQM: entradas no calzan con salidas"); end
    c = tlist( ['capa' 'tipo' 'Ni' 'No' 'Rw'],...
               tipo, Ni, No, 0 )
  case 'errEC'  // entropia cruzada
    if Ni<>No then error("errQM: entradas no calzan con salidas"); end
    c = tlist( ['capa' 'tipo' 'Ni' 'No' 'Rw'],...
               tipo, Ni, No, 0 )
  else
    error(tipo+": tipo desconocido")
  end

// calcula la salida de una capa de proceso vectorial en funcion de la entrada x
function  y = evalCapa(capa,x)
  global("red","N","D")
  select capa('tipo')
  case 'lineal'
    w = matrix(red('W')(capa('Rw')),capa('Ni'),-1)
    y = w'*x
  case 'afin'
    w = matrix(red('W')(capa('Rw')),capa('Ni')+1,-1)
    y = w(1:$-1,:)'*x + w($,:)'*ones(x(1,:))
  case 'logit'
    y = (1+exp(-x)).^(-1)
  case 'tanh'
    y = 2*(1+exp(-2*x)).^(-1) - 1
  case 'errQM'
    y = sum((x-D).^2,'r')/N/2
  case 'errEC'
    v = x + (D==0).*(1-2*x)
    y = -sum(log(v))/N
  end

// calcula el gradiente retropropagado de una capa en funcion del gradiente
// en la salida gy. Tambien calcula gradiente de pesos gW.
function gx = gradCapa(capa,x,y,gy)
  global("red","N","D")
  select capa('tipo')
  case 'lineal'
    w = matrix(red('W')(capa('Rw')),capa('Ni'),-1)
    red('gW')(capa('Rw')) = red('gW')(capa('Rw')) + matrix(x*gy',1,-1)
    gx = w*gy
  case 'afin'
    w = matrix(red('W')(capa('Rw')),capa('Ni')+1,-1)
    red('gW')(capa('Rw')) = red('gW')(capa('Rw')) + matrix([x;1]*gy',1,-1)
    gx = w(1:$-1,:)*gy
  case 'logit'
    gx = y.*(1-y).*gy
  case 'tanh'
    gx = (1-y.*y).*gy
  case 'errQM'
    gx = (x-D)/N
  case 'errEC'
    v = x - (D==0)
    gx = -v.^(-1)/N
  end

// calcula el gradiente para entradas x y salidas deseadas d
function rgrad(x,d)
  global("red","N","D")
  N = size(x,'c')
  if size(x,'r')<>red('capas')(1)('Ni')
    error("evalRed: input no calza")
  end
  if size(d,'r')<>red('capas')($)('No')
    error("evalRed: output no calza")
  end
  if size(x,'c')<>size(d,'c')
    error("evalRed: input no calza con output")
  end
  red('gW') = zeros(red('gW'))
  for k=1:size(x,'c')
    red('x')(1)(:) = x(:,k)
    for l=1:size(red('capas'))
      red('x')(l+1)(:) = evalCapa(red('capas')(l),red('x')(l))
    end
    D = d(:,k)
    gy = gradCapa(red('errf'),red('x')(l+1),0,0)
    for l=size(r('capas')):-1:1
      gy = gradCapa(red('capas')(l),red('x')(l),red('x')(l+1),gy)
    end
  end

////////////////////////////////////////////////////////////////////////////////
// INTERFAZ

// crea una red neuronal de varias capas
function r = red(varargin)
  r = tlist(['red' 'capas' 'errf' 'W' 'gW' 'x'])
  Ni = varargin(1)
  k = 2;
  capas = list()
  while k< size(varargin)
    t = varargin(k)
    No = varargin(k+1)
    capas(size(capas)+1) = capa(t,Ni,No)
    Ni = No
    k = k+2
  end
  t = varargin(k)
  if ~or(t==['errQM' 'errEC']) then  error("red: falta funcion de error"); end
  r('errf') = capa(t,Ni,No)
  w = []
  x = list()
  for k=1:size(capas)
    x(k) = zeros(1:capas(k)('Ni'))'
    capas(k)('Rw') = size(w,'c')+capas(k)('Rw')
    w = [w,zeros(capas(k)('Rw'))]
  end
  x(size(capas)+1) = zeros(1:r('errf')('Ni'))'
  r('capas') = capas
  r('W') = w
  r('gW') = w
  r('x') = x

// obtiene la salida y de la capa l de una red r con entradas x
// por omision l es la ultima capa
function y = reval(r,x,l)
  global("red")
  red = r
  [lhs,rhs] = argn(0)
  if rhs < 3 then l = size(red('capas')); end
  if size(x,'r')<>red('capas')(1)('Ni')
    error("evalRed: input no calza")
  end
  for k=1:l
    x = evalCapa(red('capas')(k),x)
  end
  y = x

// calcula el error de la red r para cada entrada x y salida deseada d
function e = rerr(r,x,d)
  global("red","N","D")
  red = r
  N = size(x,'c'); D = d
  if size(x,'r')<>red('capas')(1)('Ni')
    error("evalRed: input no calza")
  end
  if size(d,'r')<>red('capas')($)('No')
    error("evalRed: output no calza")
  end
  if size(x,'c')<>size(d,'c')
    error("evalRed: input no calza con output")
  end
  for c=red('capas')
    x = evalCapa(c,x)
  end
  e = evalCapa(red('errf'),x)

// calcula el gradiente de la red r para entradas x y salidas deseadas d
function gw = grad(r,x,d)
  global("red")
  red = r
  rgrad(x,d)
  gw = red('gW')

// calcula el gradiente de la red r para entradas x y salidas deseadas d
// con diferencias finitas (se utiliza para verificar que la retropropagacion
// este bien implementada)
function gw = gradf(r,x,d)
  delta = .0000001
  w = r('W')
  gw = zeros(w)
  e0 = sum(rerr(r,x,d))
  for i=1:size(w,'c')
    v = w(i); w(i) = v+delta; r('W') = w; e = sum(rerr(r,x,d))
    gw(i) = (e-e0)/delta
    w(i) = v
  end

// entrenamiento batch de red r con entradas x i salidas deseadas d
// f es el factor de entrenamiento, c la cantidad de ciclos de l epocas
// una epoca es una pasada sobre todos los ejemplos
// la variable global ERR contiene la secuencia de errores para cada ciclo
function r = batch(r,x,d,f,c,l)
  global("red","ERR")
  [lhs,rhs] = argn(0)
  if rhs < 6 then l = 1; end
  red = r
  ERR = zeros(1,c+1)
  ERR(1) = sum(rerr(red,x,d))
  printf("%4i %6e\n",0,ERR(1))
  for k=1:c
    for e=1:l
      rgrad(x,d)
      red('W') = red('W') - f*red('gW')
    end
    ERR(k+1) = sum(rerr(red,x,d))
    printf("%4i %6e\n",k,ERR(k+1))
  end
  r = red

// crea una red neuronal y ejemplos para un problema de codificacion
// para probar: 
//   rtest();  // se crea e inicializa la red neuronal
//   r = batch(r,x,d,1,100,20);  // en total son 2000 epocas de entrenamiento
//   reval(r,x) // normalmente el resultado es cercano a una matriz identidad
//   reval(r,x,2)  // esta es la codificacion
//   global("ERR")  //  para tener acceso a los errores
//   plot(ERR)  //  muestra la curva de error
function rtest()
  r = red(4,'afin',2,'tanh',2,'afin',4,'logit',4,'errEC')  // creacion
  r('W') = (rand(r('W'))-.5)*2/10  // inicializacion aleatoria
  x = eye(4,4)  // entradas
  d = x  // salidas deseadas
  [r,x,d] = return(r,x,d)