function [abetti,specGaps,revals] = apBetti0(X,filter,levels),
% Compute approximate betti numbers, or spectral gaps over density filter
%	[abetti,specGaps,revals] = apBetti0(X,filter,levels),
%

disp('Compute persistent approximate betti numbers via spectral gaps');
tic

%load T5000.mat X

%filter = full(sum(X,2));
[ig, rk] = sort(filter,'descend');
ig = ig/max(filter);

MAXNUMBEREIG = 20;
k = MAXNUMBEREIG; % number of eigenvalues
eps = 1e-4; % perturb the matrix X+eps*I, avoiding the ill-conditioning
opt.disp = 0;
randSurf = 1e-1;
N = length(filter);

disp(sprintf('Level\tSize\t#Comp\tB0_1\tGap_1\t\tB0_2\tGap_2\t\tB0_3\tGap_3'));
for i=1:length(levels),
	if i>1,
		m = levels(i-1)+1;
	else
		m = 1;
	end
	
	if levels(i)>1,
		n = levels(i);
	else
		n = sum(ig>=levels(i));
	end

	dispInfo = sprintf('%d\t%d',i,n);
	if n==1,
		Components{i} = 1;
		specGaps{i} = ones(1,MAXNUMBEREIG);
		aBettis{i} = [1,zeros(1,MAXNUMBEREIG-1)];
	else
		Y = X(rk(m:n),rk(m:n)) + eps*eye(n-m+1) + randSurf*ones(n-m+1,n-m+1)/(n-m+1);
		k = min(MAXNUMBEREIG,n-m+1)-1;
		[revecs{i},reval]=eigs(Y,sparse(1:n-m+1,1:n-m+1,sum(Y,2)),k,1+eps+randSurf,opt);
		revals{i} = diag(reval);
		Components{i} = sum(revals{i}>1-1e-5);
		[specGaps{i}, aBettis{i}] = sort(abs(diff(revals{i})),'descend');
	end 

	dispInfo=sprintf('%s\t%d\t%d\t%f\t%d\t%f\t%d\t%f', dispInfo,Components{i},aBettis{i}(1), specGaps{i}(1),aBettis{i}(2),specGaps{i}(2),aBettis{i}(3),specGaps{i}(3));
	disp(dispInfo);
	%toc
end
toc 
abetti=aBettis;
%gap = specGaps(1:2,:);

%filename = sprintf('gaps%d.mat',size(X,1));
%save(filename, 'gaps', 'abetti', 'revecs', 'revals');

toc
disp('Done.');
