function mcd_wrc
% writeup sanity check
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 2;
lenseq = 5;
global AA
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];

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

for T=1:1000
  for lenseq=1:NN
  %for lenseq=3
    FN = randFN(nstate,NN);
    [A,C] = initAC(erg);
    [AA,C] = initAAC(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 = [];
      AB = [];
      %for x1 = 1:nstate
      %for x1 = 1
      %for j=1
      for j=1:size(XXi,1)
        Xi = XXi(j,:);
        if (Xi(end)==1)
            YY(j) = Y(Xi);
            %YY_pr(j) = Y_pr(Xi);
	        AB(j) = Y_ab(Xi(1:end-1));
        end
      end
      if (totdif(YY,AB)>1e-10)
      U_L = max(YY)-min(YY)
      %if (abs(YY)>abs(AB)+1e-10)
	    'dont match up!!!!'
    	%'bad bound!!!'
	    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 A = initA(sc)
global nstate
while 1
  A = rand(nstate);
  A = pnormdim(A,1);
  if ( (sc(1)<=SC(A)) & (SC(A)<=sc(2)) )
    break
  end
end
return

function [AA,C] = initAAC(sc)
global NN
global nstate
AA = zeros(nstate,nstate,NN-1);
A0 = initA(sc);
for i=1:NN-1
  %AA(:,:,i) = initA(sc);
  AA(:,:,i) = A0;
end
%C(1) = 0;
%C(2) = 1;
C = pnormdim(rand(nstate,1),1);
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_ab(X)
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);
i = length(X)+1;

[alf_n, bet_n] = alfbet(n-i+1,[],a,b,X);
Y = alf_n - bet_n;
return


function [alf, bet] = alfbet(n,x,a,b,X)
global A
%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]);
  lx = length(X);
  if ~lx
      p0 = p_0(2);
  else
      p0 = A(2,X(end));
  end
  alf = p0*F([X 1 x]);
  bet = p0*F([X 2 x]);
else
%    a = p_n(1,1,n-1);
%    b = p_n(1,2,n-1);

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


function p=p_0(k)
global C
p=C(k);
return

function p=p_n(k,k1,n)
global AA
p=AA(k,k1,n);
return

function [Dalfa,Dalfb] = DalfDab(n,x,a,b)
if (n==1)
  Dalfa = 0;
  Dalfb = 0;
else
  [alf0x, bet0x] = alfbet(n-1,[1 x],a,b);

  [Dalf0xa,Dalf0xb] = DalfDab(n-1,[1 x],a,b);
  [Dbet0xa,Dbet0xb] = DbetDab(n-1,[1 x],a,b);
  
  Dalfa = a*Dalf0xa + alf0x - b*Dbet0xa;
  Dalfb = a*Dalf0xb - bet0x - b*Dbet0xb;
end
return

function [Dbeta,Dbetb] = DbetDab(n,x,a,b)
if (n==1)
  Dbeta = 0;
  Dbetb = 0;
else
  [alf1x, bet1x] = alfbet(n-1,[2 x],a,b);

  [Dalf1xa,Dalf1xb] = DalfDab(n-1,[2 x],a,b);
  [Dbet1xa,Dbet1xb] = DbetDab(n-1,[2 x],a,b);
  
  Dbeta = (a-1)*Dalf1xa + alf1x - (b-1)*Dbet1xa;
  Dbetb = (a-1)*Dalf1xb - bet1x - (b-1)*Dbet1xb;
end
return


function Dalfa2 = D2alfDa2(n,x,a,b)
if (n==1)
  Dalfa2 = 0;
else
  D2alf0x = D2alfDa2(n-1,[1 x],a,b);
  D2bet0x = D2betDa2(n-1,[1 x],a,b);

  [Dalf0xa,Dalf0xb] = DalfDab(n-1,[1 x],a,b);
  
  Dalfa2 = a*D2alf0x + 2*Dalf0xa - b*D2bet0x;
end
return


function Dbeta2 = D2betDa2(n,x,a,b)
if (n==1)
  Dbeta2 = 0;
else
  D2alf1x = D2alfDa2(n-1,[2 x],a,b);
  D2bet1x = D2betDa2(n-1,[2 x],a,b);

  [Dalf1xa,Dalf1xb] = DalfDab(n-1,[2 x],a,b);
  
  Dbeta2 = (a-1)*D2alf1x + 2*Dalf1xa - (b-1)*D2bet1x;
end
return


function Dalfb2 = D2alfDb2(n,x,a,b)
if (n==1)
  Dalfb2 = 0;
else
  D2alf0x = D2alfDb2(n-1,[1 x],a,b);
  D2bet0x = D2betDb2(n-1,[1 x],a,b);

  [Dbet0xa,Dbet0xb] = DbetDab(n-1,[1 x],a,b);
  
  Dalfb2 = a*D2alf0x - 2*Dbet0xb  - b*D2bet0x;
end
return


function Dbetb2 = D2betDb2(n,x,a,b)
if (n==1)
  Dbetb2 = 0;
else
  D2alf1x = D2alfDb2(n-1,[2 x],a,b);
  D2bet1x = D2betDb2(n-1,[2 x],a,b);

  [Dbet1xa,Dbet1xb] = DbetDab(n-1,[2 x],a,b);
  
  Dbetb2 = (a-1)*D2alf1x - 2*Dbet1xb  - (b-1)*D2bet1x;
end
return



function Dalfab = D2alfDab(n,x,a,b)
if (n==1)
  Dalfab = 0;
else
  D2alf0x = D2alfDab(n-1,[1 x],a,b);
  D2bet0x = D2betDab(n-1,[1 x],a,b);
  
  [Dalf0xa,Dalf0xb] = DalfDab(n-1,[1 x],a,b);
  [Dbet0xa,Dbet0xb] = DbetDab(n-1,[1 x],a,b);

  Dalfab = a*D2alf0x + Dalf0xb - b*D2bet0x - Dbet0xa;
end
return


function Dbetab = D2betDab(n,x,a,b)
if (n==1)
  Dbetab = 0;
else
  D2alf1x = D2alfDab(n-1,[2 x],a,b);
  D2bet1x = D2betDab(n-1,[2 x],a,b);
  
  [Dalf1xa,Dalf1xb] = DalfDab(n-1,[2 x],a,b);
  [Dbet1xa,Dbet1xb] = DbetDab(n-1,[2 x],a,b);

  Dbetab = (a-1)*D2alf1x + Dalf1xb - (b-1)*D2bet1x - Dbet1xa;
end
return

function dY = d2Yda2(n,x,a,b)
Dalfa2 = D2alfDa2(n,x,a,b);
Dbeta2 = D2betDa2(n,x,a,b);
dY = Dalfa2 - Dbeta2;
return

function dY = d2Ydb2(n,x,a,b)
Dalfb2 = D2alfDb2(n,x,a,b);
Dbetb2 = D2betDb2(n,x,a,b);
dY = Dalfb2 - Dbetb2;
return

function dY = d2Ydab(n,x,a,b)
Dalfab = D2alfDab(n,x,a,b);
Dbetab = D2betDab(n,x,a,b);
dY = Dalfab - Dbetab;
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 = 7;
global nstate
nstate = 2;
global FN
FN = randFN(nstate,NN);
maxE = 0;
for t=1:1000000
  n = ceil(rand*(NN-2))+1;
  %n = 1;
  a = rand; 
  %b = rand;
  b = randinseg(0,a);
  %x = round(rand)+1;
  x = [];
  
  if 0
    %NDALFa = Ndiff('justalf',3,{n,x,a,b});
    %NDALFb = Ndiff('justalf',4,{n,x,a,b});
    %NDBETa = Ndiff('justbet',3,{n,x,a,b});
    %NDBETb = Ndiff('justbet',4,{n,x,a,b});
    %[Dalfa,Dalfb] = DalfDab(n,x,a,b);
    %[Dbeta,Dbetb] = DbetDab(n,x,a,b);
    %err1 = totdif(NDALFa,Dalfa);
    %err2 = totdif(NDALFb,Dalfb);
    %err3 = totdif(NDBETa,Dbeta);
    %err4 = totdif(NDBETb,Dbetb);
    
    %Dalfa2 = D2alfDa2(n,x,a,b);
    %Dbeta2 = D2betDa2(n,x,a,b);
    %Dalfb2 = D2alfDb2(n,x,a,b);
    %Dbetb2 = D2betDb2(n,x,a,b);
    %err1 = totdif(Dalfa2,Ndiff('dAda',3,{n,x,a,b}));
    %err2 = totdif(Dbeta2,Ndiff('dBda',3,{n,x,a,b}));
    %err3 = totdif(Dalfb2,Ndiff('dAdb',4,{n,x,a,b}));
    %err4 = totdif(Dbetb2,Ndiff('dBdb',4,{n,x,a,b}));
  
    %errAa = totdif(NDALFa,Dalfa);
    %errAb = totdif(NDALFb,Dalfb);
    %errBa = totdif(NDBETa,Dbeta);
    %errBb = totdif(NDBETb,Dbetb);
    %toterr = errAa + errAb + errBa + errBb
    %errAa2 = totdif(NDADa,Dalfa2);
    %errBa2 = totdif(NDBDa,Dbeta2);
    %toterr = errAa2 + errBa2;

    Dalfab = D2alfDab(n,x,a,b);
    Dbetab = D2betDab(n,x,a,b);
    err1 = totdif(Dalfab,Ndiff('dAda',4,{n,x,a,b}));
    err2 = totdif(Dbetab,Ndiff('dBda',4,{n,x,a,b}));
    %sanity1 = totdif(Ndiff('dAda',4,{n,x,a,b}),Ndiff('dAdb',3,{n,x,a,b}));
    %sanity2 = totdif(Ndiff('dBda',4,{n,x,a,b}),Ndiff('dBdb',3,{n,x,a,b}));
  
    %toterr = err1 + err2 + err3 + err4;
    toterr = err1 + err2;
    %toterr = sanity1 + sanity2;
    maxE = max(maxE,toterr);
    fprintf('maxE = %1.6f \n',maxE);
  end

  if 0
    da2 = d2Yda2(n,x,a,b);
    db2 = d2Ydb2(n,x,a,b);
    dab = d2Ydab(n,x,a,b);
    HESSY = [[da2 dab]; [dab db2]];
    da2 = d2Hda2(n,a,b);
    db2 = d2Hdb2(n,a,b);
    dab = d2Hdab(n,a,b);
    HESSH = [[da2 dab]; [dab db2]];
    
    [alf, bet] = alfbet(n,x,a,b);
    Y = alf-bet;
    Hn= H(n,a,b);
    
    HESS = HESSH - HESSY;
    HW3 = HESS3;
    err = totdif(HESS,HW3)
    keyboard
    
    if (Y>=0)
      HESS = HESSH - HESSY;
    else
      HESS = HESSH + HESSY;
    end
    
    sig = eig(HESS);
    g = all(sig>=0);
    if ~g
      'bad'
      keyboard
    end
  end
  
  
  if 0
    [alf, bet] = alfbet(n,x,a,b);
    Y = alf+bet;      %************* this is wrong!!!!
    Dalf = DalfDa(n,x,a,b);
    Dbet = DbetDa(n,x,a,b);
    DH   = dHda(n,a,b);
    
    if (Y>0)
      DF = DH - Dalf + Dbet;
    else
      DF = DH + Dalf - Dbet;
    end
    if (DF<0)
      save BAD FN n x a b
      'bad'
      keyboard
    end
  end
  
  if 1
    [alf0, bet0] = alfbet(n,[1 x],a,b);
    [alf1, bet1] = alfbet(n,[2 x],a,b);
    [A,B,C] = ABC(n,a,b);
    C1 = C + (a-b)^n;
    %if ( (a-b)^n < A )
    %  'muzar'
    %  keyboard
    %end
    
    d1 = alf0 - alf1;
    d2 = bet1 - bet0;
    d3 = alf1 - bet1;
    if ( abs(d1+d3)>C1 )
      'muzar13'
      keyboard
    end
    if ( abs(alf1-bet0)>C1 )
      'gam muzar'
      keyboard
    end

    
    
    %for t1=1:10000
      %a1 = rand;
      %b1 = randinseg(0,a1);
      a1 = 0; b1 = 0;
      g00 = ( abs(a1*d1 + b1*d2 + d3) <= C );
      a1 = 1; b1 = 0;
      g10 = ( abs(a1*d1 + b1*d2 + d3) <= C );
      a1 = 1; b1 = 1;
      g11 = ( abs(a1*d1 + b1*d2 + d3) <= C );
      if ~(g00 & g10 & g11)
        'nope'
        save BAD FN n x a b A B C C1 d1 d2 d3 a1 b1
        keyboard
      end
      
    %end
  end
  
  if 0
    FN = randFN(nstate,3);
    d1 = F([1 1 1 ])-F([2 1 2]);
    d2 = F([2 2 2])-F([1 2 1]);
    d3 = F([2 1 2])-F([2 2 2]);
    z=0:.001:1;
    f1 = 1+z+z.^2-d1*z-d3;
    f2 = 2+(1-z).^2-z-d1-d2*z-d3;
    f3 = 1-d1*z-d2*z-d3;
    g = (all(f1>=0) & all(f2>=0) & all(f3>=0));
    if ~g
      'bad, dude'
      keyboard
    end
  end
  
  %fprintf('good t = %d  \n',t);
  fprintf('good t = %d n = %d \n',t,n);
  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


function h = SC(AA)
% the "stricture coefficient" of AA
[n,n,T] = size(AA);
h = 0;
for t=1:T
  A = AA(:,:,t);
  for i=1:n
    for j=i+1:n
      x = A(:,i);
      y = A(:,j);
      h = max(h,max(abs(x-y)));
    end
  end
end
return