extend_binary(x1,mem1,x2,mem2,op,tnorm,step)

%% extend_binary.m

%% This function uses Zadeh's Extension Principle (Section 7.2 in the 
%% book Uncertain Rule-Based Fuzzy Logic Systems: Introduction and 
%% New Directions, by Jerry M. Mendel, and published by Prentice-Hall, 2000.) 
%% to extend a binary operation "op" to two type-1 fuzzy sets. 
 
%% Written by Nilesh N. Karnik - July 20,1998
%% For use with MATLAB 5.1 or higher.

%% Outputs : "out" and "mem" (column vectors) are, respectively, 
%% the domain and the memberships of the result of the extended 
%% "op" operation. 

%% Inputs : "x1" and "mem1" (same size vectors) are, respectively, 
%% the domain and  memberships of the first type-1 fuzzy set; 
%% and, "x2" and "mem2" (same size vectors) are, respectively, 
%% the domain and memberships of the second type-1 fuzzy set. "op" 
%% is a string (entered in inverted commas, as " 'op' ") containing 
%% the name of the function which defines the binary operation 
%% to be extended (the function should be stored in the file 
%% "op.m"). If "tnorm" < 0 (scalar), minimum t-norm is used, 
%% otherwise product t-norm is used. "step" (scalar) is an 
%% optional parameter. If it is not provided, the two type-1 sets 
%% are assumed to have discrete domains and the output is presented 
%% just by sorting the result of the extened "op" operation. 
%% If "step" is provided, it is assumed that the two type-1 
%% sets actually have continuous domains and are discretized for 
%% computing purposes. The resulting type-1 set, in this case, has 
%% equispaced domain points : "[minx : step : maxx]", where "minx" and 
%% "maxx" are, respectively, the minimum and the maximum of the results 
%% of the "op" operation between domain elements of x1 and x2. 

%% Note 1 : The t-conorm used is maximum.
%% Note 2 : If "op" is a built-in function, you may have to use the 
%% MATLAB function "builtin" rather than "feval".
%% Note 3 : The function "op" should be written so that it can work 
%% with matrices. For an example, see the function "mult.m".
 
%% Uses function "trimvec.m" .

function[out,mem] = extend_binary(x1,mem1,x2,mem2,op,tnorm,step)

   if size(x1,2) > 1,
      x1 = x1' ;
      mem1 = mem1' ;
   end

   if size(x2,2) > 1,
      x2 = x2' ;
      mem2 = mem2' ;
   end

   o1 = ones(size(x1)) ;
   o2 = ones(size(x2)) ;
   
   ox1 = x1*o2' ;
   ox2 = o1*x2' ;

   omem1 = mem1*o2' ;
   omem2 = o1*mem2' ;

   x3 = feval(op,ox1,ox2) ;

   if tnorm >= 0, 
      mem3 = omem1.*omem2 ;   % product t-norm
   else
      mem3 = min(omem1,omem2) ;   % minimum t-norm
   end   %%  if tnorm


if nargin == 6,

    lx1 = length(x1) ;
    lx2 = length(x2) ;

    for i1 = 1 : lx1,
        x4(lx2*(i1-1)+1 : lx2*i1) = x3(i1,:) ;
        mem4(lx2*(i1-1)+1 : lx2*i1) = mem3(i1,:) ;
    end   %%% for i1
           
    [out,mem] = trimvec(x4,mem4) ;

else 
    
    minx = min(min(x3)) ;
    maxx = max(max(x3)) ;    
   
    eps = step/2 ;
    i2 = 1 ;
    for k = minx : step : maxx,
        out(i2) = k ;
        [in1,in2] = find(abs(x3-k) <= eps) ;

        if ~isempty(in1) & ~isempty(in2),
            mem(i2) = max(diag(mem3(in1,in2))) ;
        else 
            mem(i2) = mem(i2-1) ;
        end   %%% if 
        
        i2 = i2 + 1;
    end   %% for k

    if k ~= maxx,
        out(i2) = maxx ;
        [in1,in2] = find(abs(x3-maxx) <= eps) ;

        if ~isempty(in1) & ~isempty(in2),
            mem(i2) = max(diag(mem3(in1,in2))) ;
        else 
            mem(i2) = mem(i2-1) ;
        end   %%% if ~isempty
    end   % if k      


end  % if nargin

return ;