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. The accepted fields for options
%    are:
%    
%      METHOD: identification method(see below).
%      N_RULES: number of rules.
%      MIN_N_RULES: min. nb. of rules (incremental methods
%          only)
%      MAX_N_RULES: max. nb. of rules (incremental methods
%          only)
%      FUZZYNESS: fuzzyness of the model (usually greater
%          than 1)
%      TOLERANCE: minimal improvement for recursive
%          techniques before the algorithm stops.
%      SEED: seed for random generator
%      LIN_FIT: type of linear fitting: 0: normal, 1: local
%      RULE_TYPE: rule type code (see ADD_RULES for more
%          info)
%    
%    Methods of identification are:
%    
%      CLUSLMS: GK clustering followed by a LMS fit
%      FMCLUST: Robert Babuska's identification method
%      INCRSIE: Incremental identification method (Siemens)
%      TSGKLMXV: incremental method developped in IRIDIA
%      CLUSLEV: GK clust. init. and Lev.-Marq. optimisation
%      RANDLM: random init. and Lev.-Marq. optimisation
%    
%    

%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>.
%
%	The accepted fields for options are:
%<dl>
%<dt>	method   <dd>      	identification method(see below).
%<dt>	n_rules   <dd>      number of rules.
%<dt>   min_n_rules <dd>    min. nb. of rules (incremental methods only)
%<dt>   max_n_rules <dd>    max. nb. of rules (incremental methods only)
%<dt>   fuzzyness  <dd>     fuzzyness of the model (usually greater than 1)
%<dt>	tolerance <dd>      minimal improvement for recursive techniques before the algorithm stops.
%<dt>   seed     <dd>       seed for random generator
%<dt>	lin_fit	<dd>        type of linear fitting: 0: normal, 1: local
%<dt>   rule_type  <dd>     rule type code (see <CODE>add_rules</CODE> for more info)
%</dl>	
%
%	Methods of identification are:
%<dl>
%<dt>  cluslms <dd> GK clustering followed by a LMS fit
%<dt>  fmclust <dd> Robert Babuska's identification method
%<dt>  incrsie <dd> Incremental identification method (Siemens)
%<dt>  tsgklmxv  <dd> incremental method developped in IRIDIA
%<dt>  cluslev <dd> GK clust. init. and Lev.-Marq. optimisation
%<dt>  randlm  <dd> random init. and Lev.-Marq. optimisation
%</dl>
%
%MAN_PAGE_END   

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);
		m.linears=zeros(get(m,'n_out'),get(m,'n_in')+1,get(m,'n_rules'));
	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);
		options.update = 1;
		m=clustering(m,in,out,options);
		%We need to differentiate the normalised (where no variance optimisation is necessary) and the non normalised case
		if (strcmp(m.model_code{2}, 'ungaussian'))| ...
		      (strcmp(m.model_code{2}, 'ungaussian')),
			globVariance = 1;
			%Try bigger
			policy = 2;
			ttry = 0;
			
			locerror = error(fit_linears(m,in,out,options.lin_fit),in,out);
			while (abs(log(policy))>0.05),
				if (ttry == 0),
					ttry = 1;
				elseif (ttry == 1),
					policy = 1/policy;
					ttry =2;
				elseif (ttry == 2)
					policy = sqrt(1/policy);
					ttry = 0;
				end
				
				if ~((globVariance >= 8) & (policy>1)),
					locVariance = globVariance*policy;
					tm = m;
					tm.ivariances = tm.ivariances/locVariance;
					tm = fit_linears(tm,in,out,options.lin_fit);
				
					locerrorbis = error(tm,in,out);
				
					if locerrorbis < locerror,
						ttry = 0;
						locerror = locerrorbis;
						globVariance = locVariance;
					end
				end
			end
			m.ivariances = m.ivariances/globVariance;
			m = fit_linears(m,in,out,options.lin_fit);
		else
			m = fit_linears(m,in,out,options.lin_fit);
		end
	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,in,out);
		if isfield(options,'tolerance'),
			m=lev_marq(m,in,out,[],options.tolerance);
		else
			m=lev_marq(m,in,out,options.lin_fit);
		end
	else
		error(['Sorry, this identification method only works for single output systems']);
	end
else
	error([action ' is an unknown method, see help file']);
end

if ~all(all(isfinite(get(m,'limits')))),
	set(m,'limits',[min(in) min(out);max(in) max(out)]);
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,'type'),
	options.type = {'productspace' 'inversedist' 'linear'};
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(['Progress: updating ' num2str(initN) ' clusters ...\n']);
  	[thePNU,thePNCenters,thePNVars]=gkfast([in out],old_c,options.fuzzyness);
else
	fprintf(['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,2)+1 size(in,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);
m=add_rules(m,size(theClusVars,3),options.type);

opt = get(m,'opt');
opt.PNCenters=thePNCenters;
opt.PNVars=thePNVars;
opt.PNU=thePNU;

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

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

opt.IDs=theIDs;

set(m,'opt',opt);

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

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,'type'),
	options.type = {'productspace' 'inversedist' 'linear'};
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);
m=add_rules(m,options.theNumber,options.type);

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

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

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

m=denormalise(m);

