%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% For plotting polynomial and Gaussian-based 
% cluster functions (n=1 case).
%
% By: Kevin Passino
% Date: 3/24/99
%
% Note: This program provides some divide
% by zero errors, but the plots come out
% in an appropriate form.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set the cluster centers.

% Uniform spacing of clusters
v1=-5;
v2=0;
v3=5;

% Non-uniform spacing of clusters (uncomment out to get)
%v1=-3;
%v2=6;
%v3=1;

% Set cluster width parameter

m=2;

% Set sigma for the Gaussian functions

sigma1=1;
sigma2=0.7;  % For th nonuniform setting above uncomment out the next line
%sigma2=0.1;
sigma3=1;

% Define the domain 

x=-10:.1:10;

for i=1:length(x);
	% Begin with the polynomial functions
	p1(i)=(((x(i)-v1)^2/(x(i)-v1)^2)^(1/(m-1))+((x(i)-v1)^2/(x(i)-v2)^2)^(1/(m-1))+((x(i)-v1)^2/(x(i)-v3)^2)^(1/(m-1)))^(-1);
	p2(i)=(((x(i)-v2)^2/(x(i)-v1)^2)^(1/(m-1))+((x(i)-v2)^2/(x(i)-v2)^2)^(1/(m-1))+((x(i)-v2)^2/(x(i)-v3)^2)^(1/(m-1)))^(-1);
	p3(i)=(((x(i)-v3)^2/(x(i)-v1)^2)^(1/(m-1))+((x(i)-v3)^2/(x(i)-v2)^2)^(1/(m-1))+((x(i)-v3)^2/(x(i)-v3)^2)^(1/(m-1)))^(-1);
    mug1(i)=exp(-0.5*(x(i)-v1)^2/(sigma1^2));  % Calculate Gaussian
    mug2(i)=exp(-0.5*(x(i)-v2)^2/(sigma2^2));
    mug3(i)=exp(-0.5*(x(i)-v3)^2/(sigma3^2));
	sum(i)=mug1(i)+mug2(i)+mug3(i);
    mug1(i)=mug1(i)/sum(i);  % Normalize it
    mug2(i)=mug2(i)/sum(i);
    mug3(i)=mug3(i)/sum(i);
end

% Plot the results:

figure(1)
clf
subplot(211)
plot(x,p1,'k-',x,p2,'k--',x,p3,'k:')
grid
xlabel('x')
title('Polynomial-based cluster functions')

subplot(212)
plot(x,mug1,'k-',x,mug2,'k--',x,mug3,'k:')
grid
xlabel('x')
title('Gaussian-based cluster functions')

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