function m=identify(m,in,out,options)
%IDENTIFY identifies a fuzzy model from data
%    
%    M=IDENTIFY(M,IN,OUT,OPTIONS)  identify model M on the
%    basis of the data inside IN and OUT using method
%    specified in OPTIONS.
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose   identifies a fuzzy model from data.
%
%   @synopsis m=identify(m,in,out,options) 
%   @description  identify model <CODE>m</CODE> on the basis of the data
%   inside <CODE>in</CODE> and <CODE>out</CODE> using method specified in <CODE>options</CODE>.
%
%MAN_PAGE_END   

error(['Sorry, not yet operational']);

if nargin<4,
	options = struct('method','cluslms');
end

error(checkargs({m in out options},4,'taksug','numeric','numeric','struct'));

check(m);

m=rem_rules(m);

if nargin<3, options=[]; end
if ~isfield(options,'lin_fit'), options.lin_fit=[]; end
if ~isfield(options,'method'), options.method='fmclust'; end

action = options.method;

fprintf(['\n Using method ' action ' for identifying the system \n']);

if strcmp(action,'clus'),
	if get(m,'n_out') == 1,
		m=clustering(m,in,out,options);
		set(m,{1:get(m,'n_out') 'linears'},[]);
	else
		error(['Sorry, this identification method only works for single output systems']);
	end
	
elseif strcmp(action,'cluslms'),
	if get(m,'n_out') == 1,
		m=clustering(m,in,out,options);
		m=fit_linears(m,in,out,[],options.lin_fit);
	else
		error(['Sorry, this identification method only works for single output systems']);
	end
	
elseif strcmp(action,'fmclust'),
	if get(m,'n_out') == 1,
		m=fmclust_wraper(m,in,out,options);
	else
		error(['Sorry, this identification method only works for single output systems']);
	end

elseif strcmp(action,'incrsie'),
	m=incr_test(m,in,out,options);

elseif strcmp(action,'cluslev'),
	if get(m,'n_out') == 1,
		m=clustering(m,in,out,options);
		m=fit_linears(m,in,out,options.lin_fit);
		if isfield(options,'tolerance'),
			m=lev_marq(m,in,out,[],options.tolerance);
		else
			m=lev_marq(m,in,out);
		end
	else
		error(['Sorry, this identification method only works for single output systems']);
	end

	
elseif strcmp(action,'tsgklmsxv'),
	if get(m,'n_out') == 1,
		m=clustering(m,in,out,options);
		m=lms_incremental(m,in,out,options);
		m=clustering(m,in,out,struct('update',1));
		m=fit_linears(m,in,out,[],options.lin_fit);
	else
		error(['Sorry, this identification method only works for single output systems']);
	end

	
elseif strcmp(action,'randlm'),
	if get(m,'n_out') == 1,
		m=randomInit(m,in,out,options);
		m=fit_linears(m);
		if isfield(options,'tolerance'),
			m=lev_marq(m,in,out,[],options.tolerance);
		else
			m=lev_marq(m,in,out);
		end
	else
		error(['Sorry, this identification method only works for single output systems']);
	end

end


%Performs a hierarchical or flat clustering

function m=clustering(m,in,out,options);

if isfield(options,'n_rules'),
	initN=options.n_rules;
	targetN=options.n_rules;
elseif isfield(options,'min_n_rules') & isfield(options,'max_n_rules')
	initN=options.max_n_rules;
	targetN=options.min_n_rules;
else
	initN=5;
	targetN=5;
end

if ~isfield(options,'update'),
	options.update=0;
end

if ~isfield(options,'fuzzyness'),
	options.fuzzyness = 2;
end

if options.update,
	initN = m.n_rules;
	targetN=initN;
	old_c = m.centers';
	old_c = [old_c eval(m,old_c)];

	fprintf(['    Output ' num2str(output), ', progress: updating ' num2str(initN) ' clusters ...\n']);
  		[thePNU,thePNCenters,thePNVars]=gkfast([in out],old_c,options.fuzzyness);
else
	fprintf(['    Output ' num2str(output), ', progress: clustering data in ' num2str(initN) ' clusters ...\n']);
  		[thePNU,thePNCenters,thePNVars]=gkfast([in out],initN,options.fuzzyness);
  	end

	thePNU=thePNU';
	thePNCenters=thePNCenters';
	thePNVars=reshape(thePNVars,[size(in{output},2)+1 size(in{output},2)+1 initN]);

    if targetN<initN,
		fprintf(1,['                    agregating clusters ...\n']);
		[theHierarchy,thePSVars,thePSCenters]=gkreduc2(thePNCenters,thePNVars,thePNU,targetN,2,[in out]');
		[theClusVars,theClusCenters,theIDs]=get1lev(theHierarchy,thePSVars,thePSCenters,thePNVars,thePNCenters);
	else
		theHierarchy=[];
		thePSVars=[];
		thePSCenters=[];
		theClusVars=thePNVars;
		theClusCenters=thePNCenters;
		theIDs=[];
	end

theVars=[];

for j=1:size(theClusVars,3),
	theVars(:,:,j)=projclus(inv(theClusVars(:,:,j)),(1:(size(theClusVars,1)-1))');
end

theCenters=theClusCenters(1:(size(theClusCenters,1)-1),:);

m=rem_rules(m,output);
m=add_rules(m,size(theClusVars,3),{'productspace' 'inversedist' 'linear'});

m.opt.PNCenters=thePNCenters;
m.opt.PNVars=thePNVars;
m.opt.PNU=thePNU;

m.opt.Hierarchy=theHierarchy;
m.opt.PSVars=thePSVars;
m.opt.PSCenters=thePSCenters;

m.opt.ClusVars=theClusVars;
m.opt.ClusCenters=theClusCenters;

m.ivariances=theVars;
m.centers=theCenters;
m.m=options.fuzzyness;

m.opt.IDs=theIDs;

m.linears=zeros(1,size(theCenters,1)+1,size(theClusVars,3));



% Randomly initialise the rules
function m=randomInit(m,in,out,options);

if nargin<4, options=[]; end

if ~isfield(options,'theNumber')|isempty(options.theNumber),
	options.theNumber = 5;
end
if ~isfield(options,'theWidth')|isempty(options.theWidth),
	options.theWidth = 2;
end

dimension=get(m,'n_in');

m=rem_rules(m,i);
m=add_rules(m,options.theNumber,{'productspace' 'inversedist' 'linear'});

for j=1:options.theNumber,
	m.ivariances(:,:,j)=eye(dimension)*10/options.theWidth;
end;

m.opt=[];
m.centers=rand(dimension,options.theNumber);
m.linears=zeros(1,dimension+1,options.theNumber);

m.opt.original_limits = limits;
limits(1,:) = 0;
limits(2,:) = 1;

m=denormalise(m);

