function f=plot(m,varargin);
%PLOT performs a plot of the a mamdani 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 mamdani 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        


%   Copyright (c) 1998 by Antoine Duchateau			
%	$Revision: 0.1$  $Date: 13/09/98 12:02
%
%	All rights reserved

%First call the Parents method

check(m);
[f,m.famimo]=plot(m.famimo,varargin{:});	%the current axes should be theAxes

if nargin>=2,
	data=varargin{1};
end
if nargin>=3,
	options=varargin{2};
end

if (nargin<3)|(isempty(options))|(~isfield(options,'axes'))|isempty(options.axes),
	theAxes=gca;
	theParams=get(theAxes,'UserData');
	theParams=initialize(theAxes,theParams);
else
	theAxes=options.axes;
	theParams=get(theAxes,'UserData');	
end

theParams=adjust_controls(theAxes,theParams,m);
[theParams,m]=draw(theAxes,theParams,m);

if theParams.parentFig == 0,
	s=get(f,'UserData');
	s.model=m;
	set(f,'UserData',s);
else
	s=get(theParams.parentFig,'UserData');
	s.model=m;
	set(theParams.parentFig,'UserData',s);
end

set(theAxes,'UserData',theParams);
	
		

function theParams=initialize(theAxes,theParams);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	
%	
%	Inputs:
%			theAxes: reference to the graph axes
%			theParams: the parameters
%			
%
%	Outputs:
%			theParams: the modified parameters
%			
%			
%	Created 30/12/97 at 22:25
%	
%	by Antoine Duchateau
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


a=get(theAxes,'Parent');
b = uicontrol('Parent',a, ...
	'Units','normal', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','redraw_fig_hook', ...
	'Enable','Off', ...
	'Position',[0.25 0.10 0.11 0.06], ...
	'String','Antec.', ...
	'Style','checkbox', ...
	'TooltipString', 'Draw the membership functions', ...
	'Tag','Checkbox1', ...
	'Value',1);

theParams.CheckAntec=b;
theParams.toDelete=[theParams.toDelete b];

b = uicontrol('Parent',a, ...
	'Units','normal', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','redraw_fig_hook', ...
	'Enable','Off', ...
	'Position',[0.37 0.10 0.11 0.06], ...
	'String','Clusters', ...
	'Style','checkbox', ...
	'TooltipString', 'Draw the clusters', ...
	'Tag','Checkbox1', ...
	'Value',1);

theParams.CheckClusters=b;
theParams.toDelete=[theParams.toDelete b];


b = uicontrol('Parent',a, ...
	'Units','normal', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','redraw_fig_hook', ...
	'Enable','Off', ...
	'Position',[0.25 0.02 0.11 0.06], ...
	'String','Conseq.', ...
	'TooltipString', 'Draw the linear consequents of the rules', ...
	'Style','checkbox', ...
	'Tag','Checkbox1', ...
	'Value',1);

theParams.CheckConseq=b;
theParams.toDelete=[theParams.toDelete b];


function theParams=adjust_controls(theAxes,theParams,m,fVal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	
%	
%	Inputs:
%			theAxes: reference to the graph axes
%			theParams: the parameters
%			m: the TS model
%
%	Outputs:
%			theParams: the modified parameters
%			
%			
%	Created 30/12/97 at 22:25
%	
%	by Antoine Duchateau
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[in,out]=get_regressors(m);
curO=get(theParams.PopO1,'value');
theDim=size(in{curO},2);

%Antec Check box
if get(theParams.TypePanel(1),'value'),
	set(theParams.CheckAntec,'Enable','On');
else
	set(theParams.CheckAntec,'Enable','Off');
end

%Cluster Check box
if get(theParams.TypePanel(1),'value')&isfield(m.description{curO}.opt,'ClusVars')&...
		(size(m.description{curO}.opt.ClusVars,3)==m.description{curO}.n_rules),
	set(theParams.CheckClusters,'Enable','On');
else
	set(theParams.CheckClusters,'Enable','Off');
	set(theParams.CheckClusters,'Value',0);
end

%Model and Conseq Check box
if get(theParams.TypePanel(1),'value')&~isempty(m.description{curO}.linears),
	set(theParams.CheckConseq,'Enable','On');
	set(theParams.CheckModel,'Enable','On');
else
	set(theParams.CheckConseq,'Enable','Off');
	set(theParams.CheckModel,'Enable','Off');
end

%Axes color
set(theAxes,'position',get(theParams.SecAxes,'UserData'));
set(theParams.SecAxes,'position',get(theParams.SecAxes,'UserData'));
	
if get(theParams.TypePanel(1),'value')&(get(theParams.PopD,'value')<2),
	set(theAxes,'Color','none');
	set(theParams.SecAxes,'Visible','on');
else
	set(theAxes,'Color',[1 1 1]);
	set(theParams.SecAxes,'Visible','off');
end

set(theParams.TerAxes,'Visible','off','position',[0 0 1 1]);


function [theParams,m]=draw(theAxes,theParams,m);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	
%	
%	Inputs:
%			
%			
%			
%
%	Outputs:
%			
%			
%			
%	Created 31/12/97 at 15:29
%	
%	by Antoine Duchateau
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[in,out]=get_regressors(m);

saveView=0;

theClusters=get(theParams.CheckClusters,'value');
theAntec=get(theParams.CheckAntec,'value');
theConseq=get(theParams.CheckConseq,'value');
theError=get(theParams.CheckError,'value');
theChoiceMod=get(theParams.CheckModel,'value');

if get(theParams.TypePanel(1),'value'),
	theTypGraph=get(theParams.PopD,'value');
else
	theTypGraph=0;
end
theI1=get(theParams.PopI1,'value');
theI2=get(theParams.PopI2,'value');
theO1=get(theParams.PopO1,'value');

g=[];
if ~strcmp(get(theParams.Section,'String'),'Unused'),
	try,
		g=eval(get(theParams.Section,'String'));
	catch
		g=[];
	end
	if isempty(g)|(size(g,1)~=size(in{theO1},2)),
		[in,out]=get_regressors(m);
		g=mean(in{theO1})';
		signif=round(g)+(round(g)==0);
		g=round(g./signif*100).*signif/100;
		set(theParams.Section,'String',mat2str(g));
	end;
end
	
theDataMerged=[in{theO1} out{theO1}];
dimension=size(in{theO1},2);

if theTypGraph==0,
	axes(theAxes);
	[simul, m]=simulate(m);
	
	plot(simul,'g');
	plot(get(m,'data_out'),'y');
	plot(get(m,'data_in'));
	ax=[];
	if ~isempty(get(theParams.Info,'String')),
		try,
			ax=eval(get(theParams.Info,'String'));
		catch
			ax=[];
		end
	end
	if isempty(ax)|(size(ax,2)~=4),
		bounds=[-min(min(simul));max(max(simul))];
		range=sum(bounds);
		ax=[0 size(simul,1) -bounds(1)-0.1*range bounds(2)+0.1*range];
		signif=round(ax)+(round(ax)==0);
		ax=round(ax./signif*100).*signif/100;
		set(theParams.Info,'String',mat2str(ax,3));
	end;
	xlabel('Time');
	ylabel('Outputs');
	axis(ax);
	view(0,90);
	rotate3d off;

elseif theTypGraph==1,
	
	theAx1=[-min(theDataMerged(:,theI1)) -min(theDataMerged(:,dimension+1)); ...
			max(theDataMerged(:,theI1)) max(theDataMerged(:,dimension+1))	];
	theRanges=sum(theAx1);
	theMinMax=[-theAx1(1,:)-0.1*theRanges;theAx1(2,:)+0.1*theRanges];
	
	ax=[];
	if ~isempty(get(theParams.Info,'String')),
		try,
			ax=eval(get(theParams.Info,'String'));
		catch
			ax=[];
		end
	end
	if isempty(ax)|(size(ax,2)~=(get(theParams.PopD,'value')+1)*2),
		ax=theMinMax(:)';
		signif=round(ax)+(round(ax)==0);
		ax=round(ax./signif*100).*signif/100;
		set(theParams.Info,'String',mat2str(ax,3));
	end;
	
	
	try,
		if theClusters, PlotClus(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theAntec, PlotAntec(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theConseq, PlotConseq(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theError,
			PlotError(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g);
			set(theParams.TerAxes,'Visible','on');
			ylabel('Error');
			o_p=get(theParams.SecAxes,'UserData');
			set(theAxes,'position',o_p+[0 o_p(4)*4/11 0 -o_p(4)*4/11]);
			set(theParams.SecAxes,'position',o_p+[0 o_p(4)*4/11 0 -o_p(4)*4/11]);
			set(theParams.TerAxes,'position',o_p+[0 0 0 -o_p(4)*8/11]);
			xlabel(['Input ' num2str(theI1)]);
			ax_err=axis;
			axis([ax(1:2) ax_err(3:4)]);
		else
			o_p=get(theParams.SecAxes,'UserData');
			set(theAxes,'position',o_p);
			set(theParams.SecAxes,'position',o_p);
			set(theParams.TerAxes,'position',o_p);
		end	
	catch
		o_p=get(theParams.SecAxes,'UserData');
		set(theAxes,'position',o_p);
		set(theParams.SecAxes,'position',o_p);
		set(theParams.TerAxes,'position',o_p);
	end
	try,
		if theChoiceMod, PlotModel(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end

	axes(theParams.SecAxes);
	ylabel('Membership');
	axis([ax(1:2) 0 2]);
	
	axes(theAxes);
	
else
	
	theAx1=[-min(theDataMerged(:,theI1)) -min(theDataMerged(:,theI2)) ...
			-min(theDataMerged(:,dimension+1));max(theDataMerged(:,theI1)) ...
			max(theDataMerged(:,theI2)) max(theDataMerged(:,dimension+1))	];
			
	theRanges=sum(theAx1);
	theMinMax=[-theAx1(1,:)-0.02*theRanges;theAx1(2,:)+0.02*theRanges];
	
	
	ax=[];
	if ~isempty(get(theParams.Info,'String')),
		try,
			ax=eval(get(theParams.Info,'String'));
		catch
			ax=[];
		end
	end
	if isempty(ax)|(size(ax,2)~=(get(theParams.PopD,'value')+1)*2),
		ax=theMinMax(:)';
		signif=round(ax)+(round(ax)==0);
		ax=round(ax./signif*100).*signif/100;
		set(theParams.Info,'String',mat2str(ax,3));
	end;
	
	try,
		if theError&(theI1~=theI2), PlotError3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theClusters&(theI1~=theI2), PlotClus3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theAntec&(theI1~=theI2), PlotAntec3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theConseq&(theI1~=theI2), PlotConseq3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	try,
		if theChoiceMod&(theI1~=theI2), PlotModel3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,ax,g); end
	catch
	end
	
		
end





function PlotClus(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
if ~isempty(m.description{theO1}.opt),
	axes(theAxes);
	numrules=m.description{theO1}.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(inv(m.description{theO1}.opt.ClusVars(:,:,i)),[theI1;dimension+1]);
		[theRot,theHom]=eig(theInv);
		theHom=inv(sqrt(theHom));
		theRot=theRot*det(theRot);
		y=m.description{theO1}.opt.ClusCenters([theI1 dimension+1],i)*ones(1,size(x,2))+2*theRot*theHom*x;
		plot(y(1,:),y(2,:),'r');
	end
end



function PlotAntec(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theParams.SecAxes);
numrules=m.description{theO1}.n_rules;
hold on;
x=((0:0.01:1)*(theRange(2)-theRange(1))+theRange(1));
if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:theI1-1,:);x;theSec(theI1+1:size(theSec,1),:)];
end
y=membership(m,theInp',theO1)';
sumY=sum(y,1)+realmin;
y=y./sumY(ones(1,size(y,1)),:);
for i=1:numrules,
	plot(x,y(i,:),'b');
end
hold off;


function PlotModel(theAxes,theParams, m, theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
x=((0:0.01:1)*(theRange(2)-theRange(1))+theRange(1));
if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:theI1-1,:);x;theSec(theI1+1:size(theSec,1),:)];
end
y=eval(m,theInp',theO1)';
plot(x,y,'g');


function PlotError(theAxes,theParams, m, theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theParams.TerAxes);
x=theDataMerged(:,1:end-1);
[garbage,indices]=sort(theDataMerged(:,theI1));
in_range=(garbage>theRange(1)) & (garbage<theRange(2));
indices=indices(find(in_range));
y=eval(m,x(indices,:),theO1);
plot(x(indices,theI1),abs(y-theDataMerged(indices,end)),'k');

function PlotConseq(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
numrules=m.description{theO1}.n_rules;
x=((0:0.01:1)*(theRange(2)-theRange(1))+theRange(1));
if isempty(g),
	theInp=x;
else
	theSec=g(:,ones(1,size(x,2)));
	theInp=[theSec(1:theI1-1,:);x;theSec(theI1+1:size(theSec,1),:)];
end
memb=membership(m,theInp',theO1)';
for i=1:numrules,
	[maxis,ids]=max(memb);
	sub_x=x(:,find(ids == i))';
	sub_inp=theInp(:,find(ids == i))';
	try,
		y=[sub_inp ones(size(sub_x,1),1)]*m.description{theO1}.linears(:,:,i)';
		plot(sub_x,y,'r.');
	catch
	end
end


%function for plotting the model error in 3D plots
function PlotError3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
%get the data
data=theDataMerged(:,1:end-1);
%evaluate the model using the data and compare it with the real value (and add an offset) 
z=eval(m,data,theO1)-theDataMerged(:,end)+theRange(5);
%create the interpolation mesh
theI1_lin= linspace(theRange(1),theRange(2));
theI2_lin= linspace(theRange(3),theRange(4));
[X Y]=meshgrid(theI1_lin,theI2_lin);
%interpolate the data
error=griddata(data(:,theI1),data(:,theI2),z,X,Y,'cubic');
%plot the mesh
theSurf=surf(X, Y, error);
set(theSurf,'EdgeColor','none');

%plot the sparse error points
plot3(data(:,theI1),data(:,theI2),z,'.','MarkerSize',15);

%function for plotting the clusters in 3D
function PlotClus3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
numrules=m.description{theO1}.n_rules;
[theX,theY,theZ]=sphere(50);
theJoin=[theX(:)';theY(:)';theZ(:)'];

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

for i=1:numrules,
	theInv=projclus(inv(m.description{theO1}.opt.ClusVars(:,:,i)),[theI1;theI2;dimension+1]);
	[theRot,theHom]=eig(theInv);
	theHom=inv(sqrt(theHom));
	theRot=theRot*det(theRot);
	theJoinTrans=m.description{theO1}.opt.ClusCenters([theI1 theI2 dimension+1],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]);


function PlotModel3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
[x,y]=meshgrid(theRange(1):(theRange(2)-theRange(1))/40:theRange(2),...
				theRange(3):(theRange(4)-theRange(3))/40:theRange(4));
if isempty(g),
	theInp=[x(:)';y(:)'];
else
	theSec=g(:,ones(1,size(x,1)*size(x,2)));
	if theI2>theI1,
		theInp=[theSec(1:theI1-1,:);x(:)';theSec(theI1+1:theI2-1,:);y(:)';theSec(theI2+1:size(theSec,1),:)];
	else
		theInp=[theSec(1:theI2-1,:);y(:)';theSec(theI2+1:theI1-1,:);x(:)';theSec(theI1+1:size(theSec,1),:)];
	end
end;
z=eval(m,theInp',theO1)';
theSurf=surf(x,y,reshape(z,size(x)));
set(theSurf,'EdgeColor','none');

colormap(hsv);

function PlotConseq3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
numrules=m.description{theO1}.n_rules;
[x,y]=meshgrid(theRange(1):(theRange(2)-theRange(1))/40:theRange(2),...
				theRange(3):(theRange(4)-theRange(3))/40:theRange(4));

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

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)|(dimension<=2),
		theInpLoc=[edge_x';edge_y']';
	else
		theSecLoc=g(:,ones(4,1));
		if theI2>theI1,
			theInpLoc=[theSecLoc(1:theI1-1,:);edge_x';theSecLoc(theI1+1:theI2-1,:);edge_y';theSecLoc(theI2+1:size(theSec,1),:)]';
		else
			theInpLoc=[theSecLoc(1:theI2-1,:);edge_y';theSecLoc(theI2+1:theI1-1,:);edge_x';theSecLoc(theI1+1:size(theSec,1),:)]';
		end
	end
	
	try,
		z=[theInpLoc ones(4,1)]*m.description{theO1}.linears(:,:,i)';
		patch(edge_x,edge_y,z,z);
	catch
	end
end


	


function PlotAntec3(theAxes,theParams,m,theDataMerged,dimension,theI1,theI2,theO1,theRange,g);
axes(theAxes);
numrules=m.description{theO1}.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.description{theO1}.ivariances(:,:,i),[theI1;theI2]);
  [theRot,theHom]=eig(theInv);
  theHom=inv(sqrt(theHom));
  theRot=theRot*det(theRot);
  y=m.description{theO1}.centers([theI1 theI2],i)*ones(1,size(x,2))+2*theRot*theHom*x;
  z=theRange(5);
  plot3(y(1,:),y(2,:),z(1,ones(1,size(x,2))),'r');
end

