function erasure(M,N)
curd = pwd;
cd ../string_proc/combinatorics/
addpcurdir
cd(curd)

T = 1e4;
% indep
%M = 10;
%N = 6;
alf = 'ab';
p = (M-N)/M;

nseqi = [];
nseqd = [];
NI = zeros(M+1,1);
ND = zeros(M+1,1);
for t=1:T
  x = get_datum(alf,M);
  
  y = get_datum(alf,N);
  nstr = count_subseq(x,y);
  if nstr
    if nstr > length(NI);
      NI(nstr)=0;
    end
    NI(nstr) = NI(nstr) + 1;
  end
  nseqi(end+1) = nstr;
  
  y = rand_erase(x,N);
  nstr = count_subseq(x,y);
  if nstr
    if nstr > length(ND)
      ND(nstr)=0;
    end
    ND(nstr) = ND(nstr) + 1;
  end
  nseqd(end+1) = nstr;

  if 0
  figure(1)
  hist(nseqi,length(NI)/10);
  figure(2);
  hist(nseqd,length(ND)/10);
end
  
  if 1
  p1 = NI/sum(NI);
  p2 = ND/sum(ND);
  pa1 = cumsum(p1);
  pa2 = cumsum(p2);
  if 0
    pa = zeros(M,1);
    pa1 = zeros(M,1);
    pa2 = zeros(M,1);
    for a=1:M
      %pa(a) = sum(p1(a+1:end)) + sum(p2(1:a));
      pa1(a) = sum(p1(1:a));
      pa2(a) = sum(p2(1:a));
    end
    %plot(pa,'.');
  end

  figure(1)
  plot(pa1,'.');
  axis([0 length(NI) 0 1]);
  figure(2);
  plot(pa2,'.');
  axis([0 length(ND) 0 1]);
end
  
  drawnow
end

return

function y = rand_erase(x,N)
ii = randperm(length(x));
ii = ii(1:N);
ii = sort(ii);
y = x(ii);
return
y = [];
for i=1:length(x)
  if rand > p
    y(end+1) = x(i);
  end
end
return

function x = get_datum(alf,n)
na = length(alf);
x = alf(ceil(rand(1,n)*na));
return