%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Response surface method for finding optimal approximator size p
% and data set G size M for a function approximation problem.
% Program also allow you to keep either p or M constant and 
% vary the other one (to generate a 2D response surface).
%
% By: Kevin Passino
% Version: 3/22/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% Consider the follow range of sizes of the approximator
p=4:2:40; % Notice that we consider only even p for convenience since
			% it is easiest to add two parameters at a time for the 
			% type of approximator that we are using (TS fuzzy system)
			% Also, must start at p(1)>=4 due to how the spreads are defined
			% in terms of the centers.

% Consider the following range of sizes of the data set
M=100:100;

% Generate a test set

% Define a parameter to control the size of the test set
% relative to sizes of training sets

gamma=2; % Choose the size of the test set to be twice the size
		 % of the largest training set
		 
[xt,Gamma]=unknownfunction(gamma*max(M));

% Set the number of trials for each point on the RSM
% (to use an average MSE for each point on the RSM)

Ntrials=1; 
%Ntrials=100; 

% Initialize the variables that hold the training data:
x=0*ones(length(M),length(p),length(M),Ntrials);
G=0*ones(length(M),length(p),length(M),Ntrials);

% Set a variable for the number of times to run the program (for checking
% that the response surface does not change too much for additional trials)

Ntests=1;
Ntests=5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop for generating data for RSM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ll=1:Ntests
	
for ii=1:length(p)
	for jj=1:length(M)
		
		for kk=1:Ntrials
			
			[x(1:M(jj),ii,jj,kk),G(1:M(jj),ii,jj,kk)]=unknownfunction(M(jj)); % Generate training data
			
			% First, form the vector Y

			Y=G(1:M(jj),ii,jj,kk);

			% Find Phi, which involves processing x through phi (for the training data)

			% Find the centers of the membership functions
			c(1,:)=-6:(12/(0.5*p(ii)-1)):6; % Produces one center for every two parameters 
			sigma(1,:)=((c(1,2)-c(1,1))/(2))*ones(1,length(c(1,:))); % Chooses width to be half the distance
															% between adjacent centers

			% Initialize Phi

			for j=1:length(c(1,:))
				mu(j,1)=exp(-0.5*((x(1,ii,jj,kk)-c(1,j))/sigma(1,j))^2);
			end
	
			denominator(1)=sum(mu(:,1)); 
	
			for j=1:length(c(1,:))
				xi(j,1)=mu(j,1)/denominator(1);
			end

			Phi=[xi(:,1)', x(1,ii,jj,kk)*xi(:,1)'];

			% Form the rest of Phi

			for i=2:M(jj)
	
				for j=1:length(c(1,:))
					mu(j,i)=exp(-0.5*((x(i,ii,jj,kk)-c(1,j))/sigma(1,j))^2);
				end
	
				denominator(i)=sum(mu(:,i)); 
	
				for j=1:length(c(1,:))
					xi(j,i)=mu(j,i)/denominator(i);
				end

				Phi=[Phi; xi(:,i)', x(i,ii,jj,kk)*xi(:,i)'];
	
			end

			% Find the least squares estimate

			theta=Phi\Y;

		
			% Compute the approximator values for the test set (xt,Gamma)

			for i=1:length(Gamma)
	
				for j=1:length(c(1,:))
					mut(j,i)=exp(-0.5*((xt(i)-c(1,j))/sigma(1,j))^2);
				end
	
				denominatort(i)=sum(mut(:,i)); 
	
				for j=1:length(c(1,:))
					xit(j,i)=mut(j,i)/denominatort(i);
				end	
			
				phit=[xit(:,i)', xt(i)*xit(:,i)']';
				Fts(i,1)=theta'*phit;			
			
			end
		
			% Plot the data (optional - just to gain insight)
			
%			figure(1)
%			clf
%			plot(xt,Fts,'o')
%			ylabel('Approximator output F_t_s')
%			xlabel('Test data')
%			Tp=num2str(p(ii));
%			TM=num2str(M(jj));
%			T=strcat('Approximator mapping, p=',Tp,'Training set size, M=',TM);
%			title(T)
%			pause
	
			% Compute the approximation error at the test data
		
			J(jj,ii,kk)=(1/length(Gamma))*(Gamma-Fts(:,1))'*(Gamma-Fts(:,1)); % Compute the mean-squared error
		
		end % End Ntrials loop
	
		% Need to clear variables as dimensions change in two outer-most loops
		clear Y c sigma mu denominator Phi xi theta mut denominatort xit phit
		
		% Compute the average mean squared error
		Jsurf(jj,ii,ll)=(1/Ntrials)*sum(J(jj,ii,:));

end % End M loop
end % End p loop
end % Ends the Ntests loop


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the response surface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)
clf

% First consider 2D surfaces

if length(M)==1, % Plot across p

for ll=1:Ntests
	
figure(2)
plot(p,squeeze(Jsurf(1,:,ll)),'k-')
xlabel('Approximator size, p');
ylabel('Average MSE, J');
TM=num2str(M(1));
T=strcat('Response surface, Training set size, M=',TM);
title(T)
% Next line for zooming in on response surface region
%axis([10 max(p) min(Jsurf(1,:,ll)) 0.03])
hold on

end

end

%clear J 

% Below, some additional plotting possibilities

%if length(p)==1, % Plot across M
	
%figure(3)
%clf
%plot(M,Jsurf(:,1),'k-')
%xlabel('Training data set size, M');
%ylabel('Average MSE, J');
%Tp=num2str(p(1));
%T=strcat('Response surface, Approximator size, p=',Tp);
%title(T)

%end

% Next, 3D

%if length(M)>1 & length(p)>1, % If have data for a 3D surface 

%figure(4)
%clf
%surf(p,M,Jsurf)
%colormap(white);
%xlabel('Approximator size, p');
%ylabel('Training data set size, M');
%zlabel('Average MSE, J');
%title('Response surface for average MSE');
%rotate3d

%end

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