function yashatest
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 4;
lenseq = 8;
global A
global C
%C = zeros(nstate,1); C(2) = 1;
C = pnormdim(ones(nstate,1),1);

e = .1;
A = zeros(4)+e;
A(1:2,1:2) = ones(2)-e;
A(3:4,3:4) = ones(2)-e;
A = pnormdim(A,1)


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

y = Y(1)

h = strict(A);
n = lenseq;
Hn = sum(h.^(0:n-1))

if (abs(y)<= Hn + 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 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 y = F(x)
%y = length(find(x~=1));
y = length(find(x>=2));
return