function DOF = dofprod(rls,mfs,x,domain)
% DOFPROD calculates degrees of fulfillment (DOF's) of the premises,
% using product as the fuzzy logic AND operator.
%
% DOF = DOFPROD(RLS,MFS,X,DOM)
% Inputs:       RLS - rule matrix
%               MFS - membership function matrix
%               X   - input vector or matrix
%               DOM - (optional) domain for membership function
%                     defined point-wise
%
% Output:       DOF - matrix of DOF's
%
% See also TXT2RLS, RLS2TXT, MFSET, MFGET, INFER, DOFMIN, DOFTN.

% (c) Robert Babuska, December 1993

if nargin < 4, domain = []; end;
[k,ninps] = size(x);
[nrules,dummy] = size(rls);
[nmfs,npoints] = size(mfs);
if length(domain) == 0; mfany = [0 0 0 0 0];            % continuous domain
else 							% discrete domain
   mfany = ones(1,npoints); 
   domain = [domain; randn(1,npoints)]; 
end;                    
mfs = [mfs; mfany]; nmfs = nmfs + 1;
DOF = zeros(nrules,k);                                  % transposed DOF
w = rls(:,1); rls = rls(:,2:ninps+1)';                  % extract weights
aux = ones(ninps,k*nrules);                             % for rls expansion
grades = aux';                                          % prepare for grades
aux(:) = rls(:)*ones(1,k);                              % expand rls
rls = zeros(nrules,k*ninps);                            % ....
rls(:) = aux';                                          % ....
ANYmask = rls == 0;                                     % find zeros (ANY)
rls(ANYmask) = nmfs*grades(ANYmask);                    % index for ANY
mask = rls < 0;                                         % find NOT in rules
sgn = -mask + (rls >= 0);                               % not(A) = 1 - A
shift = ones(nrules,1)*[0 nmfs*(1:k*ninps-1)];          % find index shift
rls = abs(rls) + shift;
g = mgrade(x(:)',mfs,domain);                           % calculate all grades
grades(:) = mask + sgn .* g(rls);                       % select using rls
DOF(:) = prod([ones(1,k*nrules);grades']);              % calculate product
DOF = (ones(k,1)*w').*DOF';                             % weight by w

