function m=lev_marq(m,select,treshold,do_plot);
%LEV_MARQ performs a fitting of the model
%    
%    M=LEV_MARQ(M,SELECT,THRESHOLD,DO_PLOT) performs a Levenbergh Marquardt
%    optimisation of the model M. A subset of the parameters defined by SELECT
%    is taken into account. SELECT is a binary number wxyz. w is set for 
%    the parameters related to the centers. x for the covariance matrices,
%    y for the linears consequents and z for the affine terms. THRESHOLD is the
%    minimum difference between two consecutive normalised mean square error
%    before the optimisation is stopped. If DO_PLOT is greater than  0, an
%    intermediary plot is drawed every DO_PLOT steps.
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose performs a fitting of the model.
%   
%   @synopsis m=lev_marq(m,select,threshold,do_plot)
%   @description  performs a Levenbergh Marquardt
%  optimisation of the model <code>m</code>. A subset of the parameters defined by <code>select</code>
%  is taken into account. <code>select</code> is a binary number wxyz. w is set for 
%  the parameters related to the centers. x for the covariance matrices,
%  y for the linears consequents and z for the affine terms. <code>threshold</code> is the
%  minimum difference between two consecutive normalised mean square error
%  before the optimisation is stopped. If <code>do_plot</code> is greater than  0, an
%  intermediary plot is drawed every <code>do_plot</code> steps.
%
%MAN_PAGE_END



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

check(m);
for no=1:get(m,'n_out'),
	if strcmpi(m.description{no}.model_code{1},'projected'),
		error('This function only supports product space antecedents');
	end
end

if (nargin < 2)|isempty(select), select = 3; end
if (nargin < 3)|isempty(treshold), treshold = 0.001; end
if (nargin < 4)|isempty(do_plot), do_plot = 0; end

m=normalise(m);

if do_plot,
	ax=plot(m);
	option=get(ax,'userdata');
end
	
[in,out]=get_regressors(m);

for j=1:get(m,'n_out'),

	lambda = 0.0001;
	old_err = realmax;
	new_err = error(m,j);
	p = get_params(m,select,j);
	
	%allocate space
	dfdp = zeros(size(p,1),size(in{j},1));
	y = zeros(size(in{j},1),1);
	num = 0;
	
	test_m=m;
	
	while((new_err<old_err*(1-treshold))|(num<10))&(lambda<1e9),
		old_err = new_err;
		num = num+1;
		
		%compute the derivatives
		for ii = 1:size(in{j},1),
			[dfdp(:,ii),y(ii)]=jacob_params(in{j}(ii,:),test_m.description{j}.centers, test_m.description{j}.ivariances, test_m.description{j}.linears,uint8(select));
		end
		
		indiv_errors = (y - out{j})';
		
		theA = sumproduct(dfdp);
		theB = sum((indiv_errors(ones(1,size(p,1)),:).*dfdp)')';
		
		new_err = 2*old_err;
		
		while (new_err>old_err)&(lambda<1e10),
			lambda=lambda*10;
			
			test_p = p - (theA+diag(lambda(ones(size(p,1),1))))\theB;
			test_m = set_params(m,test_p,select,j);
			
			try,
				new_err = error(test_m,j);
				if ~isfinite(new_err),
					new_err = realmax;
					test_p = p;
				end
			catch,
				new_err = realmax;
				test_p = p;
			end
		end
		lambda=lambda/20;
		p = test_p;
		
		if do_plot&(rem(num,do_plot)==do_plot-1),
			plot(test_m,[],option);
			drawnow;
		end
		
	end
	
	m = set_params(m,p,select,j);
	
end

m=denormalise(m);

if do_plot,
	ax=plot(m);
end

