function newcollapse
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 4;
lenseq = 6;


global FN
global NN
NN = lenseq;
%FN = randFN(nstate,NN);

%global ksi
%ksi = 1;

i0=1;

Xi = 1;  % fix X

global A C A0 C0

ALLSEQ_n =  list_all_ind(repmat(nstate,1,lenseq));
ALLSEQ1_n = list_all_ind(repmat(nstate,1,lenseq-i0+1));
ALLSEQ_2 =  list_all_ind(repmat(2,1,lenseq));
ALLSEQ1_2 = list_all_ind(repmat(2,1,lenseq-i0+1));




for t=1:2
  
  A = pnormdim(rand(nstate,nstate,lenseq-1),1);
  %A = pnormdim(rand(nstate,nstate,1),1);
  global A00
  %A = A00;
  %A = zeros(4);
  %A = eye(4);
  %A(1:2,1:2) = ones(2);
  %A(3:4,3:4) = ones(2);
  %A = pnormdim(A+eps,1);

  %A = repmat(A,[1 1 lenseq-1]);

  %C = pnormdim(rand(nstate,1),1);
  %C = pnormdim(ones(nstate,1),1);
  C = zeros(nstate,1); C(2) = 1;
  C = C/sum(C);

  A0 = A; C0 = C;
  h0 = strict(A0);

  ALLSEQ =  ALLSEQ_n;
  ALLSEQ1 = ALLSEQ1_n;

  Y1 = Y(Xi);
  %Y2 = Y_pr(Xi)    % they match
  %Y2 = Y_gam        % they match

  C = zeros(2,1); C(2) = 1;
  C = C/sum(C);

  ALLSEQ =  ALLSEQ_2;
  ALLSEQ1 = ALLSEQ1_2;

  A = collapseA;
  %C = collapseC;
  h1 = strict(A);  

  Y2 = Y(Xi);

  good = ((abs(Y1)<=abs(Y2)+1e-10) & (h1<=h0+1e-10));

  if ~good
    'BAAAAAAAAAAAAAAAAAAAAAAAAAAAD'
    keyboard
  end
  
end
return


function B = collapseA
global A

[n,n,T] = size(A);
B = zeros(2,2,T);
for t=1:T
  B(:,:,t) = collapse_t(A(:,:,t));  
end
B = B + realmin;
return

function B = collapse_t(A)
[n,n] = size(A);
p1 = A(:,1);
P1 = repmat(p1,[1 n]);
dd = sum(abs(A-P1),1);
[a,b] = max(dd);
B(1,1) = A(1,1);
%B(2,1) = 1-B(1,1);
B(2,1) = A(b,1);
B(1,2) = A(1,b);
B(2,2) = A(b,b);
B=pnormdim(B,1);
return


function r = F(x)
global nstate
%global lenseq
if ~x(end), x=x(1:end-1); end;
r = length(find(x~=1));
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));
    pxi = p_n(x(i+1),x (i),i);
    pXi = p_n(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_gam
global lenseq
n = lenseq;
gam = gamma(n,[]);
Y = sum(gam);
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 = p_n(Xix(i+1),Xix(i),i);
  b = p_n(Xi_1x(i+1),Xi_1x(i),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);
end
return

function p = P(x)
%global A
%global C
if isempty(x)
  p = 1;
else
  %p = C(x(1));
  p = p_0(x(1));
end
for j=2:length(x)
  %p = p * A(x(j),x(j-1)); % markov case
  p = p * p_n(x(j),x(j-1),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 = p_n(x(1),Xi(end));
end
for j=2:length(x)
  p = p * p_n(x(j),x(j-1),j-1); % markov case
end
return


function gam = gamma(n,x)
%global nstate
global A
[m,m,t] = size(A);
nstate = m;
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
  GAMkx = zeros(nstate,nstate);
  for k=1:nstate
    gamkx = gamma(n-1,[k x]);
    GAMkx(:,k) = gamkx;
    for k1=1:nstate
      gam(k) = gam(k) + p_n(k,k1,n-1) * gamkx(k1);
    end
  end
  GAMkx
  %if 
  %  'your wrong buddy'
  %  keyboard
  %end
end  
return

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

function p = p_n(i,j,n)
% primitive cond prob
global A
if i
  p = A(i,j,n);
else
  p = 1;
end
return


function [h,i1,i2] = strict0(A)
[n,n,T] = size(A);
h = 0;
i1 = 0;
i2 = 0;
for t=1:T
  for i=1
    for j=i+1:n
      x = A(:,i);
      y = A(:,j);
      d = .5*sum(abs(x-y));
      if d>h
        i1=i;
        i2=j;
        h=d;
      end
      
      % this stricture has the VERY NICE property
      % of not increasing after state collapse
      %d = max(d,max(abs(x-y)));
      % this is the OLD, BAD stricture def.
      % that's the problem that Sinai was so quick to point out
      % in my terminology, this stricture might increase after collapse
  end
  end
end
return