function sketchpadtighten
global nstate
global lenseq
global ALLSEQ
global ALLSEQ1
global XXi
global Ainf Binf H
nstate = 2;
lenseq = 4;
global A
global C
global m n i0
m = nstate;
n = lenseq;

C = pnormdim(rand(m,1),1);
C = C*0; C(2)=1;


T = n-1;

adim = m^2*T;
fdim = m^n;
xdim = adim + fdim;

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


% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
%LB(adim+1:end) = LB(adim+1:end) - 1e6;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + n;  % no box constraints on phi


% normalization constraints
Aeq = zeros(m*T,xdim);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,1);


h = rand;
[Aina,Bina] = constrict(m,h,T);

[Ainf,Binf] = lipab(m,n);

Ain = assembleblock({Aina,Ainf});
Bin = [Bina;Binf];

V = 0;

while 1

  A0 = pnormdim(ones(m,m,T),1);
  A0 = A0+randn(size(A0))*1e-2;
  A0 = pnormdim(A0,1);
  phi0 = randlip(m,n);
  
  x0 = [A0(:); phi0(:)];
  [x,fval] = fmincon(@myobj1,x0,Ain,Bin,Aeq,Beq,LB,UB);
  v = -fval;
  
  if v > V
    V = v;
  end
  curmax = V
  if V > 1.2
    'bad bound'
    keyboard
  end
end

return

gdim = m^n*m;
adim = m^2;
xdim = gdim + adim;


global Ainf Binf H
[Ainf,Binf] = lipab(m,n);
i0=1;
ALLSEQ =  list_all_ind(repmat(nstate,1,lenseq));
ALLSEQ1 = list_all_ind(repmat(nstate,1,lenseq-i0+1));
XXi = list_all_ind(repmat(nstate,1,i0));


for tt=1:100

global FN
FN = randlip(m,n);
A = pnormdim(rand(m,n,n-1),1);
C = pnormdim(rand(m,1),1);
C = C*0; C(2)=1;
H = getH(A,i0);

[phi,fval] = fmincon(@myobj,FN(:),[],[],Ainf,Binf);
v = -fval

return
if 0
  %for i0=1:lenseq
  for i0=1
    %for i0=1:n-1
    %lenseq = i0+1;
    YY = [];
    %for x1 = 1:nstate
    %for x1 = 1
    %for j=1
    for j=1:size(XXi,1)
      Xi = XXi(j,:);
      YY(j) = getY(Xi);
      %YY(j) = Y_pr(Xi);
      %AB(j) = Y_ab;
    end
  %if (totdif(YY,AB)>1e-10)
  %if (abs(YY)>abs(AB)+1e-10)
  %'dont match up!!!!'
  %'bad bound!!!'
  %keyboard
  %end
  end
  YY
  d = max(YY) - min(YY);
  g = d < H
  if ~g
    'nope'
    keyboard
  end
end
end
return


function y = myobj(x)
global m n FN H
global ALLSEQ
global ALLSEQ1
global XXi
FN = x;
YY = [];
for j=1:size(XXi,1)
  Xi = XXi(j,:);
  YY(j) = Y(Xi);
  %YY(j) = Y_pr(Xi);
  %AB(j) = Y_ab;
end
d = max(YY) - min(YY);
y = -d/H;
return


function y = myobjmin(x)
global m n FN H
global ALLSEQ
global ALLSEQ1
global XXi
FN = x;
YY = [];
for j=1:size(XXi,1)
  Xi = XXi(j,:);
  %YY(j) = Y(Xi);
  YY(j) = Y_pr(Xi);
  %AB(j) = Y_ab;
end
y = min(YY);
return


function y = myobjmax(x)
global m n FN H
global ALLSEQ
global ALLSEQ1
global XXi
FN = x;
YY = [];
for j=1:size(XXi,1)
  Xi = XXi(j,:);
  %YY(j) = Y(Xi);
  YY(j) = Y_pr(Xi);
  %AB(j) = Y_ab;
end
y = max(YY);
return



function y = myobj1(x)
global m n FN H
global ALLSEQ
global ALLSEQ1
global XXi
global A i0
T = n-1;
adim = m^2*T;
fdim = m^n;
xdim = adim + fdim;

A = reshape(x(1:adim),[m m T]);
FN = x(adim+1:end);

YY = [];
for j=1:size(XXi,1)
  Xi = XXi(j,:);
  try
  %YY(j) = Y_pr(Xi);
  YY(j) = getY(Xi);
  catch
    keyboard
  end
  %AB(j) = Y_ab;
end
H = getH(A,i0);
d = max(YY) - min(YY);
y = -d/H;
return

function y = evalF(FN,x)
global nstate
%lenseq = length(x);
%FF = FN{lenseq};


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


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+eps) * pXi;
  end
  %pjoin = P(X);
  %pmarg = P(X(1:i));
  %pcond = pjoin/pmarg;
  %p = P(X);
  f = F([Xi x(i+1:end)]);
  Ef = Ef + f*p;
  %Ef = Ef + f*pcond;
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),i);
  b = pcp(Xi_1x(i+1),Xi_1x(i),i);
  p = PCOND(x,Xi(1:i-1));
  p = p / (b+eps);

  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
end
return

function p = PCOND(x,Xi)
global C
if isempty(Xi)
  p = C(x(1));
else
  p = pcp(x(1),Xi(end),length(Xi));
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 phi = randlip(m,n)
global Ainf Binf
%[A,B] = lipab(m,n);
K = genK(m,n);
[x,fval] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n);
phi = round(x');
return



function [Aina,Bina] = constrict(m,h,T)
% stricture constraints

if ~exist('T','var')
  T = 1;
end

SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina = zeros(size(SIGNS,1)*m*(m-1)/2*T,m^2*T);
rind = 0;

for t=1:T
  for i1=1:m
  %for i1=1
    for i2=i1+1:m
      for s = 1:size(SIGNS,1)
        rind = rind+1;
        ss = SIGNS(s,:);
        for j=1:m
          ind1 = coord2ind([j i1 t],[m m T]);
          ind2 = coord2ind([j i2 t],[m m T]);
          Aina(rind,ind1) =  ss(j);
          Aina(rind,ind2) = -ss(j);
        end
      end
    end
  end
end
Bina = ones(size(Aina,1),1)*h;
Bina = Bina * 2;
%keyboard


return