function m=lms_incremental(m,in,out,options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% m=lm_incremental(m,options);
% 
% Inputs:
%			
%			
%			
%
% Outputs:
%			
%			
%			
% Created 5/01/98 at 14:29
% 
% by Antoine Duchateau
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Default number of sets for the xval is 5
if ~isfield(options,'n_xval'), options.n_xval=5; end

n_out=get(m,'n_out');

t_in={};t_out={};x_in={};x_out={};

% Build the training and the validation sets
if isfield(options,'xval_data'),
	t_in{1}=in;
	t_out{1}=out;
	x_in{1} = options.xval_data.dat_in;
	x_out{1} = options.xval_data.dat_out;
else
	for i=1:options.n_xval,
		[t_in{i},t_out{i},x_in{i},x_out{i}]=splitsets(in,out,i,options.n_xval);
	end
end

theNumChildren=2;

best_m=m;
ref_m=m;
tm=m;


fprintf(1,['Cross validating architectures ...\n']);
for j=1:n_out,
	fprintf(['Output ' num2str(j), '\n']);
	
	% best_m is the best model found so far
	% ref_m is the smallest model we are trying to extend (it means that
	% if we cannot get something better than best_m, only expanding ref_m,
	% we will give up).
	% tm is the current model tested
	% The algorithm is the following:
	%	best_m=ref_m=tm=m   --Initialisation
	%   tm=new model with one split
	%   if tm is better than best_m then best_m=ref_m=tm
	%	else try another split for tm.
	%	if no split is good, keep the less bad and
	%   continue to try to expand tm.
	%   if a good (better than best_m) model is found,
	%   get rid of the unnecessary split starting from ref_m
	%   ref_m = the model obtained
	%   tm=ref_m and if ref_m is better than best_m, best_m=ref_m
	%   continue until they aren't any possible splitting anymore

	
	n_rules=m.description{j}.n_rules;
	elim=ones(n_rules,1);
	
	theNewError=0;
	for i=1:size(t_in,2),
		em=add_data(m,t_in{i},t_out{i});
		em=fit_linears(em,j,options.lin_fit);
		theNewError=theNewError+error(em,struct('dat_in',x_in{i},'dat_out',x_out{i}),j);
	end		%Compute the real xval error
	
	theBestCurrent = realmax;
	secondChoice=[];
	
	while (max(elim)>0),
		theOldError=theNewError;
		
		[r_in,r_out]=get_regressors(m);
		
		m_out = eval(m,r_in{j},j);
		theLocalErr= (r_out{j} - m_out).^2;
		
		memb=membership(m,r_in{j},j)';
		
		tt1=memb*theLocalErr;
		tt1=tt1./(sum(memb')'+realmin);
		tt1=tt1.*elim;
		[maxv,maxi]=max(tt1);
		
		elim(maxi)=0;
		
		theSplittedRule=maxi;
		theSplittedRuleID=m.description{j}.opt.IDs(maxi);
		
		if theSplittedRuleID<0,
			
			theLine=find(m.description{j}.opt.Hierarchy(:,1)==theSplittedRuleID);
			
			tm=split_rule(m,j,theSplittedRuleID,theSplittedRule,theLine,theNumChildren);
			
			
			fprintf(1,['                    best so far is ',num2str(best_m.description{j}.n_rules),' rules, now testing ',...
					num2str(tm.description{j}.n_rules),' rules...\n']);
			
			theNewError = 0;
						
			for i=1:size(t_in,2),
				em=add_data(tm,t_in{i},t_out{i});
				em=fit_linears(em,j,options.lin_fit);
				theNewError=theNewError+error(em,struct('dat_in',x_in{i},'dat_out',x_out{i}),j);
			end
			
			if theNewError<theOldError,
				
				%Get rid of unuseful splittings
				%build the list of the necessary ones
				necessary=[theSplittedRuleID];
				current=theSplittedRuleID;
				while ~isempty(current)&~any(current==ref_m.description{j}.opt.IDs),
					[row,col]=find(ref_m.description{j}.opt.Hierarchy(:,[2 3]) == current);
					current=ref_m.description{j}.opt.Hierarchy(row,1);
					necessary=[current necessary];
				end
				
				for srID=necessary,
					theLine=find(ref_m.description{j}.opt.Hierarchy(:,1)==srID);
					ref_m=split_rule(ref_m,j,srID,find(ref_m.description{j}.opt.IDs == srID),theLine,theNumChildren);
				end
				
				m=ref_m;
				
				if (size(ref_m.description{j}.opt.IDs,2)~=size(tm.description{j}.opt.IDs,2))|...
				any(ref_m.description{j}.opt.IDs ~= ref_m.description{j}.opt.IDs),
				
					theNewError = 0;
					
					for i=1:size(t_in,2),
						em=add_data(m,t_in{i},t_out{i});
						em=fit_linears(em,j,options.lin_fit);
						theNewError=theNewError+error(em,struct('dat_in',x_in{i},'dat_out',x_out{i}),j);
					end
				end
				
				if theNewError<theOldError,
					best_m=m;
				else
					theNewError=theOldError;
				end
				
				n_rules=m.description{j}.n_rules;
				theBestCurrent = realmax;
				secondChoice=[];
				secondChoiceID=0;
				elim=ones(n_rules,1);
			else
				if theNewError<theBestCurrent,
					secondChoice=tm;
					secondChoiceID=theSplittedRuleID;
					theBestCurrent=theNewError;
				end
				theNewError=theOldError;
			end
			
		end
	
		if (max(elim)==0)&~isempty(secondChoice),
			m=secondChoice;
			n_rules=n_rules+theNumChildren-1;
			theBestCurrent = realmax;
			secondChoice=[];
			secondChoiceID=0;
			elim=ones(n_rules,1);
		end
				
	end
	fprintf(1,['                    final structure is ',num2str(best_m.description{j}.n_rules),' rules !\n']);
end

m=best_m;
