function [net,sigma] = mixem(net, x, t, options, niters, alg)
%mixem  EM algorithm for mixture of experts model
%
%       NET = MIXEM(NET, X, T, OPTIONS, NITERS, ALG) uses the
%       Expectation Maximization (EM) algorithm to set the parameters of a 
%       mixtures of 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.
%       There are various alternatives for the M-step specified by ALG:
%         'irls'      = uses the iteratively reweighted least-squares 
%                       algorithm for GLMs in the M-step.
%         'heuristic' = uses the formulation as weighted least-squares
%                       problems for a single pass M-step.
%          otherwise  = standard non-linear optimisation routines (scg,
%                       gradient descent, quasi-Newton, conjugate gradient)
%                       from the Netlab toolbox.
%
%       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] = MIXEM(NET, X, T, OPTIONS, NITERS, ALG) also returns
%       the values of the expert variance parameters (default 1). 
%
%       See also
%       mix, mixfwd, mixpost, mixtrain, mixopt, mixheuristic, glmtrain, glmheuristic, kernel, scg, quasinew, graddesc, conjgrad
%
%Copyright: Perry Moerland 

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

% initialize the expert variance
sigma = ones(net.nexp,1);

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

for n=1:niters
  
 [y, y_exp, y_gate] =  mixfwd(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 alg
 case 'irls'
  for i=1:net.nexp  % IRLS training 
   net.exp(i) = mixtrain(net.exp(i),options,x,t,post(:,i));
  end;
 
 case 'heuristic'   % Least-squares heuristic 
  for i=1:net.nexp
   net.exp(i) = mixheuristic(net.exp(i),x,t,post(:,i));
  end;
   
 otherwise          % non-linear optimization with 'alg'
  switch net.tag
    case {'variance','ker&var','mfa&var','mpca&var'}
     % with a variance parameter for each expert
     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),'mixerr_exp', 'mixgrad_exp', alg, sigma(i));
     end;
   case {'standard','kernel','mfa','mpca'}
     % standard
     for i=1:net.nexp
       net.exp(i) = mixopt(net.exp(i), options, x, t, post(:,i), 'mixerr_exp', 'mixgrad_exp', alg);
     end;
  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
   switch alg
    case 'irls'        % IRLS training
     net.gate = glmtrain(net.gate,options,x,post);
    case 'heuristic'   % Least-squares heuristic
     net.gate = glmheuristic(net.gate,x,post);
    otherwise          % non-linear optimization with 'alg'
     net.gate = mixopt(net.gate, options, x, t, post, 'mixerr_gate', 'mixgrad_gate', alg);
   end;
 end;
end;

% MSE at the end of training
[y, post] =  mixfwd(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 information

for i=1:net.nout
 figure
 subplot(211)
 % plotting the target and predicted outputs
 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;
% plotting the prediction error
 subplot(212)
 plot(t(:,i)-y(:,i));
 title('Prediction error')
 grid
 subplot(111)
 drawnow
end;
