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


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

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

Xi = 1;  % fix X

for t=1:100

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


Y_old = Y(Xi)

A0=A;
A = optcollapse(A)
B = zeros(nstate);
B(1:2,1:2) = A;
A=B+realmin;

Y_new = Y(Xi)

if (abs(Y_old)<=abs(Y_new)+1e-10)
    'good'
    g=1;
    %break;
else
    g=0;
    'bad'
        keyboard
end

end

return


function B = optcollapse(A)
n = size(A,1);
B = pnormdim(ones(2),1);
maxh = 0;
for k=1:n-1
    combs = generate_combs(n,k);
    for c=1:size(combs,1)
        comb = combs(c,:);
        ii = find(comb);
        jj = find(~comb);
        ppi = mean(A(:,ii),2);
        ppj = mean(A(:,jj),2);
        p00 = sum(ppi(ii));
        Pii = sum(ppi(ii));
        Pij = sum(ppi(jj));
        Pji = sum(ppj(ii));
        Pjj = sum(ppj(jj));
        b = [[Pii Pij];[Pji Pjj]];
        h = strict(b);
        if h > maxh
            B = b;
        end
    end
end
return

function XX=collaps(XX)
global nstate
%for n=3:nstate
for j=1:size(XX,2)
    x = XX(:,j);
    ii = find(x~=1);
    x(ii) = 2*ones(size(ii));
    XX(:,j) = x;
end
return

function y = h(x)
y = (x==1) + 2*(x~=1);
return

function y = evalF(FN,x)
global nstate
lenseq = length(x);
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
y = FF(ind);
return


function r = F(x)
global nstate
%global lenseq
if ~x(end), x=x(1:end-1); end;
r = length(find(x~=1));
return
lenseq = length(x);
global FN
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
r = FF(ind);
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(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_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));
  b = pcp(Xi_1x(i+1),Xi_1x(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);
  %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)); % 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 = pcp(x(1),Xi(end));
end
for j=2:length(x)
  p = p * pcp(x(j),x(j-1)); % markov case
end
return

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



function FN = randFN(nstate,NN,d)
if ~exist('d','var')
  d=1;
end
FN = {};
for j=1:NN
  FN{j} = randF(nstate,j,d);
end
return

function F = randF(nstate,lenseq,d)
% F is a random function $\X^n \to [0,1]$ that sat. ibd w/c=1
X = list_all_ind(repmat(nstate,1,lenseq));
K = 10;
%K = 1;
%K = lenseq;
%F(1) = 0;
PROGstart(size(X,1),'generating randF...');
for i=1:size(X,1)
%for i=2:size(X,1)
  good = 0;
  while ~good
    F(i) = K*rand;
    %F(i) = round(K*rand);
    good = 1;
    for j = 1:i-1
      good = goodij(X,F,i,j,d);
      if ~good, break, end;
    end
  end
  %fprintf('randF: done %d of %d \n',i,size(X,1));
  PROGupdate(i);
end
PROGend;
return

function b = goodij(X,F,i,j,d)
x = X(i,:);
y = X(j,:);
if (length(find(x~=y))==1)
  b = (abs(F(i)-F(j))<=d);
  %b = (abs(F(i)-F(j))==1);
else
  b = 1;
end
return



function r = randinseg(a,b)
r = rand*(b-a) + a;
return


function [h,i1,i2] = strict0(A)
[n,n] = size(A);
h = 0;
i1 = 0;
i2 = 0;
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
return