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

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

p = get(m,'params');

%allocate space
ndp = size(in,1);
num = 0;
op = get(m,'optimparams');
test_m=m;

while((new_err<old_err*(1-treshold))|(num<10))&(lambda<1e20),
	old_err = new_err;
	num = num+1;
	
	theA = zeros(length(p));
	theB = zeros(1,length(p));
	
	for i = 1:size(in,1),
		%compute the derivatives
		if strcmpi(m.model_code{2},'gaussian'),
			[dTemp,y]=jacob_params(in(i,:),m.centers,m.ivariances,m.linears,1,1,1);
		elseif strcmpi(m.model_code{2},'inversedist'),
			[dTemp,y]=jacob_params(in(i,:),m.centers,m.ivariances,m.linears,1,2,m.m);
		elseif strcmpi(m.model_code{2},'ungaussian'),
			[dTemp,y]=jacob_params(in(i,:),m.centers,m.ivariances,m.linears,0,1,1);
		elseif strcmpi(m.model_code{2},'uninversedist'),
			[dTemp,y]=jacob_params(in(i,:),m.centers,m.ivariances,m.linears,0,2,m.m);
		end
		dd = dTemp(op,:)';
		theA = theA + dd(:,op)' * cost * dd(:,op);
		theB = theB + (out(i,:) - y) * cost * dd(:,op);
	end
	
	theB = theB';
	
	new_err = 2*old_err;
	
	while (new_err>old_err)&(lambda<2e21),
		lambda=lambda*10;
		test_p = p + ((theA+diag(lambda(ones(length(p),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;
	m = set(m,'params',p);
	
end


