function mcdiarmid3
% writeup sanity check
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;
erg = [0 1];

%list_totords
test_alfbet
return


NN=7;
maxR = 0;
maxA = [];
maxC = [];

for T=1:1000
  %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);
      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
    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])));
    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 = 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_gam
global lenseq
n = lenseq;
gam = gamma(n,[])
[alf_n, bet_n] = alfbet(n,[]);
Y = sum(gam);
return


function gam = gamma(n,x)
global nstate
gam = zeros(nstate,1);
if (n==1)
  for k=1:nstate
    gam(k) = p_0(k)*(F([1 x]) - F([k x]));
  end
else
  for k=1:nstate
    gam0 = gamma(n-1,[k x]);
    for k1=1:nstate
      gam(k) = gam(k) + p_n(k,k1) * gam0(k1);
    end
  end
end  
return



function [alf, bet] = alfbet(n,x,a,b)
%if (n==2)
%  alf =  a   *F([1 1 x]) - b   *F([2 1 x]);
%  bet = (a-1)*F([1 2 x]) -(b-1)*F([2 2 x]);
if (n==1)
  alf = F([1 x]);
  bet = F([2 x]);
else
  [alf0x, bet0x] = alfbet(n-1,[1 x],a,b);
  [alf1x, bet1x] = alfbet(n-1,[2 x],a,b);
  alf =    a  * alf0x -    b  * bet0x;
  bet = (a-1) * alf1x - (b-1) * bet1x;
end
return




function Hn = H(n,a,b)
nn = 0:n-1;
Hn = sum((a-b).^nn);
return



function DH = dHda(n,a,b)
nn = 1:n-1;
DH = sum(nn.*(a-b).^(nn-1));
return

function DH = dHdb(n,a,b)
nn = 1:n-1;
DH = -sum(nn.*(a-b).^(nn-1));
return

function DH = d2Hda2(n,a,b)
nn = 2:n-1;
DH = sum(nn.*(nn-1).*(a-b).^(nn-2));
return

function DH = d2Hdb2(n,a,b)
nn = 2:n-1;
DH = sum(nn.*(nn-1).*(a-b).^(nn-2));
return

function DH = d2Hdab(n,a,b)
nn = 2:n-1;
DH = -sum(nn.*(nn-1).*(a-b).^(nn-2));
return



%function mathjunk
%F[a_,b_] =  1+(a-b)+(a-b)^2 - (a*(a*f000 - b*f100) + (1-a)*(a*f001 - b*f101) + b*((1-a)*f010 - (1-b)*f110) + (1-b)*((1-a)*f011 - (1-b)*f111))
%DFDA[a_,b_] = D[F[a,b],a]
%DFDB[a_,b_] = D[F[a,b],b]
%D2FDA2[a_,b_] = D[DFDA[a,b],a]
%D2FDB2[a_,b_] = D[DFDB[a,b],b]
%D2FDAB[a_,b_] = D[DFDA[a,b],b]
%return


function H3 = HESS3
f000 = F([1 1 1]);
f001 = F([1 1 2]);
f010 = F([1 2 1]);
f011 = F([1 2 2]);
f100 = F([2 1 1]);
f101 = F([2 1 2]);
f110 = F([2 2 1]);
f111 = F([2 2 2]);

d11 = 2 - 2*f000 + 2*f001;
d22 = 2 - 2*f110 + 2*f111;
d12 = -2 + f010 - f011 + f100 - f101;

H3 = [[d11 d12]; [d12 d22]];

return


function [A,B,C] = ABC(n,a,b)
if (n==1)
  A = 1;
  B = 1;
  C = 1;
else
  [A0, B0, C0] = ABC(n-1,a,b);
  A =    a *A0 +    b *B0;
  B = (1-a)*A0 + (1-b)*B0;
  C = C0 + (a-b)^(n-1);
return

  
end

function test_alfbet
NN = 8;
global nstate
nstate = 2;
global FN
FN = randFN(nstate,NN);
%load BAD
%save BAD FN n x a b A B C C1 d1 d2 d3 a1 b1
%x0=x; a0=a;b0=b; d10=d1;d20=d2;d30=d3;
maxE = 0;
minG = 10000;
for t=1:1000000
  %n = ceil(rand*(NN-2))+1;
  n = round(randinseg(2,NN-1));
  %n = 3;
  %n = 1;
  a = rand; 
  %b = rand;
  b = randinseg(0,a);
  %x = round(rand)+1;
  x = [];
  
  if 0   % test derivatives
    FN = randFN(nstate,NN);
    minG = 10000;
    maxG = -10000;
    % check that the gradient does not vanish
    for a=.001:.01:1
      for b=0.001:.01:a	
	DHa = dHda(n,a,b);
	DHb = dHdb(n,a,b);
	[Dalfa,Dalfb] = DalfDab(n,x,a,b);
	[Dbeta,Dbetb] = DbetDab(n,x,a,b);
	DFpDa = DHa - Dalfa + Dbeta;
	DFpDb = DHb - Dalfb + Dbetb;
	normG = abs(DFpDa) + abs(DFpDb);
	%fprintf('normG = %1.5f \n',normG);
	if (normG<minG)
	  amin = a;
	  bmin = b;
	  minG = normG;
	end
	if (normG>maxG)
	  amax = a;
	  bmax = b;
	  maxG = normG;
	end

	%minG = min(minG,normG);
	%maxG = max(maxG,normG);
	%err1 = totdif(DFpDa,Ndiff('Fpl',3,{n,x,a,b}));
	%err2 = totdif(DFpDb,Ndiff('Fpl',4,{n,x,a,b}));
	%maxE = max(maxE,err1+err2);
      end
    end
    if ( (minG<1e-7) & (maxG>1) )
      'oy vey'
      keyboard
    end
    % Y3[a_,b_] = a*(a*f000 - b*f100) + (1-a)*(a*f001 - b*f101) + b*((1-a)*f010 - (1-b)*f110) + (1-b)*((1-a)*f011 - (1-b)*f111)
    % H3[a_,b_] = 1+(a-b)+(a-b)^2
    % F[a_,b_] = H3[a,b] - Y3[a,b]
    % D[F[a,b],a]
    % D[F[a,b],b]
    
  end
  
  if 1
    for a=.001:.01:1
      for b=0.001:.01:a	
	[alf0, bet0] = alfbet(n,[1 x],a,b);
	[alf1, bet1] = alfbet(n,[2 x],a,b);
    
	nn=0:n-1;
	H0 = sum((a-b).^nn);
	g00 = ( abs(alf0-bet0) <= H0+1e-10);
	g11 = ( abs(alf1-bet1) <= H0+1e-10);
	if ~(g00 & g11)
	  'truly terrible'
	  save BAD FN n x a b
	  keyboard
	end
      end
    end
    
    %INDEED, abandon this approach
    %nn=0:n+1;
    %H1 = sum((a-b).^nn);
    %g01 = ( abs(alf0-bet1) <= H1+1e-10);
    %g10 = ( abs(alf1-bet0) <= H1+1e-10);
    %if ~(g01 & g10)
    %  'abandon approach'
    %  keyboard
    %end
    
    %MY AFFINE CRAP don't work indeed
    %nn=0:n;
    %H1 = sum((a-b).^nn);
    %d1 = alf0 - alf1;
    %d2 = bet1 - bet0;
    %d3 = alf1 - bet1;
    %a1 = 0; b1 = 0;
    %g00 = ( abs(a1*d1 + b1*d2 + d3) <= H1 );
    %a1 = 1; b1 = 0;
    %g10 = ( abs(a1*d1 + b1*d2 + d3) <= H1 );
    %a1 = 1; b1 = 1;
    %g11 = ( abs(a1*d1 + b1*d2 + d3) <= H1 );
    %if ~(g00 & g10 & g11)
    %  'your affine crap dont work'
    %  save BAD FN n x a b
    %  keyboard
    %end
    
  end
  
  
  %fprintf('good t = %d  a=%1.3f b=%1.3f  FF=[%s] \n',t,a,b,num2str(FN{3}));
  fprintf('good t = %d n = %d \n',t,n);
  %fprintf('good t = %d n = %d maxE = %1.9f \n',t,n,maxE);
  %fprintf('good t = %d n = %d minG = %1.7f maxG = %1.7f \n',t,n,minG,maxG);
  %if (rand<0.005)
    FN = randFN(nstate,NN);
  %end
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 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} = '\phi[0\,1\,x#]';
s{2} = '\phi[0\,0\,x^]';
s{3} = '\phi[1\,0\,x^]';
s{4} = '\phi[1\,1\,x#]';

psiq = char(39);
psiqq = [psiq psiq];

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 ');

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

str = strrep(str,'^',psiq);
str = strrep(str,'#',psiqq);
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



function w = Wp(n,x,a,b)
[alf, bet] = alfbet(n,x,a,b);
w = H(n,a,b) - alf + bet;
return

function w = Wn(n,x,a,b)
[alf, bet] = alfbet(n,x,a,b);
w = H(n,a,b) + alf - bet;
return

function dWp = d2Wpda2(n,x,a,b)
Dalfa2 = D2alfDa2(n,x,a,b);
Dbeta2 = D2betDa2(n,x,a,b);
dWp = d2Hda2 - Dalfa2 + Dbeta2;
return

function dWp = d2Wpdb2(n,x,a,b)
Dalfb2 = D2alfDb2(n,x,a,b);
Dbetb2 = D2betDb2(n,x,a,b);
dWp = d2Hdb2 - Dalfb2 + Dbetb2;
return

function dWp = d2Wpdab(n,x,a,b)
Dalfab = D2alfDab(n,x,a,b);
Dbetab = D2betDab(n,x,a,b);
dWp = d2Hdb2 - Dalfb2 + Dbetb2;
return