function sancheck
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
nstate = 4;
lenseq = 4;
i0=1;

m = nstate;

global A C
global CURE
d = .01;
nn = 0:lenseq-1;


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

for t=1:T
  %for t=1
  
  A = pnormdim(rand(nstate,nstate,lenseq-1),1);  
  h = strict(A);
  C = pnormdim(rand(nstate,1),1);
  
  while 1
    dA = Ndiff('ER',2,{0,A,C});
    dA = reshape(dA,[m m lenseq-1]);
    %dC = Ndiff('ER',3,{0,A,C});
    %dC = reshape(dC,[m 1]);
    e = CURE;
    
    % check main condition
    H = sum(h.^nn);
    g = (sqrt(e)<=H+1e-10);
    if ~g
      'REEEEEEEEEEEEAL PROBLEM!!!!'
      keyboard
    end
    
    % perhaps gradient has stopped changing?
    %if sum(sum(abs(dA)))+sum(abs(dC)) < 1e-7
    if sum(sum(abs(dA))) < 1e-7
      brstr = 'converged';
      break
    end
    
    A = A + d*dA;
    A = A.*(A>=0);
    A = pnormdim(A+realmin,1);
    h1 = strict(A);
    if h1>h
      % we've violated the strict. constraint
      brstr = 'constraint';
      break
    end

    %C = C + d*dC;
    %C = C.*(C>=0);
    %C = pnormdim(C,1);
  end

  fprintf('verified %5d out of %5d breaktype = %s\n',t,T,brstr);
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 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