function [q,v] = martonthm2
global m n Ainf Binf
global myopt
myopt = optimset('Diagnostics','off');
myopt = optimset('Display','off');
warning off

m = 2;
n = 3;

mn = m^n;

keyboard

q = rand(mn,1); q = q/sum(q);
pp = allconddistr(q,m,n);
%phi = fillgmask(randmask(m,n),m);

K = phicoeff(pp,m,n);

% my Y decomp works
%y0 = Y(q,phi)
%y1 = phi*K

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

[q,fval] = fmincon(@myobj,q,ones(size(q')),1,[],[],0*q+1e-8,0*q+1);
v = 1/fval;

return

G = totgam(q,m,n)

[x,fval] = linprog(K,Ainf,Binf,[],[],0*K,0*K+n);
v = abs(fval)
good = v < G + 1e-7
return


function y = myobj(q)
global m n Ainf Binf myopt
pp = allconddistr(q,m,n);
K = phicoeff(pp,m,n);
G = totgam(q,m,n);
[x,fval] = linprog(K,Ainf,Binf,[],[],0*K,0*K+n,[],myopt);
v = abs(fval);
%if isempty(v)
%  keyboard
%end
y = G/v;
return


function y = Y(q,phi)
% crudest version possible
global m n
n1 = n-1;
mn = m^n;
mn1 = m^n1;
MN = repmat([m],[1 n]);
MN1 = repmat([m],[1 n1]);
Ephi = phi*q;
Ephi1 = 0;
p0 = 0;
% get marginal for P(x(1)=1)
for xind = 1:mn1
  x1 = ind2coord(xind,MN1);
  x = [1 x1];
  xind = coord2ind(x,MN);
  p0 = p0 + q(xind);
  Ephi1 = Ephi1 + phi(xind)*q(xind);
end
Ephi1 = Ephi1/p0;
y = Ephi1 - Ephi;
return


function G = totgam(q,m,n)
gg = 1;
for k=1:n-1
  g = gamk(q,m,n,k);
  gg(end+1) = g;
end
G = sum(cumprod(gg));
return

function g = gamk(q,m,n,k)
Mk = repmat([m],[1 k]);
mk = m^k;
g = 0;
for xind=1:mk
  x = ind2coord(xind,Mk);
  px = conddistr(q,m,n,x);
  for yind=xind+1:mk
    y = ind2coord(yind,Mk);
    py = conddistr(q,m,n,y);
    g = max(g,.5*sum(abs(px-py)));
  end
end
return


function p = conddistr(q,m,n,x)
k = length(x);
% we need 1 <= k < n
p = zeros(m,1);
px = margprob(q,m,n,x);
for i=1:m
  xi = [x i];
  p(i) = margprob(q,m,n,xi)/(px+realmin);
end
return

function p = margprob(q,m,n,x)
k = length(x);
% we need 1 <= k < n
p = 0;
nk = n-k;
MNk = repmat([m],[1 nk]);
MN = repmat([m],[1 n]);
mnk = m^nk;
for yind = 1:mnk
  y = ind2coord(yind,MNk);
  xy = [x y];
  xind = coord2ind(xy,MN);
  p = p + q(xind);
end
return



function pp = allconddistr(q,m,n)
pp = {};
for k=0:n-1
  mk = m^k;
  Mk = repmat([m],[1 k]);
  pk = zeros(m,mk);
  for xind = 1:mk
    x = ind2coord(xind,Mk);
    pk(:,xind) = conddistr(q,m,n,x);
  end
  pp{k+1} = pk;
end
return


function K = phicoeff(pp,m,n)
dims = repmat([m],[1 n]);
K = zeros(m^n,1);

for xind=1:m^n
  x = ind2coord(xind,dims);
  if x(1)==1
    K(xind) = (1-p0(pp,1))*ptil(pp,m,n,x);
  else
    K(xind) = -p0(pp,x(1))*ptil(pp,m,n,x);
  end
end
return

function p = p0(pp,i)
p = pp{1}(i);
return

function p = pcond(pp,m,n,x,y)
% p(y|x)
k = length(x);
dims = repmat([m],[1 k]);
xind = coord2ind(x,dims);
p = pp{k+1}(y,xind);
return


function p = ptil(pp,m,n,x)
p = 1;
for j=2:length(x)
  %p = p * A(x(j),x(j-1),j-1);
  p = p * pcond(pp,m,n,x(1:j-1),x(j));
end
return