function mcdiarmid(erg)
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 2;
lenseq = 2;
global A
global C
%global ec
global maxAas
global maxC
%global PAIRS
%PAIRS = [];
%[A,C] = initAC(erg);
%global MASK
%MASK = rand(1,100);
%MASK(1) = 0;
%MASK = ones(1,100);
%MASK = MASK*2;
global FF;
global FN;

%basecase
%list_totords
%test_alfbet
%test12(erg)
test13(erg)
%test14(erg)
%test12_pr(erg)
%test_ibd2(erg);
%test_rab;
return

NN=8;

maxR = 0;
maxA = [];
maxC = [];
global YEC
YEC=[];
for T=1:100000
  for lenseq=3:NN
  %for lenseq=3
    unix('k5init -R');
    %FF = randF(nstate,lenseq);
    FN = randFN(nstate,NN);
    [A,C] = initAC(erg);
    ec = EC(A);
    h = 1-ec;
    for i0=1:lenseq
    %for i0=1
      %for i0=1:n-1
      %lenseq = i0+1;
      ALLSEQ =  list_all_ind(repmat(nstate,1,lenseq));
      ALLSEQ1 = list_all_ind(repmat(nstate,1,lenseq-i0+1));
      XXi = list_all_ind(repmat(nstate,1,i0));
      YY = [];
      %for x1 = 1:nstate
      %for x1 = 1
      for j=1
      %for j=1:size(XXi,1)
        Xi = XXi(j,:);
        %YY(j) = Y(Xi);
        YY(j) = Y_pr(Xi);
	%AB(j) = Y_ab;
      end
      %if (totdif(YY,AB)>1e-10)
      %if (abs(YY)>abs(AB)+1e-10)
	%'dont match up!!!!'
	%'bad bound!!!'
	%keyboard
      %end
      nn = 0:lenseq-i0;
      Bn = sum(h.^nn);
      yy = max(abs(YY));
      g = (yy<=Bn+1e-10);
      rat = yy/Bn;
      maxR = max(maxR,rat);
      %if (rat>maxR)
      %  maxY = yy;
      %  maxA = A;
      %  maxC = C;
      %end
      %%fprintf('T=%d, maxY = %1.4f \n',T,maxY);
      %g = (yy<=1/ec);
      %nn = 0:lenseq-1;
      %g = (yy<=1+h);
      %g = (totdif(YY,YY_pr)<1e-9);
      fprintf('T=%d, n=%d maxR = %1.4f; Y = %1.4f ec = %1.4f g = %d\n',T,lenseq,maxR,yy,ec,g);
      if ~g
        'problem!!!!'
        keyboard
      end
      %if ~mod(T,10)
      %unix('k5init -R');
      %save tmp YEC
      %end
      %YM(i0) = yy;
    end
    %i0 = 2:2:length(PAIRS);
    %i1 = 1:2:length(PAIRS);
    %PP = (PAIRS(i0)-PAIRS(i1))
    %g = all(PP<=(2+ec));
    %if ~g
    %  'problem222222!!!!'
    %  keyboard
    %end  
  end
end
return


function y = evalF(FN,x)
global nstate
lenseq = length(x);
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
y = FF(ind);
return


function r = F(x)
global nstate
%global lenseq
if ~x(end), x=x(1:end-1); end;
lenseq = length(x);
global FN
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
r = FF(ind);
return
%function r = Fn(x,FF,ns,len)
%if ~x(end), x=x(1:end-1); end;
%ind = coord2ind(x,(repmat(ns,1,len)));
%r = FF(ind);
%return
global MASK
K = 100;
%MASK = rand(size(x));
global nstate
% f:\X^n \to \R
% is a function satisfying the idep.bdd. diff. condition
x = x-1;
%r = length(find(x));
%r = length(find(x))+10;
I = find(x);
r = sum(MASK(I));
%x = x/(nstate-1);
%r = sum(x.*MASK(1:length(x)));
%r = mod(prod(x),length(x));
%r = mod(prod(x),nstate);
%r = mod(sum(x),nstate);
%r = mod(sum(x),length(x));
%r = sum(x+1)/length(x);
%r = sum(x+1);
r = r + K;
return

function Ef = E(Xi)
% E[f(X)|Xi]
global ALLSEQ
global A
Ef = 0;
i = length(Xi);
n = size(ALLSEQ,2);
for k=1:size(ALLSEQ,1)
  x = ALLSEQ(k,:);
  p = P(x);
  if (i & (i<n))
    pxi = A(x(i+1),x(i));
    pXi = A(x(i+1),Xi(i));
    p = p / pxi * pXi;
  end
  %pjoin = P(X);
  %pmarg = P(X(1:i));
  %pcond = pjoin/pmarg;
  %p = P(X);
  f = F([Xi x(i+1:end)]);
  Ef = Ef + f*p;
  %Ef = Ef + f*pcond;
end
return


function y = Y(Xi)
% Y1 = E[f(X)|X_i] - E[f(X)|X_{i-1}]
Ef1 = E(Xi);
Ef0 = E(Xi(1:end-1));
y = Ef1 - Ef0;
return


function y = Y_pr(Xi)
% Y1 = E[f(X)|X_i] - E[f(X)|X_{i-1}]
global ALLSEQ1
global A
i = length(Xi);
%n = size(ALLSEQ1,2);
y = 0;

for j=1:size(ALLSEQ1,1)
  x = ALLSEQ1(j,:);
  x = [x 0];
  Xix = [Xi(1:i) x(2:end)];
  Xi_1x = [Xi(1:i-1) x];
  a = pcp(Xix(i+1),Xix(i));
  b = pcp(Xi_1x(i+1),Xi_1x(i));
  p = PCOND(x,Xi(1:i-1));
  p = p / b;

  f0 = F(Xi_1x);
  f1 = F(Xix);
  y = y + p*(a*f1-b*f0);
  %show(X,p,alph,f1,f0);
  %show(X,p,A(X(2),X1),A(X(2),X(1)),f1,f0,X1);
  %show(x(1:end-1),(a*f1-b*f0));
end
return

%function show(X,p,alph,f1,f0)
%function show(X,p,a,b,f1,f0,X1)
function show(X,d)
global PAIRS
global A
global C
%d = (alph*f1-f0);
%d = a*f1-b*f0;
%da = abs(d);
%e = abs(a-b);
%g = (da<=1/(e+eps));
%%fprintf('P[%s] = %1.3f; a =%1.3f f1=%1.3f f0=%1.3f d = %1.3f\n',num2str(X),p,alph,f1,f0,d);

%%fprintf('P[%s]=%1.3e; a=%1.3e f1=%1.3e f0=%1.3e d=%1.3e pd=%1.3e \n',num2str(X),p,alph,f1,f0,d,p*d);
%fprintf('P[%s]=%1.3e; a=%1.3f b = %1.3f e = %1.3f f1=%1.3f f0=%1.3f d=%1.3f pd=%1.3e g=%d\n',num2str(X),p,a,b,e,f1,f0,d,p*d,g);
%%fprintf('p=%1.9f d=%1.4f  pd=%1.9f \n',p,d,p*d);
%fprintf('X=[%s]  (c(x1)=%1.2f)[(a(x2|X1)=%1.2f)(f1=%1.2f)-(a(x2|x1)=%1.2f)(f0=%1.2f)=%1.3f]\n',num2str(X),C(X(1)),A(X(2),X1),f1,A(X(2),X(1)),f0,C(X(1))*(A(X(2),X1)*f1-A(X(2),X(1))*f0));
%fprintf('X=[%s]  (c(%d)=%1.2f)[(a(%d|%d)=%1.2f)(f(%d,%d)=%1.2f)-(a(%d|%d)=%1.2f)(f(%d,%d)=%1.2f)=%1.3f]\n',num2str(X),X(1),C(X(1)),X(2),X1,A(X(2),X1),X1,X(2),f1,X(2),X(1),A(X(2),X(1)),X(1),X(2),f0,C(X(1))*(A(X(2),X1)*f1-A(X(2),X(1))*f0));

fprintf('X=[%s] d=%1.3f \n',num2str(X),d);
PAIRS(end+1)=d;
return

function p = P(x)
global A
global C
if isempty(x)
  p = 1;
else
  p = C(x(1));
end
for j=2:length(x)
  p = p * A(x(j),x(j-1)); % markov case
  %p = p * C(x(j)); % independent case
end
return

function p = PCOND(x,Xi)
global C
if isempty(Xi)
  p = C(x(1));
else
  p = pcp(x(1),Xi(end));
end
for j=2:length(x)
  p = p * pcp(x(j),x(j-1)); % markov case
end
return

function p = pcp(i,j)
% primitive cond prob
global A
if i
  p = A(i,j);
else
  p = 1;
end
return

function [A,C] = initAC(erg)
global nstate

%n1 = floor(nstate/2);
%n2 = nstate - n1;
%A1 = rand(n1);
%A2 = rand(n2);
%A = zeros(nstate);
%A(1:n1,1:n1) = A1;
%A(1+n1:n1+n2,1+n1:n1+n2) = A2;

%A=eye(nstate);
while 1
  A = rand(nstate);
  A = pnormdim(A,1);
  %if (EC(A)>=erg)
  %if (abs(EC(A)-erg)<.05)
  if ( (erg(1)<=EC(A)) & (EC(A)<=erg(2)) )
    break
  end
end
%C = pnormdim(rand(nstate,1),1);
%A(1,1)=10000;
%A(2,2)=10000;
C(1) = 0;
C(2) = 1;

return
%A = pnormdim(A,1);
a = rand;
b = randinseg(0,a);
A(1,1) = a; A(1,2) = b;
A(2,1) = 1-a; A(2,2) = 1-b;

return
%C = pnormdim(C,1);
%A = repmat(C,[1,nstate]); % independent case
e = 1e-8;
A(1,1) = 1-e; A(2,1) = e;
C(2) = 1-e; C(1) = e;

% adjust the diagonal of A
%I = randperm(size(A,1));
%i = I(1);
%A(i,i) = 

return



function FN = randFN(nstate,NN,d)
if ~exist('d','var')
  d=1;
end
FN = {};
for j=1:NN
  FN{j} = randF(nstate,j,d);
end
return

function F = randF(nstate,lenseq,d)
% F is a random function $\X^n \to [0,1]$ that sat. ibd w/c=1
X = list_all_ind(repmat(nstate,1,lenseq));
K = 10;
%K = 1;
%K = lenseq;
%F(1) = 0;
PROGstart(size(X,1),'generating randF...');
for i=1:size(X,1)
%for i=2:size(X,1)
  good = 0;
  while ~good
    F(i) = K*rand;
    %F(i) = round(K*rand);
    good = 1;
    for j = 1:i-1
      good = goodij(X,F,i,j,d);
      if ~good, break, end;
    end
  end
  %fprintf('randF: done %d of %d \n',i,size(X,1));
  PROGupdate(i);
end
PROGend;
return

function b = goodij(X,F,i,j,d)
x = X(i,:);
y = X(j,:);
if (length(find(x~=y))==1)
  b = (abs(F(i)-F(j))<=d);
  %b = (abs(F(i)-F(j))==1);
else
  b = 1;
end
return


function d = delF(F,X)
d = 0;
for i=1:size(X,1)
  for j=i+1:size(X,1)
    x = X(i,:);
    y = X(j,:);
    if (length(find(x~=y))==1)
      d = max(d,abs(F(i)-F(j)));
    end
  end
end
PROGend;  



function checkF
global ALLSEQ
global nstate
N = size(ALLSEQ,1); N = N*(N-1)/2;
PROGstart(size(ALLSEQ,1),'checking...');
for i=1:size(ALLSEQ,1)
  for j=i+i:size(ALLSEQ,1)
    x1 = ALLSEQ(i,:);
    x2 = ALLSEQ(j,:);
    t = find(x1~=x2);
    if (t==1)
      f1 = F(x1);
      f2 = F(x2);
      if (abs(f1-f2)>1)
	fprintf('BAAAAAAAAAD i=%d j=%d \n',i,j);
      end
    end
  end
  PROGupdate(i);
end
PROGend;  

return


function e = EC(A)
% the "ergodicity coefficient" of A
[n,n] = size(A);
d = 0;
for i=1:n
  for j=i+1:n
    x = A(:,i);
    y = A(:,j);
    %d = max(d,sum(abs(x-y)));
    d = max(d,max(abs(x-y)));
  end
end
%e = 1-d/2;
e = 1-d;
%e = d;
return


function testF
global nstate
global lenseq


ALLSEQ =  list_all_ind(repmat(nstate,1,lenseq));
for t=1:10000
  %t
  X = ALLSEQ'-1;
  FF = randF(nstate,lenseq)
  modfit(X,FF);
end
return



function test13(erg)
global nstate
global lenseq
nstate = 2;
lenseq = 3;
global FN;
FN = randFN(nstate,lenseq);
global A
global C
d = 1; 					% max |f-f|
X1 = 1;
x1 = 2;
for t=1:1
  [A,C] = initAC(erg);
  h = 1-EC(A);  % stricture coeff.
  s = 0;
  MSTR = '';
  for x2=1:2
    for x3=1:2
      s=s+ A(x3,x2)*(A(x2,X1)*F([X1 x2 x3]) - A(x2,x1)*F([x1 x2 x3]));
      str = mstr13([x1,x2,x3],X1);
      MSTR = [MSTR ' + ' str];
      %fprintf('%s AB=%1.4f \n',str,(A(x2,X1)*F([X1 x2 x3]) - A(x2,x1)*F([x1 x2 x3])));
      fprintf('%s \n',str);
    end
  end
  %g = (abs(s)<=1/(1-h));
  g = (abs(s)<=(1+h+h^2+h^3));
  fprintf('t=%d, s=%1.3f h=%1.3f  g=%d \n',t,s,h,g);
  if ~g
    keyboard
  end
  MSTR
end
return


function test14(erg)
global nstate
global lenseq
nstate = 2;
lenseq = 4;
global FF;
FF = randF(nstate,lenseq);
global A
global C
d = 1; 					% max |f-f|
X1 = 1;
x1 = 2;
for t=1:1
  %[A,C] = initAC(erg);
  h = 1-EC(A);  % stricture coeff.
  s = 0;
  MSTR = '';
  for x2=1:2
    for x3=1:2
      for x4=1:2
	s=s+ A(x4,x3)*A(x3,x2)*(A(x2,X1)*F([X1 x2 x3 x4]) - A(x2,x1)*F([x1 x2 x3 x4]));
	str = mstr14([x1,x2,x3,x4],X1);
	MSTR = [MSTR ' + ' str];
        fprintf('%s AB=%1.4f \n',str,(A(x2,X1)*F([X1 x2 x3 x4]) - A(x2,x1)*F([x1 x2 x3 x4])));
      end
    end
  end
  %g = (abs(s)<=1/(1-h));
  g = (abs(s)<=(1+h+h^2+h^3+h^4));
  fprintf('t=%d, s=%1.3f h=%1.3f  g=%d \n',t,s,h,g);
  if ~g
    keyboard
  end
  MSTR
end
return



function test23(erg)
global A
global C
d = 1; 					% max |f-f|
X1 = 1;
X2 = 1
x2 = 2;
for t=1:1000
  [A,C] = initAC(erg);
  h = 1-EC(A);  % stricture coeff.
  s = 0;
  MSTR = '';
  for x3=1:2
    %s=s+ A(x2,X1)*(A(x3,X2)*F([X1 X2 x3]) - A(x3,x2)*F([X1 x2 x3]));
    s=s+ (A(x3,X2)*F([X1 X2 x3]) - A(x3,x2)*F([X1 x2 x3]));
    str = mstr23([X1,x2,x3],[X1 X2]);
    MSTR = [MSTR ' + ' str];
    %fprintf('%s AB=%1.4f \n',str,(A(x3,X2)*F([X1 X2 x3]) - A(x3,x2)*F([X1 x2 x3])));
  end
  %g = (abs(s)<=1/(1-h));
  g = (abs(s)<=(1+h));
  fprintf('t=%d, s=%1.3f h=%1.3f  g=%d \n',t,s,h,g);
  if ~g
    keyboard
  end
  %MSTR
end
return


function test12(erg)
global nstate
global lenseq
nstate = 2;
lenseq = 2;
global FN;
FN = randFN(nstate,lenseq);
global A
global C
X1 = 1;
x1 = 2;
for t=1:1
  [A,C] = initAC(erg);
  h = 1-EC(A);  % stricture coeff.
  s = 0;
  MSTR = '';
  for x2=1:2
    s=s+ (A(x2,X1)*F([X1 x2]) - A(x2,x1)*F([x1 x2]));
    str = mstr12([x1,x2],X1);
    MSTR = [MSTR ' + ' str];
    fprintf('%s AB=%1.4f \n',str,(A(x2,X1)*F([X1 x2]) - A(x2,x1)*F([x1 x2])));
  end
  %g = (abs(s)<=1/(1-h));
  g = (abs(s)<=(1+h));
  fprintf('t=%d, s=%1.3f h=%1.3f  g=%d \n',t,s,h,g);
  if ~g
    keyboard
  end
  MSTR
end
return


function test12_pr(erg)
global nstate
global lenseq
nstate = 2;
lenseq = 1;
global FF;
global A
global C
X1 = 1;
x1 = 2;
global STATS
STATS = [];
for t=1:100000
  %global DEB
  %a = DEB.a;
  %b = DEB.b;
  %FF = DEB.F;
  FF = randF(nstate,lenseq);
  a = rand;
  b = rand;
  h = abs(a-b);
  
  %alf1     =    a *(a*F([1 1 1 1]) - b*F([2 1 1 1])) +    b *((1-a)*F([1 2 1 1]) - (1-b)*F([2 2 1 1]));
  %alfF(2) =    a *(a*F([1 1 1 2]) - b*F([2 1 1 2])) +    b *((1-a)*F([1 2 1 2]) - (1-b)*F([2 2 1 2]));
  %betF(1) = (1-a)*(a*F([1 1 2 1]) - b*F([2 1 2 1])) + (1-b)*((1-a)*F([1 2 2 1]) - (1-b)*F([2 2 2 1]));
  %bet2     = (1-a)*(a*F([1 1 2 2]) - b*F([2 1 2 2])) + (1-b)*((1-a)*F([1 2 2 2]) - (1-b)*F([2 2 2 2]));

  %alf1 =    a *F([1 1 1]) -    b *F([2 1 1]);
  %bet2 = (1-a)*F([1 1 2]) - (1-b)*F([2 1 2]);
  %X = abs(alf1+bet2);
  X = abs(a*F(1)-b*F(2));
  %X = abs(alf1);
  B = 1/(1-h);
  %B = 1 + h + h^2 + h^3;
  g = (X<=B);

  fprintf('t=%d, h=%1.3f  g=%d \n',t,h,g);
  
  STATS = [STATS; [X h a b FF]];
  
  
  if ~g
    keyboard
  end

end
return


function str = mstr12_pr(x,X1)
str = '(a$x2$X1 f$X1$x2$x3 - a$x2$x1 f$x1$x2$x3)*a$x3$x2';
str = mstrrep(str,x,X1);
return


function str = mstr12(x,X1)
str = '(a$x2$X1 f$X1$x2 - a$x2$x1 f$x1$x2)';
str = mstrrep(str,x,X1);
return

function str = mstr23(x,X)
str = '(a$x3$X2  f$X1$X2$x3 - a$x3$x2  f$X1$x2$x3)';
str = mstrrep(str,x,X);
return


function str = mstr13(x,X1)
%str = 'a$x3$x2*(a$x2$X1 f$X1$x2$x3 - a$x2$x1 f$x1$x2$x3)';
%str = 'a$x3$x2*(a$x2$X1*F([$X1 $x2 $x3]) - a$x2$x1*F([$x1 $x2 $x3]))';
str = 'a$x3$x2*(a$x2$X1*f$X1$x2$x3 - a$x2$x1*f$x1$x2$x3)';
str = mstrrep(str,x,X1);
return

function str = mstr14(x,X1)
str = 'a$x4$x3*a$x3$x2*(a$x2$X1 f$X1$x2$x3$x4 - a$x2$x1 f$x1$x2$x3$x4)';
str = mstrrep(str,x,X1);
return

function str = mstrrep(str,x,X)
for j=1:length(x)
  xstr = ['$x' num2str(j)];
  %str = strrep(str,'$x1',num2str(x1));
  str = strrep(str,xstr,num2str(x(j)));
  %str = strrep(str,'$x2',num2str(x2));
  %str = strrep(str,'$x3',num2str(x3));
end
for j=1:length(X)
  xstr = ['$X' num2str(j)];
  %str = strrep(str,'$X1',num2str(X1));
  str = strrep(str,xstr,num2str(X(j)));
end

str = strrep(str,'a12','b');
str = strrep(str,'a22','(1-b)');
%str = strrep(str,'a11','1');
%str = strrep(str,'a21','0');
str = strrep(str,'a11','a');
str = strrep(str,'a21','(1-a)');
%str = strrep(str,'a12','b');
%str = strrep(str,'a22','(1-b)');
%str = strrep(str,'a11','a1');
%str = strrep(str,'a21','a2');
%str = strrep(str,'a12','b1');
%str = strrep(str,'a22','b2');

str = strrep(str,'1','0');
str = strrep(str,'2','1');
str = strrep(str,'0-','1-');


return


function r = randinseg(a,b)
r = rand*(b-a) + a;
return

function Y=Y_ab
global lenseq
% nstate = 2
% this is really $Y_1^{(n)}$
n = lenseq;
global A
a = A(1,1); b = A(1,2);
h = abs(a-b);

%[alf1, bet1] = alfbet(n,1);
%[alf2, bet2] = alfbet(n,2);
%Y = a*alf1 + (1-a)*alf2 + b*bet1 + (1-b)*bet2;
%Y = alf1 + bet1;
[alf_n, bet_n] = alfbet(n,[]);
Y = alf_n + bet_n;

%return
d=1;

%  alf =    a  * alf1x +    b  * bet1x;
%  bet = (1-a) * alf2x + (1-b) * bet2x;

%M 
% A1' =    a  * A1 +    b  * B1
% A2' =    a  * A2 +    b  * B2
% B1' = (1-a) * A1 + (1-b) * B1
% B2' = (1-a) * A2 + (1-b) * B2
%  a*A1' + (1-a)*A2' + b*B1' + (1-b)*B2'  
% (a*A1  + (1-a)*A2  + b*B1  + (1-b)*B2)*(a-b)


[alf01, bet01] = alfbet(n-1,1);
[alf02, bet02] = alfbet(n-1,2);
%Y0 = a*alf01 + (1-a)*alf02 + b*bet01 + (1-b)*bet02;
%Z = abs(Y0)*h;
%if (Y-Z>1+1e-10)
%  keyboard
%end

% nope, don't care about this
%M = max(abs(alf01+bet01),abs(alf02+bet02));
%if ( abs(Y) - M*h > d-1e-10 )
%  'bloody hell'
%  keyboard
%end

return
%ab1 = (alf1+bet1);
%ab2 = (alf2+bet2);
%aa12 = (alf1-alf2);
%bb12 = (bet1-bet2);
%aab1 = abs(alf1+bet1);
%aab2 = abs(alf2+bet2);

%if ( a^2 + 2*b -b^2 <= 1/(1-h) )
%    'good'  
%else
%  'bad'
%  keyboard
%end


if (n==3)
  Y0 = alf1 + bet1;
else
  [alf01, bet01] = alfbet(n-1,1);
  [alf02, bet02] = alfbet(n-1,2);
  Y0 = a*alf01 + (1-a)*alf02 + b*bet01 + (1-b)*bet02;
end

%if ((a>b) & (Y0>0) & (Y>0) & (Y>Y0) )
  %if (max(ab1,ab2)>d/(1-a+b)+1e-10) % FAAAAAAAAAALLLSEEEE
  %if (max(aa12,bb12)>d+1e-10)  % FAAAAAAAAAALLLSEEEE
  %if (abs(a*aa12+b*bb12)>d/(1-a+b)+1e-10)
%  if (abs( a*(alf1-alf2)+b*(bet1-bet2)+(alf2+bet2)    )>d/(1-a+b)+1e-10)
%    keyboard
%  else
%    'good'  
%  end
%end

if 0
  
  if (n>3)
    %if ((a>b) & (Y0>0) & (Y>0) )
    %if (max(ab1,ab2)>1+1e-10)
    if (abs( (1-a)*ab2+b*ab1 ) > 1+1e-10)
      keyboard
    end
    %end
  end
end
return


function [alf, bet] = alfbet(n,x)
global A
a = A(1,1); b = A(1,2);
if (n==2)
  alf =    a *F([1 1 x]) -   b *F([2 1 x]);
  bet = (1-a)*F([1 2 x]) -(1-b)*F([2 2 x]);
else
  [alf1x, bet1x] = alfbet(n-1,[1 x]);
  [alf2x, bet2x] = alfbet(n-1,[2 x]);
  alf =    a  * alf1x +    b  * bet1x;
  bet = (1-a) * alf2x + (1-b) * bet2x;
end
return


function [dA, dB, dAB] = delAB(n,a,b)
if (n==1)
  dA = 1;
  dB = 1;
  dAB = 0;
else
  [dA0, dB0, dAB0] = delAB(n-1,a,b);
  dA =    a *dA0 +    b *dB0;
  dB = (1-a)*dA0 + (1-b)*dB0;
  
  h = abs(a-b);
  %nn = 0:n-2;
  %nn = 0:n-1;
  nn = 0:n+1;
  Hn = sum(h.^nn);
  H = 1/(1-h);

  %dAB = 2*abs(a*dA0-b*dB0) + Hn;  %not tight enough!!
  %dAB = abs( 2*(a*dA0-b*dB0) + Hn ); % plain wrong!
  %dAB = h*(dA0+dB0) + Hn; %not tight enough!!
  dAB = Hn; % too tight
  %dAB = H; 
  
end
return


function test_rab
global nstate
nstate = 2;

d=1;

for t=1:100000
  a = rand;
  b = rand;
  h = abs(a-b);
  M = d*(1+h+h^2);

  %ALF = randFN(nstate,1,d*(a+b));
  %BET = randFN(nstate,1,d*(2-a-b));
  %alf0 = evalF(ALF,1);
  %alf1 = evalF(ALF,2);
  %bet0 = evalF(BET,1);
  %bet1 = evalF(BET,1);
  %g = ( abs(a*alf0 + b*bet0 + (1-a)*alf1 + (1-b)*bet1) <= M + 1e-10 );
  %if ~g
  %  'not quite'
  %  keyboard
  %end
  
  g = ( d*(1+h)*max(a+b,2-a-b) <= M + 1e-10 );
  if ~g
    'nope'
    keyboard
  end
  t
end
return


function test_ibd2(erg)
global nstate
nstate = 2;
NN=3;
K=10;
global FN;
global A;

global STATS
%STATS = [];
d=1;

maxR = 0;

for t=1:1000000
  [A,C] = initAC(erg);
  a = A(1,1); b = A(1,2);
  h = abs(a-b);
  if 1
    n = ceil(rand*7)+1;
    %n = 2;
    NN = 8;
    nn = 0:n-1;
    Hn = sum(h.^nn);
    H1 = Hn + h^n;
    
    while 1
      d = rand;

      %phi00 = randinseg(-K,K);
      %phi01 = randinseg(-K,K);
      %phi10 = randinseg(-K,K);
      %phi11 = randinseg(-K,K);
      alf0 = randinseg(-K,K);
      alf1 = randinseg(-K,K);
      bet0 = randinseg(-K,K);
      bet1 = randinseg(-K,K);
      g = 1; 

      %g = g & (abs(phi00-phi10)<=d*(1+h));
      %g = g & (abs(phi01-phi11)<=d*(1+h));
      %g = g & (abs(phi00-phi01)<=d*(a+b));
      %g = g & (abs(phi10-phi11)<=d*(2-a-b));
      %g = g & (abs(phi00-phi11)<=d*(2+h));
      %g = g & (abs(phi01-phi10)<=d*(2+h));

      %g = g & (abs(phi00-phi10)<=d);
      %g = g & (abs(phi01-phi11)<=d);
      %g = g & (abs(phi00-phi01)<=d);
      %g = g & (abs(phi10-phi11)<=d);

      [dA, dB, dAB] = delAB(n,a,b);
      
      g = g & (abs(alf0+bet0)<=d*Hn);
      g = g & (abs(alf1+bet1)<=d*Hn);
      g = g & (abs(alf0-alf1)<=d*dA);
      g = g & (abs(bet0-bet1)<=d*dB);
      g = g & (abs(alf0+bet1)<=d*dAB);
      g = g & (abs(alf1+bet0)<=d*dAB);
      %g = 1;
      if g, break, end;
    end
    %g = ( abs(a*phi00 - b*phi10 + (1-a)*phi01 - (1-b)*phi11) <= d*M+1e-10 );
    %g = ( abs(a*phi00 - b*phi10 + (1-a)*phi01 - (1-b)*phi11) <= d*(1+h)+1e-10 );
    g = ( abs(a*alf0 + b*bet0 + (1-a)*alf1 + (1-b)*bet1) <= d*H1+1e-10  );
    %STATS = [STATS; [a b phi00 phi01 phi10 phi11 g]];
    if ~g
      'kol ha zayin!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
      %'gotcha mofo'
      save BAD a b d alf0 bet0 alf1 bet1
      keyboard
    end
    %alf0 =    a *F([1 1 1]) -    b *F([2 1 1]);
    %alf1 =    a *F([1 1 2]) -    b *F([2 1 2]);
    %bet0 = (1-a)*F([1 2 1]) - (1-b)*F([2 2 1]);
    %bet1 = (1-a)*F([1 2 2]) - (1-b)*F([2 2 2]);
    
  end
  
  if 1
    NN = 7;
    %NN = 3;
    n = ceil(rand*(NN-2))+1;
    %n = 3;
    nn = 0:n-1;
    Hn = sum(h.^nn);
    Hn1 = Hn + h^n;

    FN = randFN(nstate,NN,d);

    %alf0 =    a *F([1 1 1]) -    b *F([2 1 1]);
    %alf1 =    a *F([1 1 2]) -    b *F([2 1 2]);
    %bet0 = (1-a)*F([1 2 1]) - (1-b)*F([2 2 1]);
    %bet1 = (1-a)*F([1 2 2]) - (1-b)*F([2 2 2]);
    [alf0, bet0] = alfbet(n,1);
    [alf1, bet1] = alfbet(n,2);
    %[Alf0, Bet0] = alfbet(n+1,1);
    %[Alf1, Bet1] = alfbet(n+1,2);
    %[Alf , Bet] = alfbet(n+1,[]);
 
    [dA, dB, dAB] = delAB(n,a,b);
    
    g = 1;
    g = g & (abs(alf0+bet0)<=d*Hn1);
    g = g & (abs(alf1+bet1)<=d*Hn1);
    g = g & (abs(alf0-alf1)<=d*dA);
    g = g & (abs(bet0-bet1)<=d*dB);
    g = g & (abs(alf0+bet1)<=d*dAB);
    g = g & (abs(alf1+bet0)<=d*dAB);
    if ~g
      %'gotcha'
      'vse idet v pizdu'
      save BAD a b d alf0 bet0 alf1 bet1
      keyboard
    end
    Dab = max(abs(alf0+bet1),abs(alf1+bet0));
    maxR = max(maxR,Dab/d/dAB);
    
    %alf =    a *F([1 1]) -    b *F([2 1]);
    %bet = (1-a)*F([1 2]) - (1-b)*F([2 2]); 
    
    %dab = max(abs(alf0-bet1),abs(alf1-bet0));
    %STATS = [STATS; [Dab h a b alf0 alf1 bet0 bet1 d]];
  end
  fprintf('GOOD: t=%d  n=%d  maxR = %1.3f \n',t,n,maxR);
end

return


function test_ibd(erg)
global nstate
nstate = 2;
NN=3;
K=3;
global FN;
global A;

global STATS
STATS = [];
d=1;

II = list_all_ind([2 2 2 2 2 2]) - 1;

for t=1:100000
  FN = randFN(nstate,NN+K);
  [A,C] = initAC(erg);
  a = A(1,1); b = A(1,2);
  h = abs(a-b);
  M = (1+h+h^2);

  for n=2:NN-1
    ell = NN-n+K;
    x0 = round(rand(1,ell));
    i = ceil(rand*ell);
    x1 = x0;
    x1(i) = 1-x0(i);
    x0=x0+1; x1=x1+1;
    
    [alf0, bet0] = alfbet(n,x0);
    [alf1, bet1] = alfbet(n,x1);

    
    TERMS = [a*alf0, -a*alf1, b*bet0, -b*bet1, alf1, bet1];
    g = (abs(sum(TERMS)) <= M+1e-10 );
    if ~g   % sanity check
      'major problem 1'
      keyboard
    end

    ii = 13; % bB
    ptn = II(ii,:);
    i1 = find(ptn);
    i2 = find(~ptn);
    g13 = ( b*d*(2-a-b) + abs(sum(TERMS(i2))) <= M+1e-10 );
    
    
    ii = 49; % aA
    ptn = II(ii,:);
    i1 = find(ptn);
    i2 = find(~ptn);
    g49 = ( a*d*(a+b) + abs(sum(TERMS(i2))) <= M+1e-10 );
    
    if ~(g13 | g49)
      'abandon approach'
      keyboard
    end
    
    
    %stat = [];
    %for i=1:size(II,1)
    %  ptn = II(i,:);
    %  i1 = find(ptn);
    %  i2 = find(~ptn);
    %  g = ( abs(sum(TERMS(i1))) + abs(sum(TERMS(i2))) <= M+1e-10 );
    %  stat(end+1) = g;
    %end
    
      

    %da = abs(alf1-alf2);
    %db = abs(bet1-bet2);
    %M = abs(a*alf1+b*bet1-a*alf2-b*bet2);
    % 
    %if (da>d*(a+b)+1e-10)
    %  'bad da'
    %  keyboard
    %end
    %if (db>d*(2-(a+b))+1e-10)
    %  'bad db'
    %  keyboard
    %end
    %
    %if ( a*(a+b) - b*(2-a-b) + (1+h) > (1+h+h^2)+1e-10 )  % FALSE
    %if (abs( a*(alf1-alf2) + b*(bet1-bet2) + (alf2+bet2)  ) > (1+h+h^2)+1e-10 )  % TRUE
    %if ( a*abs(alf1-alf2) + b*abs(bet1-bet2) + abs(alf2+bet2) > (1+h+h^2)+1e-10 )
    %if ( dab > a^2+2*b-b^2 + 1e-10)
    %if ( abs( a*(alf1-alf2) + b*(bet1-bet2)) + abs(alf2+bet2)   > (1+h+h^2)+1e-10 )  
      %'oh well'
      %'real bad'
      %'kinda bad'
      %'pretty bad'
      %keyboard
    %else
    %  'so far so good'
    %end
    
    %if ( M <= h^2 + 1e-10)
    %  'candy'
    %else
    %  
    %end

    
    %STATS = [STATS; [h dab a b alf1 alf2 bet1 bet2]];
    %STATS = [STATS; stat];
  end
  t
end
  
return


function test_alfbet
global nstate
global lenseq
nstate = 2;
lenseq = 3;
global FF;

global STATS
STATS = [];
global DEB

maxR = 0;
k = 1;
for t=1:100000
  %a = DEB.a;
  %b = DEB.b;
  %FF = DEB.F;
  a = rand;
  %a=1;
  %b = randinseg(0,a);
  b = rand;
  %b=0;
  FF = randF(nstate,lenseq);

  h = abs(a-b);
  
  if 0
    K = k*max(abs(FF));
    %B = K*rand;
    %B = rand; 				% initialize the bound
    %B = 1+abs(a-b);
    while 1
      %K = 1+h;
      B0 = 1/(1-h);
      %B0 = K; %baaaaaaaaad
      %B0 = 1+h+h^2;
      alf(1) = randinseg(-K,K);
      alf(2) = randinseg(-K,K);
      bet(1) = randinseg(-K,K);
      bet(2) = randinseg(-K,K);

      %B1 = max([abs(alf+bet) abs(alf(end:-1:1)+bet)]);
      B1 = max(abs(alf+bet));
      
      %SS = 2*(list_all_ind([2 2 2 2]) - 1) -1;
      %for j=1:size(SS,1)
      %s = SS(:,j);
      %
      %bet(1) = B - alf(1);
      %alf(2) = alf(1) + (a+b);
      
      
      %alf(2) = alf(1) + (a+b);
      %alf(2) = alf(1) + (a+b);
      %bet(1) = B - alf(1);
      %bet(2) = bet(1) + (a+b);
      g = goodalfbet(alf,bet,K,1,B0,B1,a,b);
      %g = 1;
      if g, break, end;
    end
    % pure form
    X = abs(a*alf(1) + (1-a)*alf(2) + b*bet(1) + (1-b)*bet(2));
    	
    Y = B1*h + 1;
    %Y = 1+h+h^2+h^3;
    %X1 = abs( a*alf(1) + (1-a)*(alf(1)+(a+b)) + b*bet(1) + (1-b)*(bet(1)+(a+b)) );
    %M F[a1_,b1_]=(B*(a-b)+1)^2 - (a*a1 + (1-a)*a2 + b*b1-(1-b)*b2 )^2
    %M a*a1 + (1-a)*(a1+(a+b)) + b*b1 + (1-b)*(b1+(a+b)) 
    %Y = B*abs(a-b)+1;
    g = (X<=Y);
    %g = (X<=X1);
    if ~g
      'bad1'
      keyboard
    end
  end
  %% FALSE??


  alf(1) =    a *F([1 1 1]) -   b *F([2 1 1]);
  alf(2) =    a *F([1 1 2]) -   b *F([2 1 2]);
  bet(1) = (1-a)*F([1 2 1]) -(1-b)*F([2 2 1]);
  bet(2) = (1-a)*F([1 2 2]) -(1-b)*F([2 2 2]);
  
  
  %alf(1) =    a *(a*F([1 1 1 1]) - b*F([2 1 1 1])) +    b *((1-a)*F([1 2 1 1]) - (1-b)*F([2 2 1 1]));
  %alf(2) =    a *(a*F([1 1 1 2]) - b*F([2 1 1 2])) +    b *((1-a)*F([1 2 1 2]) - (1-b)*F([2 2 1 2]));
  %bet(1) = (1-a)*(a*F([1 1 2 1]) - b*F([2 1 2 1])) + (1-b)*((1-a)*F([1 2 2 1]) - (1-b)*F([2 2 2 1]));
  %bet(2) = (1-a)*(a*F([1 1 2 2]) - b*F([2 1 2 2])) + (1-b)*((1-a)*F([1 2 2 2]) - (1-b)*F([2 2 2 2]));
  
  
  %da = abs(alf(1)-alf(2));
  %db = abs(bet(1)-bet(2));
  %if (max(da,db)>(a+b))
  %  % FAAAAAAAAAALSEEEEEEEEEEEEEe
  %  keyboard
  %end
  
  Mc = max(abs(alf+bet));                           % coupled max
  m1 = abs(alf(1)+bet(1));
  m12 = abs(alf(1)+bet(2));
  mc = min(abs(alf+bet));                           % coupled min
  Mi = max([abs(alf+bet) abs(alf([2 1])+bet)]);   % indep max
  mi = max([abs(alf+bet) abs(alf([2 1])+bet)]);   % indep min
  DEL = abs(a*alf(1) + (1-a)*alf(2) + b*bet(1) + (1-b)*bet(2));
  BND = m12*h + 1;  %*****
  g = (DEL<=BND+1e-10);  
  %g = (Mi <= Mc*h + 1);    % FAAAAAAAAAAALLLSSSSSSEEEE
  if ~g
    keyboard
  end
  %m = max([alf bet]);
  %MAH = [MAH;[m,h]];  
  
  %g = ( a*alf(1) + (1-a)*alf(2) + b*bet(1) + (1-b)*bet(2) <= B*(a-b) + 1 );
  %X = abs(a*alf(1) + (1-a)*alf(2) + b*bet(1) + (1-b)*bet(2));
  %Y = B*(a-b) + 1;
  %X = a*alf(1) - a*bet(2) - b*alf(2) + b*bet(1);
  %B1 = B*(a-b);
  %g = (X<=B1);
  %X = a*alf(1) + (1-a)*alf(2) + b*bet(1) + (1-b)*bet(2);
  %Y = a*alf(1) + alf(2) - a*(bet(2)-B) + b*bet(1) + bet(2) - b*(alf(2)-B);
  %Y = a*(alf(1)-alf(2)) + b*(alf(2)-alf(1)) + B;
  %g = (X<=Y);
  %maxR = max(maxR,X/Y);
  %if (maxR>1)
  %  keyboard
  %end
  %F = 1 + B*(a-b) - (a*alf(1) + (1-a)*alf(2) + b*(B-alf(1)) + (1-b)*(B-alf(2)));
  %g = (F>=0);
  %fprintf('t=%d  alf=[%s]  bet=[%s]  a=%1.3f  b=%1.3f  g=%d maxR = %1.3f\n',t,num2str(alf),num2str(bet),a,b,g,maxR);
fprintf('t=%d  alf=[%s]  bet=[%s]  a=%1.3f  b=%1.3f  \n',t,num2str(alf),num2str(bet),a,b);
end
return


function g = goodalfbet(alf,bet,K,d,B0,B1,a,b)
g = 1;
%g = g & (max(abs(alf+bet))<=B);
g = g & (max(abs(alf+bet))<=B0);
%g = g & (B1<=B0);  %*********
%g = g & (B1<=K);
%ab = max(a+b,2-(a+b));
%g = g & (max(alf) <= K);
%g = g & (max(bet) <= K);
g = g & (abs(alf(1)-alf(2)) <= d*( a*(a+b) + b*(2-a-b)  ));
g = g & (abs(bet(1)-bet(2)) <= d*( (1-a)*(a+b) + (1-b)*(2-a-b) ));
return


function basecase
nstate = 2;
NN = 2;
PERMS = {};
global FN;
for t= 1:10000
  a = rand;
  b = rand;
  d = 10*rand;
  FN = randFN(nstate,NN,d);
  %[tmp,inds] = sort(FN{2});
  %PERMS{end+1} = num2str(inds);
  %PERMS = unique(PERMS);
  %nperms = length(PERMS)
  %if (nperms==24)
  %  break
  %end
  d1 = F([1 1]) - F([1 2]);
  d2 = F([2 2]) - F([2 1]);
  d3 = F([1 2]) - F([2 2]);
  g = ( abs(a*d1+b*d2+d3) <= d*(1+abs(a-b)) );
  if ~g
    'major fuck'
  end
  fprintf('good %d \n',t);
end
return


function list_totords
perms = generate_perms(4);
% want: 1<2, 3<4
for j=1:size(perms,1)
    prm = perms(j,:);
    i1 = find(prm==1);
    i2 = find(prm==2);
    i3 = find(prm==3);
    i4 = find(prm==4);
    %if ((i1<i2)&(i3<i4)) % good perm
        fprintf('%s \n',prm2str(prm));
        %end
end
return

function str = prm2str(prm)
for j=1:4
    i(j) = find(prm==j);
end

s={};
%s{1} = '\phi[1\,0]';
%s{2} = '\phi[1\,1]';
%s{3} = '\phi[0\,1]';
%s{4} = '\phi[0\,0]';
s{1} = 'phi00';
s{2} = 'phi01';
s{3} = 'phi10';
s{4} = 'phi11';

str = num2str(prm);
str = strrep(str,' ','');
k=length(str);
tmp=repmat('<',1,k);
tmp = [str;tmp];
str = tmp(:)';
str = str(1:end-1);
%str = strrep(str,'<',' \leq ');
str = strrep(str,'<',' < ');

for j=1:4
    str=strrep(str,num2str(j),s{j});
end

return

%-------obsolete----------

function y = Y_old(X1)
% Y1 = E[f(X)|X1] - E[f(X)]
global ALLSEQ_old
global ALLSEQ1_old

EfX = 0;
for j=1:size(ALLSEQ_old,1)
  X = ALLSEQ_old(j,:);
  p = P(X);
  f = F(X);
  EfX = EfX + f*p;
end
EfX

EfX1 = 0;
for j=1:size(ALLSEQ1_old,1)
  X2 = ALLSEQ1_old(j,:);
  X = [X1 X2];
  pjoin = P(X);
  pmarg = P(X1);
  pcond = pjoin/pmarg;
  f = F(X);
  EfX1 = EfX1 + f*pcond;
end
EfX1

y = EfX1 - EfX;
return

function y = Y_pr_old(X1)
% Y1 = E[f(X)|X1] - E[f(X)]
global ALLSEQ
global A

y = 0;

for j=1:size(ALLSEQ,1)
  X = ALLSEQ(j,:);
  %alph = A(X(2),X1)/A(X(2),X(1));
  a = A(X(2),X1);
  b = A(X(2),X(1));
  p = P(X)/b;
  f0 = F(X);
  f1 = F([X1 X(2:end)]);
  %y = y + p*(alph*f1 - f0);
  y = y + p*(a*f1-b*f0);
  %y = y + p*abs(a*f1-b*f0);
  %show(X,p,alph,f1,f0);
  show(X,p,A(X(2),X1),A(X(2),X(1)),f1,f0,X1);
end
return


%function F = populateF
% f:\X^n \to \R
% is a function satisfying the idep.bdd. diff. condition
%global nstate
%global lenseq
%global ALLSEQ
%global F
%seq = ALLSEQ(1,:)
%F = zeros(size(ALLSEQ,1),1);
%function populate_sub(seq,val)
%return


%function [N,K1,K2,A] = modfit(X,F)
function modfit(X,F)
%function [N,K,A] = modfit(X,F)
global nstate
global lenseq
good = 0;
%AA = list_all_ind(repmat(2*(nstate-1)+1,1,lenseq))-nstate;
AA = list_all_ind(repmat(7,1,lenseq))-4;
%A = repmat(nstate+1,1,lenseq);
mf = max(abs(F))+1;
MF = 2*mf;
PP = primes(MF);
ii = find(PP>=mf);
p0 = PP(ii(1));

for N=1:p0
%for N=1:2
  %for N=nstate
  for K1=-max(F):max(F)
    %for K2=-max(F):max(F)
      for a=1:size(AA,1)
	A = AA(a,:);
	%good = all(F==s*mod(A*X+K,N));
	%good = all(F==max(mod(A*X+K1,N),0)+K2);
	good = all(F==max(mod(A*X+K1,N),0));
	%good = all(F==max(mod(A*X,N),0));
	if good,break,end;
      end
      %if good,break,end;
    %end
    if good,break,end;
  end
  if good,break,end;
end
if ~good
  'cant do it!!!'
  keyboard
end
fprintf('A=[%s]  N = %d \n',num2str(A),N);
return


function test12_old(erg)
f11=10;f21=11;f12=11;f22=12;
%In[1]:= a*f11-b*(f11+d)+(1-a)*f21-(1-b)*(f21+d)
%Out[2]= -d + (a - b) (f11 - f21)
%
%In[5]:= (a*f11-b*f21)+((1-a)*f12-(1-b)*f22)
%Out[6]= a (f11 - f12) + f12 - b f21 - f22 + b f22

d = 1; 					% max |f-f|
for t=1:10000
  [A,C] = initAC(erg);
  a = A(1,1); b = A(1,2);
  h = EC(A);
  %y = (A(1,1)*f11-A(1,2)*f21)+(A(2,1)*f12-A(2,2)*f22);
  y = (a*f11-b*f21)+((1-a)*f12-(1-b)*f22);
  y = abs(y);
  %g = (abs(y)<=1/h);
  y1 = a*(f11 - f12) + f12 - b*f21 - f22 + b*f22;
  y1 = abs(y1);
  %|y| = |y1| check!
  %g0 = (abs(y-y1)<1e-9);
  %if ~g0
  % g0
  %  keyboard
  %end
  %g = (y1<=1/h);
  d1 = rand;
  %g = (y<=1/(1-abs(a-b)));
  y2 = abs(a*d1-b*d1+d1);
  y2 = abs(y2);
  %g = (y2<=abs(d1)/(1-abs(a-b)));
  %g = (abs(a-b+1)<=1/(1-abs(a-b)));
  g = (abs(a-b+1)<=1/(1-abs(a-b)));
  fprintf('t=%d, y=%1.3f h=%1.3f  g=%d \n',t,y,h,g);
  if ~g
    keyboard
  end
end
return