function y = getY(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),i);
    pXi = A(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_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 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),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 A
if i
  p = A(i,j,t);
else
  p = 1;
end
return




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

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


function r = F(x)
global nstate
if ~x(end), x=x(1:end-1); end;
lenseq = length(x);
global FN
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
%y = FF(ind);
try
r = FN(ind);
catch
  keyboard
end
return