function y=mgrade(x,mfs,domain)
% MGRADE(X,MFS,DOMAIN) returns the membership grade of elements of row
% vector X to fuzzy sets with membership functions defined in matrix MFS.
% MFS is a matrix containing membership functions definition data.
% Its size is M x 5, where M is the number of membership functions.
% Number in the first column determines the membership function type:
%
%       0 ... one everywhere,        1 ... Trapezoidal
%       2 ... Exponential,           3 ... Sigmoidal
%
% The remaining four columns define the membership function position
% and shape. If we denote them [p1 p2 p3 p4] then membership degree Ax
% of numbers x is as follows:
%
%       Ax = 0                  for x <= p1 or  x >= p2
%            1                  for x >= p2 and x <= p3
%            in (0,1)           for x >  p1 and x <  p2 or
%                                   x >  p3 and x <  p4
%
% The shape of the membership function curve for x between p1 and p2
% and between p3 and p4 is determined by membership function type,
% i.e. trapezoidal, exponential and sigmoidal. For curve equations
% see the User's Guide.

% Copyright (c) Robert Babuska, 1992-94

[k,l]=size(mfs);
x = x(:)';
lx = length(x);
mx = ones(1,lx);
if nargin < 3, domain = []; end;

if length(domain) == 0,

  if ~all(mfs(:,2)<=mfs(:,3) & mfs(:,3)<=mfs(:,4) & mfs(:,4)<=mfs(:,5)),
    error('Illegal membership function parameters');
  end;
  xx = ones(k,1)*x;
  y = 1.0*(xx >= mfs(:,3)*mx & xx <= mfs(:,4)*mx);
  for j=1:k
    if mfs(j,1) == 0,                    % ANY
      y(j,:) = mx;
    elseif mfs(j,1) == 1,                % trapezoidal
      if mfs(j,3) > mfs(j,2),  y(j,:) = y(j,:) + ...
	(x > mfs(j,2) & x < mfs(j,3)).*(x - mfs(j,2))/(mfs(j,3) - mfs(j,2));
      end;
      if mfs(j,5) > mfs(j,4),  y(j,:) = y(j,:) + ...
	(x > mfs(j,4) & x < mfs(j,5)).*(mfs(j,5) - x)/(mfs(j,5) - mfs(j,4));
      end;
    elseif mfs(j,1) == 2,                % exponential
      wl = mfs(j,3) - mfs(j,2);
      wr = mfs(j,5) - mfs(j,4);
      if wl > 0
	    y(j,:)=y(j,:)+(x<mfs(j,3)).*exp(-7*(x-mfs(j,3)).^2/wl^2);
      end;
      if wr > 0
	    y(j,:)=y(j,:)+(x>mfs(j,4)).*exp(-7*(x-mfs(j,4)).^2/wr^2);
      end;
    elseif mfs(j,1)==3,                  % sigmoidal
      wl = (mfs(j,3) - mfs(j,2))/2;
      wr = (mfs(j,5) - mfs(j,4))/2;
      if wl > 0
	y(j,:)=y(j,:)+0.5*((x>mfs(j,2)) & (x<=mfs(j,2)+wl)).*...
		      (((x-mfs(j,2))/wl).^2)+...
		      ((x>mfs(j,2)+wl) & (x<mfs(j,3))).*...
		      (1-0.5*(((mfs(j,3)-x)/wl).^2));
      end;
      if  wr > 0
	y(j,:)=y(j,:)+((x>mfs(j,4)) & (x<=mfs(j,4)+wr)).*...
		      (1-0.5*(((mfs(j,4)-x)/wr).^2))+...
		      0.5*((x>mfs(j,4)+wr) & (x<mfs(j,5))).*...
		      (((mfs(j,5)-x)/wr).^2);
      end;
    end;
  end;
else
  y=zeros(k,lx); [md,nd] = size(domain);
  if md == 1, domain = ones(k,1)*domain; end;
  kmfs = ([mfs(:,2:nd) mfs(:,nd)] - mfs) ./ ...
          [domain(:,2:nd)-domain(:,1:nd-1) ones(k,1)];
  for j=1:k
    xx = x;
    maxd = max(domain(j,:)); mask = x > maxd; xx(mask) = maxd*mx(mask);
    mind = min(domain(j,:)); mask = x < mind; xx(mask) = mind*mx(mask);
    i = sum( (domain(j,:)')*ones(1,lx) <= ones(nd,1)*xx + eps );
    y(j,:) = mfs(j,i) + kmfs(j,i).*(xx - domain(j,i));
  end;
end;

