%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For training Takagi-Sugeno fuzzy systems using the
% Levenberg-Marquardt method
%
% By: Kevin Passino
% Version: 2/16/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% Specify the maximum number of iterations:
Nlm=15;

% First, generate the training data, G

% For the M=121 case

x=-6:0.1:6;

M=length(x)

for i=1:M,
	z(i)=0.15*(rand-0.5)*2;  % Define the auxiliary variable
	G(i)=exp(-50*(x(i)-1)^2)-0.5*exp(-100*(x(i)-1.2)^2)+atan(2*x(i))+2.15+...
	0.2*exp(-10*(x(i)+1)^2)-0.25*exp(-20*(x(i)+1.5)^2)+0.1*exp(-10*(x(i)+2)^2)-0.2*exp(-10*(x(i)+3)^2);
	if x(i) >= 0
		G(i)=G(i)+0.1*(x(i)-2)^2-0.4;
	end
	Gz(i)=G(i)+1*z(i); % Adds in the influence of the auxiliary variable (or can multiply by 0)
end

% Next, we provide another option for specifying the data set for M=13

%clear x G z Gz

%x=-6:1:6;

%M=length(x)

%for i=1:M,
%	z(i)=0.15*(rand-0.5)*2;  % Define the auxiliary variable
%	G(i)=exp(-50*(x(i)-1)^2)-0.5*exp(-100*(x(i)-1.2)^2)+atan(2*x(i))+2.15+...
%	0.2*exp(-10*(x(i)+1)^2)-0.25*exp(-20*(x(i)+1.5)^2)+0.1*exp(-10*(x(i)+2)^2)-0.2*exp(-10*(x(i)+3)^2);
%	if x(i) >= 0
%		G(i)=G(i)+0.1*(x(i)-2)^2-0.4;
%	end
%	Gz(i)=G(i)+1*z(i); % Adds in the influence of the auxiliary variable (or can multiply by 0)
%end


% Next, we need to initialize the appoximator

c(1,:,1)=-5:1:5; % The first 1 is the index for the first input (n=1 case)
						% and the second one is for the time index
						
R=length(squeeze(c(1,:,1)))

sigma(1,:,1)=0.5*ones(1,R);
%sigma(1,:,1)=0.2*ones(1,R);  % Another way to initialize
%sigma(1,:,1)=1*ones(1,R);  % Another way to initialize

% We initialize all the parameters that enter linearly to be zero
thetalin(:,1)=0*ones(2*R,1);

% With this the a_i,0 parameters are first in the thetalin vector,
% followed by the a_i,1 parameters, from i=1,2,...,R (thetalin
% is the vector of parameters that enter linearly)

% Next, form the whole parameter vector is

theta(:,1)=[c(1,:,1)'; sigma(1,:,1)'; thetalin(:,1)];

p=length(theta(:,1)); % Set the number of parameters for this case

% Next, pick the parameter for the update method (remember you want lambda
% bigger to make the updates *smaller*)

lambda=0.5;

Lambda=lambda*eye(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, we start the main loop:

for j=1:Nlm

	j % Print to the screen the interation number
	
% First, check if parameters are in range, and fix if not:

for jj=1:R
	if sigma(1,jj,j)<0.1
		sigma(1,jj,j)=0.1;
	end
	if sigma(1,jj,j)>1
		sigma(1,jj,j)=1;
	end
	if c(1,jj,j)>6
		c(1,jj,j)=6;
	end
	if c(1,jj,j)<-6
		c(1,jj,j)=-6;
	end
end


% Next, we need to form the vector (epsilon) and matrix (del_epsilon) that are used in the
% Levenberg-Marquardt update formula.  

for i=1:M

	% Define the phi vector, for each data pair in Gz

	for jj=1:R
		mu(jj,i,j)=exp(-0.5*((x(i)-c(1,jj,j))/sigma(1,jj,j))^2);
	end

	denominator(i)=sum(mu(:,i,j)); 

	xi(:,i,j)=mu(:,i,j)*(1/denominator(i));

	phi(:,i,j)=[xi(:,i,j)', x(i)*xi(:,i,j)']';

	Fts(i)=thetalin(:,j)'*phi(:,i,j);
	
	epsilon(i,j)=Gz(i)-Fts(i);

	for jj=1:R
		
		g(jj,i,j)=thetalin(jj,j)+thetalin(jj+R,j)*x(i);

		% First, compute it for the centers, then the spreads, then the consequent parameters

		temp=(g(jj,i,j)-Fts(i))/denominator(i);
		
		del_epsilon(jj,i,j)=-temp*mu(jj,i,j)*(x(i)-c(1,jj,j))*(1/(sigma(1,jj,j)^2));

		del_epsilon(jj+R,i,j)=-temp*mu(jj,i,j)*((x(i)-c(1,jj,j))^2)*(1/(sigma(1,jj,j)^3));
		
		del_epsilon(jj+2*R,i,j)=-xi(jj,i,j); % For the a_i,0
		
		del_epsilon(jj+3*R,i,j)=-x(i)*xi(jj,i,j);  % For the a_i,1
		
	end
end
		
	
% Next, we seek to generate the next estimate of the parameter vector theta for 
% the Levenberg-Marquardt update formula:

theta(:,j+1)=theta(:,j)-...
      inv(del_epsilon(:,:,j)*del_epsilon(:,:,j)'+Lambda)*del_epsilon(:,:,j)*epsilon(:,j);

% Load the components into the relevant vectors:

c(1,:,j+1)=theta(1:R,j+1)';
sigma(1,:,j+1)=theta((R+1):(2*R),j+1)';
thetalin(:,j+1)=theta((2*R+1):(4*R),j+1);

%end  % End of main loop for Levenberg-Marquardt update formula

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, compute the approximator values (i.e., at a test set)

% Define the points at which it will be tested

xt=-6:0.01:6; 
Mt=length(xt);

% First, define the phi vector, for each data pair in Gz
% (note that j+1 is the last time at which we have generated a theta
% and we use it here - first, we have to compute phi)

for i=1:Mt
	
	for jj=1:R
		mut(jj,i)=exp(-0.5*((xt(i)-c(1,jj,j+1))/sigma(1,jj,j+1))^2);
	end

	denominatort(i)=sum(mut(:,i)); 

	xit(:,i)=mut(:,i)*(1/denominatort(i));

	phit(:,i)=[xit(:,i)', xt(i)*xit(:,i)']';

	Fts(i)=thetalin(:,j+1)'*phit(:,i); % Just call the jth one the last one
end

% Next, plot the basis functions, data and the approximator to compare

%figure(1)
figure(j)
clf
plot(x,Gz,'ko',xt,xit,'k--',xt,Fts,'k-')
xlabel('x')
T=num2str(j);
T=strcat('Basis functions, data, and approximator mapping, j=',T);
ylabel(T)
title('Takagi-Sugeno fuzzy system, 11 rules')
grid
axis([min(x) max(x) 0 max(G)])

% Could omit the "end" on the main loop and add the following to see the mapping 
% shape at each iteration (in this case change figure(1) above to figure(j) if you want)
pause  % Can add this if you want to view the mappings at each iteration while
        % the program is running
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%