%This program uses 2_D datasets to check fuzzy C-means performance.
%Stop if cost function improvement over previous iteration is zero.

%Ref Roger Jang, 12-13-94
%Carson KAO, May 19, 1997

clear all;
case_num = 1;		%case number
cluster_n = 4;		%the numbers of cluster center
expo = 2;		%Exponent for U
max_iter = 1000;	%maximum iteration
display =0;		%display iteration counts or not
plot_opt = 0;		%plot data & center position or not
epoch_opt = 0;		%display epoch train & elapsed time or not
%min_impro = 1.0E-5;	%minimum amount of improvement
stop_impro = 0;		%stop amount of improvement
load_data = 0;		%load desired data file or default data pairs
save_center = 1;	%save center & Uij or not

t0 = clock;
%load data pair file
if load_data == 1,
   datafile = 'normal.dat';
   index = findstr(datafile, '.');
   eval(['load ', datafile]);
   eval(['data = ', datafile(1:index-1) ';']);
end

%load default dataset within the unit cube [-1,1]x[-1,1]
if load_data == 0,
   datafile = 'default.dat';
   index = findstr(datafile, '.');
   eval(['load ', datafile]);
   eval(['data = ', datafile(1:index-1) ';']);
end

[data_n data_dim] = size(data);
obj_fcn = zeros(max_iter, 1);	
epoch_n = zeros(case_num,1);
%=====================================
for m = 1:case_num,
 U = initfcm(cluster_n, data_n);			% Initial fuzzy partition
 [U, center, obj_fcn(1)] = stepfcm(data, U, cluster_n, expo);
 if plot_opt == 1,
   center_prev = center;
   figure;
   plot(data(:, 1), data(:, 2), 'go');
   headH = line(ones(2,1)*center(:,1)', ones(2,1)*center(:,2)',...
	'erase', 'none', 'linestyle', '.', ...
	'markersize', 30, 'tag', 'headH');
   trajH = line(zeros(2,cluster_n), zeros(2,cluster_n), ...
		'erase', 'none', 'linewidth', 3, ...
		'tag', 'trajH');
 end
 %main loop 
 for i = 1:max_iter,
    [U, center, obj_fcn(i)] = stepfcm(data, U, cluster_n, expo);
    if display, 
	fprintf('Iteration count = %d, obj. fcn = %f\n', i, obj_fcn(i));
    end
    % check termination condition
    if i > 1,
	%if abs(obj_fcn(i) - obj_fcn(i-1)) < min_impro, break; end,
        if abs(obj_fcn(i) - obj_fcn(i-1)) == stop_impro, break; end,
    end
    if plot_opt == 1,
       for j = 1:cluster_n,
	   set(headH(j), 'xdata', center(j, 1), 'ydata', center(j, 2));
           set(trajH(j), 'xdata', [center_prev(j, 1) center(j, 1)], ...
			'ydata', [center_prev(j, 2) center(j, 2)]);
       end
       drawnow;
       center_prev = center;			
    end
 end
 iter_n = i;	% Actual number of iterations 
 obj_fcn(iter_n+1:max_iter) = [];
 epoch_n(m) = length(obj_fcn);
 %show clustering centers position
 if plot_opt == 1,
   figure;
   plot(data(:, 1), data(:, 2), 'go');
   headH = line(ones(2,1)*center(:,1)', ones(2,1)*center(:,2)',...
	'erase', 'none', 'linestyle', '.', ...
	'markersize', 30, 'tag', 'headH');
   for j = 1:cluster_n,
	set(headH(j), 'xdata', center(j, 1), 'ydata', center(j, 2));
   end
 end
 [tmp index_center] = sort(center(:,1));
 center = center(index_center,:);
 U = U(index_center,:);
 if plot_opt == 1,
    fprintf('FCM');
    center
 end
 eval(['center',int2str(m),'=center;']);	%save center_case1,case_2
 eval(['U',int2str(m),'=U;']);			%save U_case1,case_2

 %find the amount of data sets in cluster C
 maxU = max(U);
 for k = 1:cluster_n,
       tk = int2str(k);
       eval(['index',tk,'=find(U(k,:) == maxU);']);
       if plot_opt == 1,
          eval(['fprintf(''FCM Cluster %d contains %d datasets\n'',k,length(index'...
          ,tk,'));']);
       end
 end
 %point out all data pairs belong to which group in different color
 if plot_opt == 1,
   figure;
   plot(data(:, 1), data(:, 2), 'o');
   alllinestyle = ['+';'*';'o';'x'];
   allcolor = ['r';'y';'m';'c';'b';'g'];
   %line(data(index1,1),data(index1,2), 'linestyle','alllinestyle(1)','color','r');
   line(data(index1,1),data(index1,2), 'linestyle','o','color','r');
   line(data(index2,1),data(index2,2), 'linestyle','o','color','y');
   line(data(index3,1),data(index3,2), 'linestyle','o','color','m');
   line(data(index4,1),data(index4,2), 'linestyle','o','color','c');
 end
end		%end for

if epoch_opt == 1,
   ave_epoch = sum(epoch_n)/case_num;
   fprintf('Average epoch_n = %.2f epochs\n',ave_epoch);
   elapsed_time = etime(clock,t0);
   fprintf('Elapsed time = %f\n',elapsed_time);
   figure;
   plot(1:case_num, epoch_n, '-', 1:case_num, epoch_n, 'o');
   ylabel('Number of Epochs');
end
center_sum = 0;
U_sum = zeros(size(U));
for i = 1:case_num,
    eval(['center_sum = center_sum+center',int2str(i),';']);
    eval(['U_sum = U_sum+U',int2str(i),';']);
end 
center_avg = center_sum/case_num;
U_avg = U_sum/case_num;
%save cluster center & cost function
if save_center == 1,
   if load_data == 1,
      save fncentr.dat center_avg -ascii;
      save fncost.dat U_avg -ascii;
   else
      save fdcentr.dat center_avg -ascii;
      save fdcost.dat U_avg -ascii;
   end
end
