function cout = amutual(cin, len)

% core/amutual
%
% amutual(cin, len)
%
% cin=core object
% len=maximal lag?
%
% mutual information of a time series against itself, with increasing lag
%
% cmerk Maerz 1998

partitionen = 64;
epsilon = 1e-10;

N = dlens(cin,1);
points = data(cin);

points = points - min(points);
points = 1+floor(points / (max(points)/(partitionen-epsilon)));
% Aus den Werten sind jetzt fertige Indizes fuer die Partititonen geworden

amf = zeros(len,1);

for lag=0:len-1
	overlap = N-lag;
	increment = 1/overlap;
	
	one = ones(overlap,1);
	histA= sparse(points(1:overlap), one, increment);
	histB= sparse(one, points(1+lag:N), increment);
	histAB = sparse(points(1:overlap), points(1+lag:N), increment);
	imagesc(histAB);

	[ind1, ind2, value] = find(histAB);
	amf(lag+1) = sum(value .* log2(value ./ (histA(ind1) .* histB(ind2)')));
end

cout = core(amf);

	
% 	s = 0;		
% 	for i=1:length(ind1)
% 		s = s + value(i) * log2(value(i)/(histA(ind1(i))*histB(ind2(i))));
% 	end
	
% 	sum(histA)  
% 	sum(histB)
% 	sum(sum(histAB))

% 	sp2 = sparse(histA * histB);
% 	index = find(sp2);
% 	quot = histAB;
% 	quot(index) = quot(index) ./  sp2(index);	% in dieser Matrix stehen die Quotienten aus P(AB)/(P(A)*P(B))
% 									
% 	index = find(quot); % Nur Elemente logarithmieren, die ungleich Null sind
% 	quot(index) = log2(quot(index));	
% 
% 	amf(lag+1) = sum(sum(histAB(index) .* quot(index)));	% Gewichtung und Aufsummierung
% 	amf(lag+1) = s;
