function plot(m,varargin)
%PLOT performs a plot of the a taksug model
%    
%    F=PLOT(M)  opens a new figure if needed and plots the
%    model M. A complete user interface is set to allow
%    advanced visualisation options to be selected.
%    
%    F=PLOT(M,OPTIONS)  opens a new figure if needed and plots
%    the model M. OPTIONS are passed to the drawing engine.
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose performs a plot of the a taksug model.
%
%   @synopsis f=plot(m) 
%   @description opens a new figure if needed and plots the model
%   <CODE>m</CODE>. A complete user interface is set to allow advanced visualisation
%   options to be selected.
%
%   @synopsis f=plot(m,options) 
%   @description opens a new figure if needed and plots the model
%   <CODE>m</CODE>. <CODE>options</CODE> are passed to the drawing engine.
%
%MAN_PAGE_END        

% By Antoine Duchateau. All Rights Reserved.
% Written: 7/06/99

ni = nargin;
if (ni<1)|(ni/2 == floor(ni/2)),
	error('Wrong number of inputs');
end

f = gcf;
cura = gca;
cla;

ud = get(f,'userdata');
udbis = get(cura, 'userdata');

udbis.m = m;
if ~isfield(udbis, 'secaxes'), udbis.secaxes = []; end

% Defaults values for parameters
select_p = 1:6;
axis_p = [];
data = [];

if get(m,'n_in') <= 1,
	subset_p = [1 1];
else
	subset_p = [1 2 1];
end

g_p = zeros(get(m,'n_in'),1);

for i=1:(ni - 1)/2,
	name = varargin{i*2-1};
	value = varargin{i*2};
	
	switch name,
			
	case 'data',
		data = value;		

	case 'g',
		g_p = value;		

	case 'select',
		select_p = value
		
	case 'axis',
		axis_p = value;		

	case 'subset',
		subset_p = value;
		if subset_p(1) > get(m,'n_in'), subset_p(1) = []; end
		if subset_p(end-1) > get(m,'n_in'), subset_p(end-1) = []; end
	end
end

if length(subset_p) == 2,
	
	win = subset_p(1);
	wout = subset_p(2);
	
	lims = get(m,'limits');
	theAx1 = lims(:,[win wout+get(m, 'n_in')]);
	if isempty(data),
		theAx1(find(theAx1 == inf)) = 1;
		theAx1(find(theAx1 == -inf)) = -1;
	else
		theAx1(find(~isfinite(theAx1))) = 0;
		theAx1(1,:) = min(theAx1(1,:),min(data(:,[win wout+get(m, 'n_in')])));
		theAx1(2,:) = max(theAx1(2,:),max(data(:,[win wout+get(m, 'n_in')])));
	end
	
	theRanges=sum(theAx1);
	theMinMax=[theAx1(1,:)-0.1*theRanges;theAx1(2,:)+0.1*theRanges];
	
	ax=[];
	
	if ~isempty(axis_p),
		ax = axis_p;
	end
	
	if isempty(ax)|(size(ax,2) ~= 2 * length(subset_p)),
		ax = theMinMax(:)';
		signif = round(ax) + (round(ax) == 0);
		ax = round(ax./signif*100) .* signif/100;
		if isfield(ud,'controls'),
			set(ud.controls.axis,'String',mat2str(ax,3));
		end
	end;
	
	if any(select_p == 6),
		try, PlotData(m, subset_p, ax, g_p, data); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 1), 
		try, PlotModel(m, subset_p, ax, g_p); catch, end
		axis(ax(:)');
	end
		
	if any(select_p == 3),
		try, PlotClus(m, subset_p); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 4),
		try, PlotConseq(m, subset_p, ax, g_p); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 2), 
		if ~(isfield(udbis.secaxes, 'memb')),
			udbis.secaxes.memb = axes(	'Position', get(cura,'Position'), ...
										'YAxisLocation', 'right', ...
										'ylim',[0 2], ...
										'SelectionHighlight', 'on', ...
										'XTick',0, ...
										'XLim', get(cura,'XLim'));
		end
		
		set(cura,'Color','none');
		set(udbis.secaxes.memb,'visible','on');
		axes(udbis.secaxes.memb);
		cla;
		try, PlotAntec(m, subset_p, ax, g_p); catch, end
		set(udbis.secaxes.memb,'XLim', get(cura,'XLim'));
		axes(cura);
	else
		if (isfield(udbis.secaxes, 'memb')),
			axes(udbis.secaxes.memb);
			cla;
			set(udbis.secaxes.memb,'visible','off');
			axes(cura);
			set(cura,'Color',[1 1 1]);
		end
	end
	
	if any(select_p == 5), 
		if ~(isfield(udbis.secaxes, 'error')),
			siz = get(cura,'Position');
			set(cura,'position',siz.*[1 1 1 0.48] + [0 siz(4)*0.52 0 0])
			if (isfield (udbis.secaxes, 'memb')),
				set(udbis.secaxes.memb, 'Position', get(cura,'Position'));
			end
			udbis.secaxes.error = axes(	'Position', siz.*[1 1 1 0.48], ...
										'SelectionHighlight', 'on');
			udbis.secaxes.oldsize = siz;
		end
		
		cura = cura;
		axes(udbis.secaxes.error);
		try, PlotError(m, subset_p, ax, g_p, data); catch, end
		set(udbis.secaxes.error,'XLim', get(cura,'XLim'));
		axes(cura);
	else
		% Need to clean up
		if (isfield(udbis.secaxes, 'error')),
			set(cura, 'position', udbis.secaxes.oldsize);
			if (isfield (udbis.secaxes, 'memb')),
				set(udbis.secaxes.memb, 'Position', udbis.secaxes.oldsize);
			end
			delete(udbis.secaxes.error);
			udbis.secaxes = rmfield(udbis.secaxes, 'error');
		end
	end

		
elseif length(subset_p) == 3,
	
	win = subset_p(1:2);
	wout = subset_p(3);
	
	lims = get(m,'limits')
	theAx1 = lims(:,[win wout+get(m, 'n_in')]);
	if isempty(data),
		theAx1(find(theAx1 == inf)) = 1;
		theAx1(find(theAx1 == -inf)) = -1;
	else
		theAx1(find(~isfinite(theAx1))) = 0;
		theAx1(1,:) = min(theAx1(1,:),min(data(:,[win wout+get(m, 'n_in')])));
		theAx1(2,:) = max(theAx1(2,:),max(data(:,[win wout+get(m, 'n_in')])));
	end
	
	theRanges=sum(theAx1);
	theMinMax=[theAx1(1,:)-0.1*theRanges;theAx1(2,:)+0.1*theRanges];
	
	ax=[];
	
	if ~isempty(axis_p),
		ax = axis_p;
	end
	
	if isempty(ax)|(size(ax,2) ~= 2 * length(subset_p)),
		ax = theMinMax(:)';
		signif = round(ax) + (round(ax) == 0);
		ax = round(ax./signif*100) .* signif/100;
		if isfield(ud,'controls'),
			set(ud.controls.axis,'String',mat2str(ax,3));
		end
	end;
	axis(ax(:)');
	
	if any(select_p == 6),
		try, PlotData3(m, subset_p, ax, g_p, data); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 1), 
		try, PlotModel3(m, subset_p, ax, g_p); catch, end
		axis(ax(:)');
	end
		
	if any(select_p == 2), 
		try, PlotAntec3(m, subset_p, ax, g_p); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 3),
		try, PlotClus3(m, subset_p); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 4),
		try, PlotConseq3(m, subset_p, ax, g_p); catch, end
		axis(ax(:)');
	end
	
	if any(select_p == 5), 
		try, PlotError3(m, subset_p, ax, g_p); catch, end
		axis(ax(:)');
	end	
	
	
end


% Deal with axes dependencies
cbk = '';
if (isfield(udbis.secaxes, 'error')),
	cbk = [cbk 'delete(getfield(getfield(get(gcbo,''userdata''),''secaxes''),''error''));'];
end
if (isfield(udbis.secaxes, 'memb')),
	cbk = [cbk 'delete(getfield(getfield(get(gcbo,''userdata''),''secaxes''),''memb''));'];
end
set(cura,'deletefcn',cbk);


set(f,'userdata',ud);
set(cura,'userdata',udbis);


function PlotClus(m, subset_p);
if isfield(get(m,'opt'),'ClusVars')&isfield(get(m,'opt'),'ClusCenters'),
	hold on;

	numrules=m.n_rules;
	
	x(1,:)=cos([0:0.1:2*pi+0.1]);
	x(2,:)=sin([0:0.1:2*pi+0.1]);
	
	opt = get(m,'opt');
	
	for i=1:numrules,
		theInv=projclus(inv(opt.ClusVars(:,:,i)),[subset_p(1);get(m,'n_in')+1]);
		
		[theRot,theHom]=eig(theInv);
		theHom=inv(sqrt(theHom));
		theRot=theRot*det(theRot);
		
		y=opt.ClusCenters([subset_p(1) get(m,'n_in') + 1],i) * ones(1,size(x,2)) + 2*theRot*theHom*x;
		
		plot(y(1,:),y(2,:),'r');
	end
	hold off;
end


function PlotAntec (m, subset_p, r, g);
numrules=m.n_rules;

hold on;

x= ((0:0.01:1) * (r(2)-r(1)) + r(1));
if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:subset_p(1)-1,:);x;theSec(subset_p(1)+1:size(theSec,1),:)];
end

y=membership(m,theInp')';

if strcmpi(m.model_code{2},'gaussian')|strcmpi(m.model_code{2},'inversedist'),
	sumY=sum(y,1)+realmin;
	y=y./sumY(ones(1,size(y,1)),:);
end

for i=1:numrules,
	plot(x,y(i,:),'b');
end

hold off;

function PlotData(m, subset_p, r, g, data);

hold on;

plot(data(:,subset_p(1))',data(:,subset_p(2)+get(m, 'n_in'))','y+');

hold off;

function PlotModel(m, subset_p, r, g);

hold on;

x=((0:0.01:1)*(r(2)-r(1))+r(1));
if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:subset_p(1)-1,:);x;theSec(subset_p(1)+1:size(theSec,1),:)];
end
y=eval(m,theInp',subset_p(2))';
plot(x,y,'g');

hold off;


function PlotError (m, subset_p, r, g, data);
x = data(:,1)';
[x, idx] = sort(x);

if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:subset_p(1)-1,:);x;theSec(subset_p(1)+1:size(theSec,1),:)];
end
y=(eval(m,theInp',subset_p(2)) - data(idx,2))';
plot(x,y,'b');


function PlotConseq (m, subset_p, r, g);
hold on;

numrules=m.n_rules;
x=((0:0.01:1)*(r(2)-r(1))+r(1));
if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:subset_p(1)-1,:);x;theSec(subset_p(1)+1:size(theSec,1),:)];
end

memb=membership(m,theInp')';
[maxis,ids]=max(memb);

lims = find(ids(1:end-1) - ids(2:end) ~= 0);
lims = [0 lims length(ids)];

for i=2:length(lims),
	sub_x=x(:,[lims(i-1)+1 lims(i)])';
	sub_inp=theInp(:,[lims(i-1)+1 lims(i)])';
	try,
		y=[sub_inp ones(size(sub_x,1),1)] * m.linears(subset_p(2),:,ids(lims(i)))';
		plot(sub_x,y,'r');
	catch
	end
end

hold off;

%function for plotting the model error in 3D plots
function PlotError3(m, subset_p, r, g, data);
%
%

%function for plotting the clusters in 3D
function PlotClus3(m, subset_p, r, g);

hold on;

numrules=m.n_rules;
ni = get(m,'n_in');
[theX,theY,theZ]=sphere(50);
theJoin=[theX(:)';theY(:)';theZ(:)'];

theLocX=theX;
theLocY=theY;
theLocZ=theZ;

opt = get(m,'opt');
for i=1:numrules,

	theInv=projclus(inv(opt.ClusVars(:,:,i)),[subset_p(1);subset_p(2);subset_p(3)+ni]);
	[theRot,theHom]=eig(theInv);
	theHom=inv(sqrt(theHom));
	theRot=theRot*det(theRot);
	theJoinTrans=opt.ClusCenters([subset_p(1) subset_p(2) subset_p(3)+ni],i)*ones(1,size(theJoin,2))+2*theRot*theHom*theJoin;
	theLocX(:)=theJoinTrans(1,:);
	theLocY(:)=theJoinTrans(2,:);
	theLocZ(:)=theJoinTrans(3,:);
	theCluster=surf(theLocX,theLocY,theLocZ);
	set(theCluster,'EdgeColor','none');
end

colormap(gray);
light('Position',[0 -1 2]);
hold off;

function PlotModel3(m, subset_p, r, g);

hold on;
[x,y]=meshgrid(r(1):(r(2)-r(1))/40:r(2), r(3):(r(4)-r(3))/40:r(4));
if isempty(g),
	theInp=[x(:)';y(:)'];
else
	theSec=g(:,ones(1,size(x,1)*size(x,2)));
	if subset_p(2)>subset_p(1),
		theInp=[theSec(1:subset_p(1)-1,:);x(:)';theSec(subset_p(1)+1:subset_p(2)-1,:);y(:)';theSec(subset_p(2)+1:size(theSec,1),:)];
	else
		theInp=[theSec(1:subset_p(2)-1,:);y(:)';theSec(subset_p(2)+1:subset_p(1)-1,:);x(:)';theSec(subset_p(1)+1:size(theSec,1),:)];
	end
end;
z=eval(m,theInp',subset_p(3))';
theSurf=surf(x,y,reshape(z,size(x)));
set(theSurf,'EdgeColor','none');
colormap(hsv);
hold off;

function PlotConseq3(m, subset_p, r, g);
hold on;
numrules=m.n_rules;

[x,y]=meshgrid(r(1):(r(2)-r(1))/40:r(2), r(3):(r(4)-r(3))/40:r(4));

if isempty(g),
	theInp=[x(:)';y(:)'];
else
	theSec=g(:,ones(1,size(x,1)*size(x,2)));
	if subset_p(2)>subset_p(1),
		theInp=[theSec(1:subset_p(1)-1,:);x(:)';theSec(subset_p(1)+1:subset_p(2)-1,:);y(:)';theSec(subset_p(2)+1:size(theSec,1),:)];
	else
		theInp=[theSec(1:subset_p(2)-1,:);y(:)';theSec(subset_p(2)+1:subset_p(1)-1,:);x(:)';theSec(subset_p(1)+1:size(theSec,1),:)];
	end
end;
memb=membership(m,theInp')';

for i=1:numrules,
	[maxis,ids]=max(memb);
	sub_x=x(find(ids == i))';
	sub_y=y(find(ids == i))';
	[min_x,min_x_i]=min(sub_x);
	[max_x,max_x_i]=max(sub_x);
	[min_y,min_y_i]=min(sub_y);
	[max_y,max_y_i]=max(sub_y);
	edge_x=[min_x max_x max_x min_x]';
	edge_y=[min_y min_y max_y max_y]';
	if isempty(g),
		theInpLoc=[edge_x';edge_y']';
	else
		theSecLoc=g(:,ones(4,1));
		if subset_p(2)>subset_p(1),
			theInpLoc=[theSecLoc(1:subset_p(1)-1,:);edge_x';theSecLoc(subset_p(1)+1:subset_p(2)-1,:);edge_y';theSecLoc(subset_p(2)+1:size(theSec,1),:)]';
		else
			theInpLoc=[theSecLoc(1:subset_p(2)-1,:);edge_y';theSecLoc(subset_p(2)+1:subset_p(1)-1,:);edge_x';theSecLoc(subset_p(1)+1:size(theSec,1),:)]';
		end
	end
	
	try,
		z=[theInpLoc ones(4,1)]*m.linears(:,:,i)';
		patch(edge_x,edge_y,z,z);
	catch
	end
end
hold off;

function PlotData3(m, subset_p, r, g, data);

hold on;
plot3(data(:,subset_p(1))',data(:,subset_p(2))',data(:,subset_p(3)+get(m, 'n_in'))','y+');
hold off;


function PlotAntec3(m, subset_p, r, g);
hold on;
numrules=m.n_rules;
x(1,:)=cos([0:0.1:2*pi+0.1]);
x(2,:)=sin([0:0.1:2*pi+0.1]);

for i=1:numrules,
  theInv=projclus(m.ivariances(:,:,i),[subset_p(1);subset_p(2)]);
  [theRot,theHom]=eig(theInv);
  theHom=inv(sqrt(theHom));
  theRot=theRot*det(theRot);
  y=m.centers([subset_p(1) subset_p(2)],i)*ones(1,size(x,2))+2*theRot*theHom*x;
  z=r(5);
  plot3(y(1,:),y(2,:),z(1,ones(1,size(x,2))),'r');
end
hold off;
