source: trunk/toolbox/MLPfit.m @ 35

function  [w1,w2,errtot,rx] = MLPfit(Xv,Yv,Xa,Ya,w1,w2,F1,F2,OptimArg,Visu)

%MLPfit performs iterations of the MLP by the Harris method

%

%       [W1,W2,errtot,Yr] = MLPfit(Xv,Yv,Xa,Ya,W1,W2,F1,F2,OptimArg,Visu)

%

% Xv  input validation data

% Yv  output validation data 

% Xa the input learning data

% Ya the output leraning data

% W1 the initial parameter matrix from the input to the hidden layer

% W2 the initial parameter matrix from the hidden layer to the output

% F1 (default 'tah') the activation function of the hidden layer units.

%    Choose in {'tah', 'sig'}

% F2 (default 'lin') the activation function of the output layer units

%    Choose in {'tah', 'sig', 'lin'}

%

% OptimArg (default [1000, 0.000001]) Optimisation parameters: NbIter or

%         [NbIter, Threshold]

% Visu (default [ -, 1]) Visulalisation parameters: DisplayFrequency or

%          [DisplayFrequency, DisplaYdevice]

%

% W1 and W2 Final parameter matrix

% errtot the training error

% Yr the the output of the network at point x

%

% modified on Nov. 15th, 1999 (Evry, France)



if nargin < 6;

   help MLPfit

   error(sprintf('\n *** MLPfit error: invalid call***\n\n\t[W1,W2,errtot,Yr] =MLPfit(Xa,Ya,W1,W2,F1,F2,OptimArg,Visu);\n\n'));

end;



nbre_err=10 ;%nbre d'erreur de suites maxi
i_err=0; %compteur du nbre d'erreur

% Optimization control parameters



nbitemax = 1000; % default values

seuil = 10^(-6);



if nargin >= 9

        nbitemax = OptimArg(1);

        lO = length(OptimArg);

        if lO > 1

                seuil = OptimArg(2);

        end;

end;



%%plot(Xa,MLPval(Xa,w1,w2,F1,F2),'b-',X,Yb,'r-')

%plot(Xa,MLPval(Xa,w1,w2,F1,F2),'b-')



%pause(0.1)





% Optimisation Process monitoring visualisation parameters



df = nbitemax+2;                % Defaults is NO output

fid = 1;                            % default is standart output                

if nargin == 10

 df = Visu(1);                  % display frequency

 if df == 0; df =  nbitemax+2; end;

 lV = length(Visu);

 if lV > 1

        fid = Visu(2);      % Display Device

 end;

end;



% Check for F1 and F2



  if (strcmp(lower(F1),'tah') == 0) & (strcmp(lower(F1),'sig') == 0)

    F1 = 'tah'; disp('Unknown activation function 1. - Set as hyperbolictangent');

  end



  if (strcmp(lower(F2),'tah') == 0) & (strcmp(lower(F2),'sig') == ...
                                       0) & (strcmp(lower(F2),'lin') ...
                                             == 0) & (strcmp(lower(F2),'exp')==0)

    F2 = 'lin'; disp('Unknown output activation function - Set as linear');

  end



% initialisation



  ell=size(Xa,1);

  onell = ones(ell,1);

  [t nout] = size(w2);

  errold = 1e16;

  nbite = 0;

  w1p = w1;               % sauvegardes

  w2p = w2;

  w1pp = zeros(size(w1));                 % sauvegardes

  w2pp = zeros(size(w2));

  gradw1p = zeros(size(w1));

  gradw2p = zeros(size(w2));

  descw1 = gradw1p;

  descw2 = gradw2p;

  pas1 = 0.1 * ones(size(w1))*2/ell;

  pas2 = 0.1 * ones(size(w2))*2/ell;

  dim = 0.5;              % diminution du pas en cas d'augmentation de l'erreur

  a = 1.5;

  b = 1/a;

  alpha = .9/(1+.9);                    % momentum



 lambda = 0.*ones(size(w2));

 lambda(t,:) = zeros(1,nout);

  errtmin = 10000000000000000000000;



  ovf1 = 1e3 * ones(size(w1));

  ovf2 = 1e3 * ones(size(w2));



  errtot = [];                          % memoire erreur pour sortie



% boucle principale



continumongars=1;



while (nbite < nbitemax) & continumongars ;  %
  
  
  
  %Prop
  
  
  
  a1 = [Xa onell]*w1;
  
  
  
  if strcmp(lower(F1),'tah')
    
    x1 = tanh(a1);
    
  elseif strcmp(lower(F1),'sig')
    
    x1 = phi(a1);
    
  else
    
    error('Unknown activation function 1.')
    
  end
  
  
  
  a2 = [x1 onell]*w2;
  
  
  
  if strcmp(lower(F2),'tah')
    
    y = tanh(a2);
    
  elseif strcmp(lower(F2),'sig')
    
    y = phi(a2);
    
  elseif strcmp(lower(F2),'lin')
    
    y = a2;
    
    
  elseif strcmp(lower(F2),'exp')
    
    y = exp(a2);
  else
    
    error('Unknown output activation function')
    
  end
  
  
  
  err =  (y - Ya);
  
  errnew = sum(sum(err'.*err')) + sum(sum(lambda.*w2.^2));
  
  
  if ~isempty(Xv)
    Yvi=MLPval(Xv,w1,w2,F1,F2);
    err2=(Yv-Yvi);
    errval = sum(sum(err2'.*err2'));
    % keyboard
  else
    errval=0;
  end
  
  errtotnew=[errnew,errval];
  
  errtot = [errtot ;errtotnew];
  
  
  
  

  if errnew >= errold ;                 % ca merde => un pas en arriere
    
    
    i_err=i_err+1;
    
    if i_err>nbre_err
      disp('erreur dans l''apprentissage')
      rx = y;
      return
    end
    
    w1 = w1p ;                          
    
    w2 = w2p ;
    
    gradw1 = gradw1p ;
    
    gradw2 = gradw2p ;
    
    descw1 = gradw1;
    
    descw2 = gradw2;
    
    pas1 = dim * pas1 ;
    
    pas2 = dim * pas2 ;
    
    w1 = w1 - pas1 .* gradw1;
    
    w2 = w2 - pas2 .* gradw2;
    
    
    
  else                                  % ca marche => un pas en avant
    
    
    i_err=0;
    % BP
    
    
    
    dJdy = 2*err;
    
    
    
    if strcmp(lower(F2),'tah')
      
      dJda2 = dJdy.*(1-y.*y);
      
    elseif strcmp(lower(F2),'sig')
      
      dJda2 = dJdy.*(y-y.*y);
      
    elseif strcmp(lower(F2),'lin')
      
      dJda2 = dJdy;
      
    elseif strcmp(lower(F2),'exp')
      
      dJda2 = dJdy.*y;
      
    else

      error('Unknown output activation function')
      
    end
    
    
    
    gradw2 = (([x1 onell]'*dJda2) + (2*lambda.*w2))./ell;
    
    
    
    if strcmp(lower(F1),'tah')
      
      dJdx1 =  (w2(1:t-1,:) * dJda2')' .*(1-x1.*x1);
      
    elseif strcmp(lower(F1),'sig')
      
      dJdx1 =  (w2(1:t-1,:) * dJda2')' .*(x1-x1.*x1);

    else
      
      error('Unknown activation function 1.')
      
    end
    
    

    gradw1 = [Xa onell]' * dJdx1  ./ell;
    
    
    
    errold = errnew;                    
    
    w1p = w1 ;
    
    w2p = w2 ;
    
    gradw1p = gradw1 ;
    
    gradw2p = gradw2 ;
    
    descw1 = (1-alpha) * gradw1 + alpha * descw1;
    
    descw2 = (1-alpha) * gradw2 + alpha * descw2;
    
    
    
    test1 = (gradw1 .* descw1) >= 0;
    
    pas1 = ((test1 * a) + ((~test1) * b)) .* pas1;
    
    pas1 = (pas1 <= ovf1) .* pas1 + (pas1>ovf1) .* ovf1;
    
    
    
    test2 = (gradw2 .* descw2) >= 0;
    
    pas2 = ((test2 * a) + ((~test2) * b)) .* pas2;
    
    pas2 = (pas2 <= ovf2) .* pas2 + (pas2>ovf2) .* ovf2;
    
    
    
    w1 = w1 - pas1 .* descw1;
    
    w2 = w2 - pas2 .* descw2;
    
    
    
    nbite = nbite + 1;       % il n'y a que les bonnes iterations qui comptent
    
    
    
    if (rem(nbite,df) == 0),
      
      if (rem(nbite,20*df) == df),
        
        
        
        disp(['|  #epoch | Errorapp|Errorval| Mean Gw1 |  Max Gw1 | Mean Gw2 |  Max Gw2 |Max: ' num2str(nbitemax)])
        
        
        
        
        
      end

      
      
      fprintf(fid,'| %7.0f | %7.4f  | %7.4f | %6.6f | %6.6f | %6.6f | %6.6f|\n',[nbite errnew/size(Xa,1)  ...
                          errval/size(Xv,1),mean(mean(abs(gradw1))) max(max(abs(gradw1))) mean(mean(abs(gradw2))) max(max(abs(gradw2))) ])
      
      
      
      %figure(2)
        %hold on
        %%plotplot(Xa,MLPval(Xa,w1,w2,F1,F2),'b-',X,Yb,'r-')
        %plot(Xa,MLPval(Xa,w1,w2,F1,F2),'b-')
        %S=sprintf('en rouge fonction utilisée pour simuler les données en bleu courbe apprise par le MLP');
        %title(S);
      
      %pause(0.1)
      
      

      w1min = w1;
      
      w2min = w2;
      
      
      
      bougeti=max([max(max((abs(w1pp-w1p)./max(abs(w1p),1)))) max(max((abs(w2pp-w2p)./max(abs(w2p),1))))]);
      
      if bougeti< seuil
        
        continumongars=0;
        
      end;
      
      %        disp([nbite/100 errnew-sum(sum(lambda.*w2.^2)) mean(mean(pas1)) bougeti sum(mean(abs(errtest)))]);
      
      if (rem(nbite,1000) == 0),
        
        w1pp=w1p;w2pp=w2p;
        
        pas1 = mean(mean(pas1)) * ones(size(w1))*2/ell;
        
        pas2 = mean(mean(pas2)) * ones(size(w2))*2/ell;
        
      end;
      
    end;
    
    
    %disp(int2str(nbite))
    
    
  end % fin de l'adaptation des pas.
  
  
  
  
end;            % fin de la boucle principale



%w1min = w1p;w2min = w2p;



rx = y;