function m=lms_incremental(m,in,out,options)
%MAN_PAGE_BEGIN
%
%   @purpose	incrementally build a TS model adding rules one by one.
%
%   @synopsis	m=lms_incremental(m,in,out,options)
%   @description  
%
%   @see 
%
%MAN_PAGE_END

% By Antoine Duchateau. All Rights Reserved.
% Written: 05/01/98

% Default number of sets for the xval is 10
if ~isfield(options,'n_xval'), options.n_xval=10; 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']);

% 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.n_rules;
elim=ones(n_rules,1);

theNewError=0;
for i=1:size(t_in,2),
	theNewError=theNewError+flat_error(fit_linears(m,t_in{i},t_out{i},options.lin_fit),x_in{i},x_out{i});
end		%Compute the real xval error

theBestCurrent = realmax;
secondChoice=[];

globVariance = 1;

while (max(elim)>0),
	theOldError=theNewError;
	
	m_out = eval(m,in);
	theLocalErr= (out - m_out).^2;
	
	memb=membership(m,in)';
	
	tt1=memb*theLocalErr;
	tt1=tt1./(sum(memb')'+realmin);
	tt1=tt1.*elim;
	[maxv,maxi]=max(tt1);
	
	elim(maxi)=0;
	
	theSplittedRule=maxi;
	theSplittedRuleID=m.mapping.opt.IDs(maxi);
	
	if theSplittedRuleID<0,
		
		theLine=find(m.mapping.opt.Hierarchy(:,1)==theSplittedRuleID);
		
		tm=split_rule(m,theSplittedRuleID,theSplittedRule,theLine,theNumChildren);
		
		
		fprintf(1,['                    best so far is ',num2str(best_m.n_rules),' rules, now testing ',...
				num2str(tm.n_rules),' rules...\n']);
		
		theNewError = 0;
		
		for i=1:size(t_in,2),
			%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')),
				%Try bigger
				policy = 2;
				ttry = 0;
				
				locerror = error(fit_linears(tm,t_in{i},t_out{i},options.lin_fit),t_in{i},t_out{i});
				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;
						ttm = tm;
						ttm.ivariances = ttm.ivariances/locVariance;
						ttm = fit_linears(ttm,in,out,options.lin_fit);
					
						locerrorbis = error(ttm,in,out);
					
						if locerrorbis < locerror,
							ttry = 0;
							locerror = locerrorbis;
							globVariance = locVariance;
						end
					end
				end
			else
				ttm = fit_linears(tm,t_in{i},t_out{i},options.lin_fit);
			end
			
			theNewError=theNewError+flat_error(ttm,x_in{i},x_out{i});
		end
		
		% clf;
		% plot(ttm,'axis',[0 10 -2 2],'data',[(0:0.1:10)' sin(0:0.1:10)'])
		% drawnow;
 		
		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.mapping.opt.IDs),
				[row,col]=find(ref_m.mapping.opt.Hierarchy(:,[2 3]) == current);
				current=ref_m.mapping.opt.Hierarchy(row,1);
				necessary=[current necessary];
			end
			
			for srID=necessary,
				theLine=find(ref_m.mapping.opt.Hierarchy(:,1)==srID);
				ref_m=split_rule(ref_m,srID,find(ref_m.mapping.opt.IDs == srID),theLine,theNumChildren);
			end
			
			m=ref_m;
			
			if (size(ref_m.mapping.opt.IDs,2)~=size(tm.mapping.opt.IDs,2))|...
			any(ref_m.mapping.opt.IDs ~= ref_m.mapping.opt.IDs),
			
				theNewError = 0;
				
				for i=1:size(t_in,2),
					theNewError=theNewError+flat_error(fit_linears(m,t_in{i},t_out{i},options.lin_fit),x_in{i},x_out{i});
				end
			end
			
			if theNewError<theOldError,
				best_m=m;
			else
				theNewError=theOldError;
			end
			
			n_rules=m.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.n_rules),' rules !\n']);

m=best_m;
