function mcd_naux
% writeup sanity check
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 3;
%lenseq = 2;
%global A
%global ec
global FN
global NN
NN=2;
maxR = 0;
global ksi

global AA
global C

for T=1:10000
  FN = randFN(nstate,NN);
  ak = pnormdim(rand(nstate,1),1);
  bk = pnormdim(rand(nstate,1),1);
  h = max(abs(ak-bk));
  s = 0;
  for k=1:nstate
    s = s+ak(k)*F([1 k])-bk(k)*F([2 k]);
  end
  g = (abs(s)<=1+h+1e-10);
  if ~g
    'problem!!!!'
    keyboard
  end
end
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 AA
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 = AA(x(i+1),x(i),i);
    pXi = AA(x(i+1),Xi(i),i);
    p = p / pxi * pXi;
  end
  f = F([Xi x(i+1:end)]);
  Ef = Ef + f*p;
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),i);
  b = pcp(Xi_1x(i+1),Xi_1x(i),i);
  p = PCOND(x,Xi(1:i-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 p = P(x)
global AA
global C
if isempty(x)
  p = 1;
else
  p = C(x(1));
end
for j=2:length(x)
  p = p * AA(x(j),x(j-1),j-1); % markov case
  %p = p * C(x(j)); % independent case
end
return

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

function p = pcp(i,j,t)
% primitive cond prob
global AA
if i
  p = AA(i,j,t);
else
  p = 1;
end
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
C = pnormdim(rand(nstate,1),1);
%AA = initA(sc);
%return
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;
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 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

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 s=test12(sc)
global nstate
global lenseq
lenseq = 2;
global FN;
%FN = randFN(nstate,lenseq);
%FN = randFN(nstate,NN);
global AA
global C
X1 = 1;
%X1 = 0;
%x1 = 2;
for t=1:1
  %[AA,C] = initAAC(sc);
  %[A,C] = initAC(erg);
  %h = SC(AA);
  s = 0;
  MSTR = '';
  for x1=1:nstate
    for x2=1:nstate
      s=s+ p_0(x1)*(p_n(x2,X1,1)*F([X1 x2]) - p_n(x2,x1,1)*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])));
      %fprintf('%s \n',str);
    end
  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 p=p_0(k)
global C
p=C(k);
return

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


function Y=Y_gam
global lenseq
n = lenseq;
gam = gamma(n,[]);
%Y = gam(1)-sum(gam(2:end));
Y=sum(gam);

% p_0(1)*(F(1)-F(1)) + p_0(2)*(F(1)-F(2)) + p_0(3)*(F(1)-F(3)) 

return


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


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






function Y=Y_ab
global lenseq
% nstate = 2
% this is really $Y_1^{(n)}$
n = lenseq;

%[alf_n, bet_n] = alfbet(n,[]);
%Y = sum(alf_n - bet_n);
[alf, bet] = alfbet0(n,[]);
Y = alf + bet;
%[alf, bet] = alfbet(1,[]);
return


function [alf, bet] = alfbet(n,x)
global nstate
alf = 0;
bet = zeros(nstate,1);
if (n==1)
  alf = p_0(1)*F([1 x]);
  for k=2:nstate
    bet(k) = p_0(k)*F([k x]);
  end
else
  [alf0x, bet0x] = alfbet(n-1,[1 x]);
  alf = p_n(1,1,n-1) * alf0x;
  for k=1:nstate
    [alfkx, betkx] = alfbet(n-1,[k x]);
    for k1=2:nstate
        alf    = alf                          - p_n(1,k1,n-1)*bet0x(k1);
        bet(k) = bet(k) + p_n(k1,1,n-1)*alfkx - p_n(k,k1,n-1)*betkx(k1);
    end
  end
  bet = -bet;
end
return


function [alf, bet] = alfbet00(n,x)
global nstate
alf = zeros(nstate,1);
bet = zeros(nstate,1);
if (n==1)
  for k=1:nstate
    alf(k) = p_0(k)*F([1 x]);
    bet(k) = p_0(k)*F([k x]);
  end
else
  [alf0x, bet0x] = alfbet(n-1,[1 x]);
  for k=1:nstate
    [alfkx, betkx] = alfbet(n-1,[k x]);
    for k1=1:nstate
        for k2=1:nstate
            alf(k) = alf(k) + p_n(1,1,n-1)*alf0x(k1) - p_n(1,k2,n-1)*bet0x(k2);
            bet(k) = bet(k) + p_n(k,1,n-1)*alfkx(k1) - p_n(k,k2,n-1)*betkx(k2);
        end
    end
  end
  bet = -bet;
end
return



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