function m = lev_marq(m, in, out, cost, treshold)
%LEV_MARQ Performs a Levenbergh Marquardt optimisation of
%    the mapping
%    
%    FUNCTION M = LEV_MARQ(M, IN, OUT, COST) Optimises the
%    mapping M with respect to input/outpu couples
%    <CODE>in/out</CODE>. This method is mainly used
%    internaly. Identify should be prefered.
%    
%    See also: IDENTIFY
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose	Performs a Levenbergh Marquardt optimisation of the mapping
%
%   @synopsis	function m = lev_marq(m, in, out, cost)
%   @description  
%	Optimises the mapping <CODE>m</CODE> with respect to input/outpu couples
%	<CODE>in/out</CODE>.
%	This method is mainly used internaly. Identify should be prefered.
%
%   @see identify
%
%MAN_PAGE_END

% By Antoine Duchateau. All Rights Reserved.
% Written: 20/08/99

check(m);

if strcmpi(m.model_code{1},'projected'),
	error('This function only supports product space antecedents');
end

if (nargin < 4)|isempty(cost), cost = eye(get(m,'n_out')); end
if (nargin < 5)|isempty(treshold), treshold = 0.001; end

l = get(m,'limits');

for i=1:get(m,'n_in'),
	if isinf(l(1,i)),
		l(1,i) = min(in(:,i));
	end
	if isinf(l(2,i)),
		l(1,i) = max(in(:,i));
	end
end

m=normalise(m);

lambda = 0.0001;
old_err = realmax;
errs = error(m,in,out);
new_err = errs'*cost*errs;

p = get(m,'params');

%allocate space
ndp = size(in,1);
dfdp = zeros(get(m,'n_out'),size(p,1),ndp);
y = zeros(ndp,get(m,'n_out'));
num = 0;

test_m=m;

while((new_err<old_err*(1-treshold))|(num<10))&(lambda<1e9),
	old_err = new_err;
	num = num+1;
	
	theA = zeros(size(p,1));
	theB = zeros(size(p,1),1);
	
	for j = 1:size(in,2),
		%compute the derivatives
		if strcmpi(m.model_code{2},'gaussian'),
			[dfdp(:,:,i),y(i,:)]=jacob_params(x,m.centers,m.ivariances,m.linears,1,1,uint8(15))';
		elseif strcmpi(m.model_code{2},'inversedist'),
			[dfdp(:,:,i),y(i,:)]=jacob_params(x,m.centers,m.ivariances,m.linears,2,m.m,uint8(15))';
		end
		theA = theA + dfdp(:,:,i)' * cost * dfdp(:,:,i)';
		theB = theB + (out(i,:) - y(i,:)) * cost * dfdp(:,:,i);
	end
	
	deltas = (y - out)';
	
	theB = theB';
	
	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(m,'params',test_p);
		
		try,
			new_err = error(test_m,in,out);
			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;
	
end

m = set(m,'params',p);
m=denormalise(m);
