function erasediscr(M,N)
curd = pwd;
cd ../string_proc/combinatorics/
addpcurdir
cd(curd)
alf = 'ab';

E0 = indepExp(M,N);

T=100;
acc = 0;
ntr = M;

while 1
  r = rand;
  k = r*M;
  for t=1:T
    b = round(rand); % indep?
  
    X = 0;
    for i=1:ntr
      x = get_datum(alf,M);
      if b
        y = get_datum(alf,N); 
      else
        y = rand_erase(x,N);
      end
      X = X + count_subseq(x,y);
    end
    p = (X/ntr <  k*E0);
  
    acc = acc + (p==b);
  end
  acc = acc/T;
  hold on
  plot(k,acc,'.')
  axis([0 M 0 1])
  drawnow
  
end

return

maxr = 0;
for m=1:M
  for n=1:m
    K = total_subsec_cnt(m,n);
    r = choose(m,n) * 2^(m-n) / K * (m+1);
    maxr = max(r,maxr);
    fprintf('m = %2d   n = %2d    rat = %1.12f  mrat = %1.12f \n',m,n,r,maxr);
  end
end

return

T = 1e10;

%E = indepExp(M,N);
E = eraseEmp(M,N);
alf = 'ab';
p = (M-N)/M;
nn = 0;

global m n P1
P0 = eraseproc(M,N);
P1 = P0*0;
Cemp = P1;
m = 2;
n = M + N;


for t=1:T
  x = get_datum(alf,M);
  %y = get_datum(alf,N); 
  %y = rand_erase(x,p);
  y = rand_erase(x,N);
  
  %if veq(x,'abb') & veq(y,'ba')
  %  keyboard
  %end

  xi = c2i([x y]-'a');
  Cemp(xi) = Cemp(xi)+1;
  P1 = Cemp/sum(Cemp);
  tv = .5*totdif(P0,P1);

  nstr = count_subseq(x,y);
  nn = nn + nstr;
  fprintf('E = %1.8f     A = %1.8f  tv = %1.8f \n',E,nn/t,tv);
end
return



function E = indepExp(M,N)
E = exp(logchoose(M,N) + N*log(1/2));
return

function E = eraseEmp(M,N)
global m n
m = 2;
n = M + N;
E = 0;
Csum = 0;
for xi = 1:m^n
  xy = i2c(xi);
  x = xy(1:M);
  y = xy(M+1:M+N);
  Csum = Csum + count_subseq(x,y)^2;
end
E = (1/2)^M * Csum / choose(M,N);
return

function E = eraseExp(M,N)
p = 0;
for f=0:M-length(y)
  %x = logchoose(M,N-f) - logchoose(M,N) + f*log(1/2);
  x = logchoose(M,N-f) + f*log(1/2);
  p = p + exp(x);
end
E = p;
return




function x = get_datum(alf,n)
na = length(alf);
x = alf(ceil(rand(1,n)*na));
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