function [net,sigma] = mixmlpem(net, x, t, options, niters, alg)
%mixmlpem  EM algorithm for gated experts 
%
%       NET = MIXMLPEM(NET, X, T, OPTIONS, NITERS, ALG) uses the
%       Expectation Maximization (EM) algorithm to set the parameters of a 
%       gated experts model NET on a data set {X,T}. Each row of X 
%       corresponds to one input vector and each row of T corresponds to 
%       one target vector. The M-step of the EM algorithm is performed 
%       seperately for each of the expert networks and the gating network.
%       The M-step is performed using  standard non-linear optimisation 
%       routines (scg, gradient descent, quasi-Newton, conjugate gradient).
%
%       The optional parameters OPTIONS are passed on to these optimization 
%       routines and their complete interpretation can be found there. Here 
%       it is useful to know that OPTIONS(14) specifies the number of iterations
%       in the M-step. NITERS specifies the maximum number of outer
%       iterations (EM-steps).
%
%       If the gating network is a Gaussian mixture model (NET.TAG = 'kernel'|'ker&var')
%       the M-step for the gating network requires only a single pass. 
%
%       [NET,SIGMA] = MIXMLPEM(NET, X, T, OPTIONS, NITERS, ALG) also returns
%       the values of the expert variance parameters (default 1). 
%
%       See also
%       mixmlp, mixmlpfwd, mixpost, mixopt, kernel, scg, quasinew, graddesc, conjgrad
%
%       Copyright: Perry Moerland 


% Check that inputs are consistent
errstring = consist(net, 'mixmlp', x);
if ~isempty(errstring)
  error(errstring);
end

sigma = ones(net.nexp,1);

mse =(2*size(t,1));

for n=1:niters
  
 [y, y_exp, y_ehidden, y_gate] =  mixmlpfwd(net,x);
 
 % E-step: calculation of the posteriors
 switch net.tag
   case {'variance','ker&var','mfa&var','mpca&var'}
      post = mixpost(net,y_exp,y_gate,t,sigma); 
   case {'standard','kernel','mfa','mpca'}
      post = mixpost(net,y_exp,y_gate,t);
 end;
  

 E_me = sum(sum((y-t).^2))/mse;
 
 fprintf('iteration # %i   MSE = %f\r', n, E_me)  % Print on-line inform
 
 % M-step for the expert networks

 if ~(options(14))
  % the default number of M-steps
  options(14) = 3;
 end;
 
  switch net.tag
   case {'variance','ker&var','mfa&var','mpca&var'}
     for i=1:net.nexp
       n1 = t-y_exp(:,:,i);
       s = (sum(post(:,i)));
       sigma(i) =  (1/net.nout)*sum(post(:,i).*sum(n1.^2,2))/s;
       net.exp(i) = mixopt(net.exp(i), options, x, t, post(:,i),'mixmlperr_exp', 'mixmlpgrad_exp', alg, sigma(i));
     end;
   case {'standard','kernel','mfa','mpca'}
     for i=1:net.nexp
       net.exp(i) = mixopt(net.exp(i), options, x, t, post(:,i), 'mixmlperr_exp', 'mixmlpgrad_exp', alg);
     end;
  end;
 
 % M-step for the gating network
 switch net.tag
  case {'kernel','ker&var'}
    % gating network is a Gaussian mixture model
    net.gate = kernel(net.gate, x, post, options); 
  case {'mfa','mfa&var'}
    % gating network is a mixture of FAs
    net.gate = mfa_kernel(post,x,net.gate,1);
  case {'mpca','mpca&var'}
    % gating network is a mixture of PCAs
    net.gate = mpca_kernel(post,x,net.gate);
  otherwise
    net.gate = mixopt(net.gate, options, x, t, post, 'mixmlperr_gate','mixmlpgrad_gate', alg);
 end;
end;

% MSE at the end of training
[y] =  mixmlpfwd(net,x);
E_me = sum(sum((y-t).^2))/mse;

fprintf('iteration # %i   MSE on training set = %f\r\n', niters+1, E_me)  % Print on-line inform

% for i=1:net.nout
%  figure
%  subplot(211)
%  plot(t(:,i),'bx'); hold on
%  plot(y(:,i),'mo');hold off
%  if net.nout==1
%   title('Measurement data = x     Network output=o')
%  else
%   title(['Measurement data = x     Network output=o      (output # ',num2str(i) ')'])
% end;
%  grid;
% 
%  subplot(212)
%  plot(t(:,i)-y(:,i));
%  title('Prediction error')
%  grid
%  subplot(111)
%  drawnow
% end;


