function direct_bound
global nstate
nstate = 2;
global lenseq
lenseq = 2;

global A C
A = rand(nstate,nstate);
C = zeros(nstate,1); C(2) = 1;

global ALLSEQ
ALLSEQ =  list_all_ind(repmat(nstate,1,lenseq));

N = 4;
n = 2; 
G = 0;
Xk = zeros(nstate,N);
for k=1:nstate
  x = randx(N);
  Xk(k,:) = x;
  gam = gamma(n,x);
  G = G + gam(k);
end
Xk
B = bound(A,Xk,n);
if (abs(G) <= B + 1e-10)
  'good'
else
  'bad'
  keyboard
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 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
  f = F([Xi x(i+1:end)]);
  Ef = Ef + f*p;
end
return


function Y=Y_gam
global lenseq
n = lenseq;
gam = gamma(n,[]);
Y=sum(gam);
return


function B = bound(A,Xk,n)
global nstate
h = strict(A);
R = 0;
for k1=1:nstate
  for k2=k1+1:nstate
    R = R + rho(Xk(k1,:),Xk(k2,:));
  end
end
Hn = sum(h.^(0:n-1));
B = Hn + h^(n-1)*R;
return

function x = randx(N)
global nstate
x = ceil(nstate*rand(1,N));
return

function r=rho(x,y)
x=colps(x);
y=colps(y);
r = length(find(x~=y));
return

function x=colps(x)
ii = find(x~=1);
x(ii)=2;
return

%function gam = gamma(n,x,MC,FN)
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 y = F(x)
y = length(find(x~=1));
return

function p = p_0(x)
global A
global C
p = C(x);
return

function p = p_n(k,k1,n)
global A
global C
p = A(k,k1);
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