weighted_avg(Zout,Zmem,Wout,Wmem,tnorm,step)

%% weighted_avg.m 

%% Function to extend the weighed average of crisp numbers to
%% the case where all the quantities involved are type-1 sets.
%% It implements the generalized centroid in Equation (9-11) in 
%% the book Uncertain Rule-Based Fuzzy Logic Systems: Introduction 
%% and New Directions, by Jerry M. Mendel, and published by 
%% Prentice-Hall, 2000. 

%% Written by Nilesh N. Karnik - July 22,1998
%% For use with MATLAB 5.1 or higher.

%% NOTE : This program uses an algorithm similar to 
%% "extend_nary2.m"; therefore, it is slow, but uses less memory
%% than "extend_nary1.m".
 
%% Outputs : "out" and "mem" (column vectors) are, respectively, 
%% the domain and the memberships of the result of the extended 
%% weighted average operation. 
 
%% Inputs : "Zout" and "Zmem" are "M x Nz" matrices, containing the 
%% domain elements  and memberships of the "M" type-1 sets "Z_i"s - see 
%% Eq. (9-11). "Wout" and "Wmem" are "M x Nw" matrices, containing 
%% the domain elements and memberships of the "M" type-1 "W_i"s - see 
%% Eq. (9-11). The domain of each "Z_i" is assumed to contain "Nz" 
%% elements and that of each "W_i" is assumed to contain "Nw" elements. 
%% If "tnorm" (scalar) is < 0, minimum t-norm is used, otherwise 
%% product t-norm is used. Here, "step" (scalar) is a required 
%% parameter. In order to save memory, the domain of the output type-1 
%% set is taken to be "[minz : step : maxz]", where 
%% "minz = min(min(Zout))" and "maxz = max(max(Zout))". Note that 
%% since this is a weighted average of the "Z_i"s, the smallest possible 
%% domain point in the output set is "minz" and the largest possible is 
%% "maxz".
 
%% Note 1 : The t-conorm used is maximum.
%% Note 2 : The notation in this program is somewhat different than 
%% that in "extend_nary1.m" and "extend_nary2.m", but it conforms with
%% the notation in Eq. (9-11) .
%% Note 3 : Since type-reduction calculations can be represented as 
%% extended weighted average calculations, this program can be used 
%% for type-reduction calculations in Chapter 9. 

function[out,mem] = weighted_avg(Zout,Zmem,Wout,Wmem,tnorm,step)

[M,Nz] = size(Zout) ;
Nw = size(Wout,2) ;

Zout1 = Zout ;
Zmem1 = Zmem ;
Wout1 = Wout ;
Wmem1 = Wmem ;

if tnorm < 0,
   tnorm_op = 'min' ;
else
   tnorm_op = 'prod' ;
end    %%% if tnorm

                                %% Points for the output
minz = min(min(Zout)) ;
maxz = max(max(Zout)) ;

out = [minz : step : maxz] ;
mem = zeros(size(out)) ;

count_Z = Nz^(M-1) ;
count_W = Nw^(M-1) ;
temp1 = zeros(1,Nz) ;
temp2 = temp1 ;
temp3 = zeros(1,Nw) ;
temp4 = temp3 ;
temp5 = temp3 ;


for i6 = 1 : count_Z,

 for i4 = 1 : Nz,
    Zpoint = Zout1(:,i4)*ones(1,Nw) ;
    memZ = feval(tnorm_op,Zmem1(:,i4)) ;
    num_mat = Zpoint .* Wout1 ;
    
    for i1 = 1 : count_W,
        answer = sum(num_mat)./sum(Wout1) ;    
        ind = round((answer - out(1))/step) + 1 ;
    
        memW = feval(tnorm_op,Wmem1) ;
        if tnorm < 0,,
           memout = min(memW,memZ) ;
        else
           memout = memZ.*memW ;
        end  %% if minflag
    
        for i3 = 1 : Nw,
           mem(ind(i3)) = max(mem(ind(i3)) , memout(i3)) ;   % max t-conorm
        end  %% if i3

        temp5(1:Nw-1) = num_mat(M,2:Nw) ;
        temp5(Nw) = num_mat(M,1) ;
        num_mat(M,:) = temp5 ;    
    
        temp3(1:Nw-1) = Wout1(M,2:Nw) ;
        temp3(Nw) = Wout1(M,1) ;
        Wout1(M,:) = temp3 ;    
        
        temp4(1:Nw-1) = Wmem1(M,2:Nw) ;
        temp4(Nw) = Wmem1(M,1) ;
        Wmem1(M,:) = temp4 ;    
    
        for i2 = 1 : M-1,
            if mod(i1,Nw^i2) == 0,
    
               temp5(1:Nw-1) = num_mat(M-i2,2:Nw) ;
               temp5(Nw) = num_mat(M-i2,1) ;
               num_mat(M-i2,:) = temp5 ;    

               temp3(1:Nw-1) = Wout1(M-i2,2:Nw) ;
               temp3(Nw) = Wout1(M-i2,1) ;
               Wout1(M-i2,:) = temp3 ;    
               
               temp4(1:Nw-1) = Wmem1(M-i2,2:Nw) ;
               temp4(Nw) = Wmem1(M-i2,1) ;
               Wmem1(M-i2,:) = temp4 ;    
           end %% if mod(i1,Nw^i2)
        end   %%% for i2
    
    end   %%% for i1
    
  end %% for i4
    
    temp1(1:Nz-1) = Zout1(M,2:Nz) ;
    temp1(Nz) = Zout1(M,1) ;
    Zout1(M,:) = temp1 ;    
       
    temp2(1:Nz-1) = Zmem1(M,2:Nz) ;
    temp2(Nz) = Zmem1(M,1) ;
    Zmem1(M,:) = temp2 ;    
    
    
    for i5 = 1 : M-1,
        if mod(i6,Nz^i5) == 0,
 
           temp1(1:Nz-1) = Zout1(M-i5,2:Nz) ;
           temp1(Nz) = Zout1(M-i5,1) ;
           Zout1(M-i5,:) = temp1 ;    
               
           temp2(1:Nz-1) = Zmem1(M-i5,2:Nz) ;
           temp2(Nz) = Zmem1(M-i5,1) ;
           Zmem1(M-i5,:) = temp2 ;    
    
        end %% if mod(i1,N^i5)
    end   %%% for i5

end %%% for i6


return ;