function [ym,yl,ylm,delta] = sugval(p,x,f,def,S)
% SUGVAL calculates the output of the Sugeno-Takagi model.
%
% [Ym,Yl,Ylm] = SUGVAL(P,X,F)
%    Input:
%       P .....	consequents parameters for every cluster
%               (obtained for instance with SUGLMS).
%       X ..... input data matrix.
%       F ..... fuzzy partition matrix (membership degrees),
%		optional, defaults to ones(size(x,1),1) for
%		which SUGVAL is a standard linear model
%       DEF ... default value returned when the sum of grades
%               equals to one (optional, defaults to 0)
%	S   ... optional matrix generated by SUGLMS to produce error 
%		estimates on  predictions, see SUGLMS for details. 
%
%    Output:
%	Ym ....	global model output for the given input data
%	Yl ....	output of local submodels (corresponding to clusters)
%	Ylm ...	output of local submodels with data corresponding to
%		degrees < 0.2 masked (for plots)
%	delta . error estimates, y +/- delta contains at least 50% 
%		of the predictions.

% (c) Robert Babuska, January 1994

if nargin < 4, def = 0; end;              	% no default supplied
if isempty(def), def = 0; end;
if nargin < 3, f = ones(size(x,1),size(p,1)); end;          % no memberships supplied
if isempty(f), f = ones(size(x,1),size(p,1)); end;
[mx,nx] = size(x); [mf,nf] = size(f);
sumDOF = sum([zeros(1,mf);f'])';		% sum of DOFs
sumDOF = sumDOF(:,ones(1,nf));
NoRule = sumDOF == 0;                           % no rule applicable
sumDOF = sumDOF + NoRule;                       % set zeros to one
yl = x*p';                                      % local models
ym = sum([zeros(1,mx); ...
         (yl.*f./sumDOF)'])' + def.*NoRule(:,1); % global model
ylm = yl;
mask = find(f < 0.2);				% find membership degrees < 0.2
ylm(mask) = NaN*ones(size(mask));		% mask with NaN's for plots

if nargin > 4 & nargout > 3,
    nc = nf*nx;
    [ms,ns] = size(S);
    if (ms ~= ns+2) | (nc ~= ns)
        error('S matrix must be n+2-by-n where n = length(p(:))')
    end
    R = S(1:nc,1:nc);
    df = S(nc+1,1);
    normr = S(nc+2,1);

    x1 = zeros(mx,nf*nx); 	 			% auxilliary variables
    xx = zeros(nx,nf*mx);                   		% auxilliary variable
    x = x'; f1 = x(:); xx(:) = f1(:,ones(1,nf));
    f1 = x1;
    x1(:) = xx';                                        % reshape data matrix
    xx = f(:)./sumDOF(:);
    f1(:) = xx(:,ones(1,nx));                           % reshape partition matrix
    x1 = f1.*x1;

    E = x1/R;
    e = sqrt(1+sum((E.*E)')');
    if df == 0
        disp('Warning: zero degrees of freedom implies infinite error bounds.')
        delta = Inf*e;
    else
        delta = normr/sqrt(df)*e;
    end
end;

