function mcd7
% 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];

test_alfbet
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 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 [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);
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 = [];
  h=a-b;
  
  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);
      H0 = sum(h.^nn);
      g00 = ( abs(alf0-bet0) <= H0+1e-10);
      g11 = ( abs(alf1-bet1) <= H0+1e-10);
      g10 = ( abs(alf1-bet0) <= H0+1+1e-10);
      g01 = ( abs(alf0-bet1) <= H0+1+1e-10);
      if ~(g00 & g11)
        'truly terrible'
        save BAD FN n x a b
        keyboard
      end
      if ~(g10 & g01)
        'kinda bad'
        save BAD FN n x a b
        keyboard
      end
    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