function newcollapse1
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 3;
lenseq = 2;


global FN
global NN
NN = lenseq;
i0=1;

Xi = 1;  % fix X

global A C
global AApairs CC BAR
AApairs = {};
CC = {};
BAR = [];

for t=1:10
%for t=1

ALLSEQ =  list_all_ind(repmat(nstate,1,lenseq));
ALLSEQ1 = list_all_ind(repmat(nstate,1,lenseq-i0+1));

  
A = pnormdim(rand(nstate,nstate,lenseq-1),1);
%A = pnormdim(rand(nstate,nstate,1),1);
global 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]);
%A = A00;


%C = pnormdim(rand(nstate,1),1);
%C = pnormdim(ones(nstate,1),1);
C = zeros(nstate,1); C(2) = 1;
C = C/sum(C);
global A0 C0
A0 = A; C0 = C;
%C = zeros(nstate,1); C(2)=1;
global h0
h0 = strict(A0);

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

%A = pnormdim(rand(2,2,lenseq-1),1);
A = collapseA;
C = pnormdim(rand(2,1),1);
%C = zeros(2,1); C(2) = 1;
C = C/sum(C);

global CURE
d = .1;
%global PI PI0 PI1
%PI = pnormdim(rand(nstate,lenseq-1),1);
%PI = pnormdim(ones(nstate,lenseq-1),1);
%PI0 = PI;
%PI1 = pnormdim(ones(2,1),1);
ee = [];
ALLSEQ =  list_all_ind(repmat(2,1,lenseq));
ALLSEQ1 = list_all_ind(repmat(2,1,lenseq-i0+1));

h = strict(A0);
[A,C,fval] = mylp(h);
if abs(fval)+1e-10>abs(Y1)
  'success!'
  AApairs{t,1} = A0;
  AApairs{t,2} = A;
  CC{t} = C;
  %break;
end
end
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] = strict(A)
[n,n,T] = size(A);
h = 0;
i1 = 0;
i2 = 0;
for t=1:T
  %for i=1:n
  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
    end
  end
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