function out=jacob_inputs(m,varargin)
%JACOB_INPUTS computes the jacobian of the model
%    
%    The procedure for computing the jacobian depends on the
%    underlying representation  of the model. The result OUT
%    is an array of size number of outputs * number of inputs
%    containing the derivatives of the outputs with respect to
%    the inputs.
%    
%    JACOB_INPUTS(M,X) computes for each output the jacobian
%    of the model with respect to the x (the input of the
%    linear) at point X.
%    
%    JACOB_INPUTS(M,X,J) computes for output J the jacobian of
%    the model with respect to the x (the input of the linear)
%    at point X.
%    
%    

%MAN_PAGE_BEGIN
%
%   @purpose    computes the jacobian of the model.
%   The procedure for computing the jacobian depends on the underlying representation 
%   of the model. The result <code>out</code> is an array of size number of outputs * number of inputs
%	containing the derivatives of the outputs with respect to the inputs.
%
%   @synopsis   jacob_inputs(m,x)
%   @description computes for each output the jacobian of
%  the model with respect to the x (the input of the
%  linear) at point <CODE>x</CODE>.
%
%   @synopsis   jacob_inputs(m,x,j)
%   @description computes for output J the jacobian of
%  the model with respect to the x (the input of the
%  linear) at point <CODE>x</CODE>.
%
%MAN_PAGE_END

%   Copyright (c) 1998 by Antoine Duchateau			
%	$Revision: 0.1$  $Date: 10/12/1997
%
%	All rights reserved

error(checkargs([{m} varargin],[2 3],'taksug','numeric','numeric'));

check(m);

x = varargin{1};
if (nargin == 2)|isempty(varargin{2}),
	j = 1:get(m,'n_out');
else
	j = varargin{2};
end

%Can't use the mex file so do it by hand
nr = get(m,'n_rules');
ni = get(m,'n_in');
no = get(m,'n_out');

idxs = find(eye(ni) == 1);	% Which are the indices of a diagonal in a matrix ?

out = zeros(no,ni,size(x,1));

% Compute the degrees of fullfilment of the rules
for i=1:ni,
	mg(:,:,i) = mgrade(x(:,i), m.mfs{i}(m.rls(:,i),:));
	dmg(:,:,i) = dmgrade(x(:,i), m.mfs{i}(m.rls(:,i),:));
end

dof = prod(mg,3);
sdof = sum(dof,2) + realmin;
ndof = dof ./ sdof(:,ones(nr,1));
	
	
% This is awfull, I have to hack a faster version of this in the future
for nex=1:size(x,1),
	for i=1:nr,
		fsets = squeeze(mg(nex,zeros(1,ni)+i,:))';
		fsets(idxs) = dmg(nex,i,:);
		dAidx(i,:) = prod(fsets,1);
	end
	
	lpart = zeros(no,ni);
	oval = zeros(no,1);
	
	for i=1:nr,
		lpart = lpart + m.linears(:,:,i) * [x(nex,:)';1] * dAidx(i,:) + dof(nex,i) * m.linears(:,1:end-1,i);
		oval = oval + dof(nex,i) * m.linears(:,:,i) * [x(nex,:)';1];
	end
	
	out(:,:,nex) = (lpart - oval * sum(dAidx,1) / sdof(nex)) / sdof(nex);
end

