function [e,y,a] = mixerr_exp(net, x, t, post, varargin)
%mixerr_exp Evaluate error function for one expert network in a mixture of experts model trained by EM
%
%       E = MIXERR_EXP(NET, X, T, POST) takes an expert network data
%       structure NET (which is a GLM) together with a matrix X of input 
%       vectors and a matrix T of target vectors, the posterior
%       probabilities POST, and evaluates the error function with respect 
%       to the network weights. The error function is the expert contribution
%       to the EM error fucntion.
%
%       E = MIXERR_EXP(NET, X, T, POST, SIGMA) if the 
%       conditional  density for the experts is an isotropic Gaussian (with 
%       one variance parameter for each expert), SIGMA contains these
%       variances. 
%
%       [E,Y,A] = MIXERR_EXP(NET, X, T, POST, VARARGIN) also returns a
%       matrix Y giving the outputs of the gating network and a matrix
%       A giving the summed inputs to each output unit, where each row
%       corresponds to one pattern. 
%
%       See also
%       mixerr_em, mix, mixem, mixopt, glmfwd
%
%       Copyright: Perry Moerland

% Check arguments for consistency
errstring = consist(net, 'glm', x, t);
if ~isempty(errstring);
  error(errstring);
end

[y,a] = glmfwd(net, x);

switch net.actfn

  case 'linear'         % Linear outputs

    if size(varargin)
     s = varargin{1};
     e = 0.5*sum(post.*(sum((y - t).^2,2)/s));
    else
     e = 0.5*sum(post.*sum((y - t).^2,2));
    end;
  

  case 'logistic'       % Logistic outputs

    % Ensure that log(1-y) is computable: need exp(a) > eps
    maxcut = -log(eps);
    % Ensure that log(y) is computable
    mincut = -log(1/realmin - 1);
    a = min(a, maxcut);
    a = max(a, mincut);
    y = 1./(1 + exp(-a));
    e = - sum(post.*sum(t.*log(y) + (1 - t).*log(1 - y),2));

  case 'softmax'        % Softmax outputs
  
    nout = size(a,2);
    % Ensure that sum(exp(a), 2) does not overflow
    maxcut = log(realmax) - log(nout);
    % Ensure that exp(a) > 0
    mincut = log(realmin);
    a = min(a, maxcut);
    a = max(a, mincut);
    temp = exp(a);
    y = temp./(sum(temp, 2)*ones(1,nout));
    % Ensure that log(y) is computable
    y(y<realmin) = realmin;
    e = - sum(post.*sum(t.*log(y),2));

end
