function psisketch(m0,n0)
global Ainf Binf m n mylinopt
m = m0; n=n0;
mylinopt = optimset('Diagnostics','off');
mylinopt = optimset('Display','off');

global K w ptil phi alf bet 

while 1
  K = randn(m^n,1);
  w = rand(n,1);
  w = w*0+1;
  
  % main claim
  [Ainf,Binf] = lipabw(m,n,w);
  g=tryfail(phinorm(K),psiwnorm(K,w),'phi<psi fails');
  if~g, keyboard, end
  
  w1 = w(1:n-1);
  %[x,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+sum(w),0,mylinopt);
  [x,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+sum(w1),0,mylinopt);
  phi = x';

  % explicit opt phi construction
  phat = optphi(K,w,m,n);
  
  lhs = phi*K;
  rhs = phat*K;
  
  g=tryfail(abs(lhs-rhs),0,'optphi fails');
  if~g, keyboard, end
  
  'good'
end
return


function phat = optphi(K,w,m,n)
global mylinopt
if n==1
  ksig = (K>0)';
  phat = w * ksig;
else
  w1 = w(1:n-1);
  [Ainf,Binf] = lipabw(m,n-1,w1);
  KSI = {};
  kk = [];
  for y=1:m
    Ky = getky(K,y,m,n);
    %phy = optphi(Ky,w1,m,n-1);
    [x,fvaln] = linprog(-Ky,Ainf,Binf,[],[],0*Ky,0*Ky+sum(w1),0,mylinopt);
    phy = x';
    
    KSI{y} = phy;
    kk(y) = sum(Ky);
  end
  %vv = optvv(KSI,kk,w(n),m,n-1);
  %phat = phat + max(vv);
  %vs = sum(vy) * (sum(K)>0);
  phat = zeros(1,m^n);
  for y=1:m
    phat = putky(phat,KSI{y},y,m,n);
  end
  %phat = phat + w(n)*(sum(K)>0);
end
return

function perturb(phat,K,w)
global Ainf Binf m n mylinopt
if sum(K) > 0
  % positives dominate
  
else
  
end
return

function vv = optvv(KSI,kk,w,m,n)
B = zeros(m);
%MN = repmat([m],[1 n]);
%mn = m^(n);
for y=1:m
  for z=setdiff(1:m,y)
    bmin = min(KSI{y}-KSI{z});
    B(y,z) = bmin;
  end
end
B = B + w;
kk = kk(:)';
Jplu = find(kk>=0);
Jmin = find(kk<0);
vv = zeros(1,m);
if sum(kk)>0
  % positives dominate
  vv(Jplu) = vv(Jplu) + w;
  for y=Jmin
    u = min(B(Jplu),y);
    vv(y) = w - u;
  end
else
  %negatives dominate
  for y=Jplu
    u = min(y,B(Jmin));
    vv(y) = u;
  end
end
return

function K = putky(K,Ky,y,m,n)
MN = repmat([m],[1 n]);
mn = m^(n);
MN1 = repmat([m],[1 n-1]);
mn1 = m^(n-1);
for ind = 1:mn1
  z = ind2coord(ind,MN1);
  zy = [z y];
  zyind = coord2ind(zy,MN);
  K(zyind) = Ky(ind);
end
return



function psisketch4(m0,n0)
global Ainf Binf m n mylinopt
m = m0; n=n0;
mylinopt = optimset('Diagnostics','off');
mylinopt = optimset('Display','off');

global K w ptil phi alf bet 

while 1
  K = randn(m^n,1);
  w = rand(n,1);
  
  if 1
    e = lemma43(K,w);
    g=tryfail(e,1e-9,'Psy[y] decomp fails');
    if~g, keyboard, end
  end
  
  % debugging
  %w(1)=0.01
  %w = n*w/sum(w); % this is important!!!
  %w = w*0+1;
  if 1
    % main claim
    LBK = -ones(size(K));
    UBK =  ones(size(K));
    LBf =   zeros(size(K));
    UBf =  n*ones(size(K));
    LB = [LBK];
    UB = [UBK];
    [Ainf,Binf] = lipabw(m,n,w);
    g=tryfail(phinorm(K),psiwnorm(K,w),'phi<psi fails');
    if~g, keyboard, end
  end

  [x,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+sum(w),0,mylinopt);
  phi = x';
  %if phi*K < 0
  %  K = -K;
  %end
  K1 = getK1(K,m,n);
  ksig = (K>0)';
  ptil = pl(phi - w(1)*ksig);
  
  
  [A,B0,B1,LB,UB0,UB1] = alfab(K,m,n,w);
  [x,fval] = linprog(-K,A,B0,[],[],LB,UB0,0,mylinopt);
  alf = x';
  %[x,fval] = linprog(-K,A,B1,[],[],LB,UB1,0,mylinopt);
  %bet = x';
  
  %b = polytopcheck(K,w,alf);
  %b = polytopcheck(K,w,ptil);
  
  %w1 = suff(w);
  %w1 = w1 * max((n-1)/(n-w(1)),1);
  %w1 = w1 * (n-1)/(n-w(1));

  
  if 1 % 1
    % lemma 4.7 attempt
    %failstr = 'lemma 4.7(1)'; 
    failstr = 'ptil'; 
    fprintf('%s\n',failstr);
    lhs = ptil * K;
    rhs = Psiw(K1,m,n-1,suff(w));
    %rhs = Psiw(K1,m,n-1,w1);
    g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
  end
  
  if 1 % 1
    failstr = 'lemma 4.7'; fprintf('%s\n',failstr);
    lhs = alf*K;
    rhs = Psiw(K1,m,n-1,suff(w));
    %rhs = Psiw(K1,m,n-1,w1);
    g=tryfail(lhs,rhs,'lemma 4.7 fails');
    if~g, keyboard, end
  end
  
  if 0 % 0
    % lemma 4.6 -- the "core" lemma
    failstr = 'lemma 4.6'; fprintf('%s\n',failstr);
    lhs = bet*K;
    w1 = w(1);
    %w1 = max(w(1),1);
    %w1 = w(2);
    %w1 = max(w);
    %rhs = w1*pl(sum(K)) + alf*K;
    rhs = w(1)*pl(sum(K))+Psiw(K1,m,n-1,suff(w));
    
    g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
  end
  
  if 1 % 1
    % Ky decomp
    failstr = 'Ky decomp'; fprintf('%s\n',failstr);
    wn = w(end);
    for y=1:m
      Ky = getky(K,y,m,n);
      Ky1 = getK1(Ky,m,n-1);
      bet = getky(alf,y,m,n);
      lhs = bet*Ky;
      rhs = wn*pl(sum(Ky)) + Psiw(Ky1,m,n-2,w(2:n-1));
      g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
    end
  end
  
  if 1
    % new Ky lemma
    failstr = 'new Ky lemma'; fprintf('%s\n',failstr);
    for ly=0:n-1
    %for ly=0:2
      ML  = repmat([m],[1 ly]);
      ml = m^ly;
      for yi = 1:ml
        y = ind2coord(yi,ML);
        [lhs, rhs] = gkycheck(K,w,y,m,n);
        g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
      end
    end
  end

  if 0 % 0
    % lemma 4.6
    % note -- abs is NO GOOD!
    % makes psi-decomp fail!
    failstr = 'lemma 4.6(0)'; fprintf('%s\n',failstr);
    wn = w(end);
    %wn = max(w);
    %wn = max(w(1),w(end));
    for y=1:m
      Ky = getky(K,y,m,n);
      Ky1 = getK1(Ky,m,n-1);
      bet = getky(alf,y,m,n);
      [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,w(1:n-1));
      [x,fval] = linprog(-Ky,A,B0,[],[],LB,UB0,0,mylinopt);
      gam1 = x';
      [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,w(1:n-1));
      [x,fval] = linprog(Ky,A,B0,[],[],LB,UB0,0,mylinopt);
      gam2 = x';
      gamK = max(abs(gam1*Ky),abs(gam2*Ky));
      %gamK = max(pl(gam1*Ky),pl(gam2*Ky));
      lhs = bet*Ky;
      %rhs = wn*pl(sum(Ky)) + gam*Ky;
      %rhs = wn*pl(sum(Ky)) + pl(gam1*Ky);
      %rhs = wn*pl(sum(Ky)) + abs(gam1*Ky);
      rhs = wn*pl(sum(Ky)) + gamK;
      %lhs = pl(lhs);      rhs = pl(rhs);
      g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
    end
  end

  if 0 % 0
    % lemma 4.6
    % makes psi-decomp fail!
    failstr = 'lemma 4.6(1)'; fprintf('%s\n',failstr);
    wn = w(end);
    for y=1:m
      Ky = getky(K,y,m,n);
      [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,w);
      [x,fval] = linprog(-Ky,A,B1,[],[],LB,UB1,0,mylinopt);
      bet1 = x';
      [x,fval] = linprog( Ky,A,B1,[],[],LB,UB1,0,mylinopt);
      bet2 = x';
      %lhs = max(abs(bet1*Ky),abs(bet2*Ky));
      lhs = abs(bet1*Ky);
      
      [x,fval] = linprog(-Ky,A,B0,[],[],LB,UB0,0,mylinopt);
      gam1 = x';
      [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,w(1:n-1));
      [x,fval] = linprog(Ky,A,B0,[],[],LB,UB0,0,mylinopt);
      gam2 = x';
      gamK = max(abs(gam1*Ky),abs(gam2*Ky));
      rhs = wn*pl(sum(Ky)) + gamK;
      g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
    end
  end
  
  
  
  if 0 % 0
    % the w(n+1) approach
    failstr = 'the w(n+1) approach'; fprintf('%s\n',failstr);
    %w1 = max(w(1),w(end));
    w1 = w(1);
    for y=1:m
      Ky = getky(K,y,m,n);
      %alfy = getky(alf,y,m,n);
      [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,w);
      [x,fval] = linprog(-Ky,A,B0,[],[],LB,UB0,0,mylinopt);
      alfy = x';
      [x,fval] = linprog(-Ky,A,B1,[],[],LB,UB1,0,mylinopt);
      bet = x';
      lhs = bet*Ky;
      %rhs = w1*pl(sum(Ky)) + alfy*Ky;
      Ky1 = getK1(Ky,m,n-1);
      rhs = w1*pl(sum(Ky)) + Psiw(Ky1,m,n-2,w(2:n-1));
      g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
    end
  end

  if 0
    % the w(2:n) approach
    failstr = 'the w(2:n) approach'; fprintf('%s\n',failstr);
    global iw
    for i=2:n
      iw = i;
      for y=1:m
        Ky = getkyi(K,y,m,n,i);
        %alfy = getky(alf,y,m,n);
        [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,w);
        [x,fval] = linprog(-Ky,A,B0,[],[],LB,UB0,0,mylinopt);
        alfy = x';
        [x,fval] = linprog(-Ky,A,B1,[],[],LB,UB1,0,mylinopt);
        bet = x';
        lhs = bet*Ky;
        rhs = w1*pl(sum(Ky)) + alfy*Ky;
        g=tryfail(lhs,rhs,failstr); if~g, keyboard, end
      end
    end
  end

    
  if 0
    alf = lemma46(K,w,bet);
    g=tryfail(bet*K,w(1)*pl(sum(K))+alf*K,'bad alf');
    if~g, keyboard, end
  end
  
  'good'
end
return

function [lhs, rhs] = gkycheck(K,w,y,m,n)
global mylinopt
l = n - length(y);
Ky = getkyg(K,y,m,n);
[A,B,LB,UB] = Awlt1(Ky,m,n,w,l,y);
w(end+1)=0;
Ky1 = getK1(Ky,m,l);
[x,fval] = linprog(-Ky,A,B,[],[],LB,UB,0,mylinopt);
alf = x';
lhs = alf * Ky;
ww = sum(w(l+1:end));
%rhs = w(l+1)*pl(sum(Ky)) + Psiw(Ky1,m,l-1,w(2:l));
rhs = ww*pl(sum(Ky)) + Psiw(Ky1,m,l-1,w(2:l));
g=tryfail(lhs,rhs,'gky'); if~g, keyboard, end
return

function Ky = getkyg(K,y,m,n)
MN = repmat([m],[1 n]);
Ky = [];
ly = length(y);
mn = m^n;
for ind = 1:mn
  z = ind2coord(ind,MN);
  if veq(z(n-ly+1:n),y)
    Ky(end+1) = K(ind);
  end
end
Ky = reshape(Ky,size(K(1:length(Ky))));
return


function Ky = getkyi(K,y,m,n,i)
MN = repmat([m],[1 n]);
Ky = [];
mn = m^n;
for ind = 1:mn
  z = ind2coord(ind,MN);
  if z(i)==y
    Ky(end+1) = K(ind);
  end
end
Ky = reshape(Ky,size(K(1:length(Ky))));
return


function alf = lemma46(K,w,bet)
bet = round(bet);
alf = pl(bet-w(1));
return



function psisketch00(m0,n0)
global m n
m = m0; n=n0;

K = randn(m^n,1);
w = randn(n,1);
%w = n*w/sum(w); % this is important!!!
%not important for the psi decomp
e = lemma43(K,w)

return

function e = lemma43(K,w)
global m n
S1 = Psiw(K,m,n,w);
S2 = altPsiw(K,m,n,w);
e = abs(S1-S2);
return


function psisketch0(m0,n0)
global Ainf Binf m n mylinopt maxr
%warning off
warning on

%m = 2;
%n = 3;
m = m0; n=n0;

mylinopt = optimset('Diagnostics','off');
mylinopt = optimset('Display','off');

K = randn(m^n,1);
phi = randlip(m,n)';

LBK = -ones(size(K));
UBK =  ones(size(K));

LBf =   zeros(size(K));
UBf =  n*ones(size(K));

LB = [LBK];
UB = [UBK];


global alpha

maxr = 0;
while 1
  
  alpha = rand(n,1);
  alpha = n*alpha/sum(alpha); % this is important!!!
  %alpha = 0*alpha+1;
  %alphasum = sum(alpha)
  [Ainf,Binf] = lipabw(m,n); 
  %keyboard
  
  x = randn(m^n,1);
  [x,fval] = fmincon(@myobj_ratio,x,[],[],[],[],LB,UB);
  r = abs(fval);
  if r > maxr
    global maxalf maxK
    maxalf = alpha;
    maxK = x;
    maxr = r;
  end
  fprintf('m=%d   n=%d                      maxr = %1.9f \n',m,n,maxr);
end
return


function y = myobj_ratio(x)
global m n maxr alpha
%phi = x(1:m^n)';
%K = x(m^n+1:end);
K = x;

yPhi = phinorm(K);
yPsi = psinorm(K,m,n)+realmin;

%[yPhi yPsi]

r = yPhi/yPsi;
if r > maxr
  global maxalf maxK
  maxalf = alpha;
  maxK = x;
  maxr = r;
  fprintf('m=%d   n=%d                      maxr = %1.9f \n',m,n,maxr);
end

y = -r;
return


function y = phinorm(K)
global Ainf Binf m n mylinopt w
W = sum(w);

[x,fvalp] = linprog( K,Ainf,Binf,[],[],0*K,0*K+W,0,mylinopt);
[x,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+W,0,mylinopt);
y = max(abs(fvalp),abs(fvaln));
return

function y = psinorm(K,m,n)
y = max(Psi(K,m,n),Psi(-K,m,n));
return

function y = psiwnorm(K,w)
global m n
y = max(Psiw(K,m,n,w),Psiw(-K,m,n,w));
return


function S = Psi(K,m,n)
global alpha
S = 0;
for t=1:n
  S = S + alpha(t)*sum(K.*(K>0));
  K = getK1(K,m,n-t+1);
end
%S = S / sum(alpha);
%S = S+pl(K);
return


function Ky = getky(K,y,m,n)
MN = repmat([m],[1 n]);
Ky = [];
mn = m^n;
for ind = 1:mn
  z = ind2coord(ind,MN);
  if z(end)==y
    Ky(end+1) = K(ind);
  end
end
Ky = reshape(Ky,size(K(1:length(Ky))));
return


function S = Psiw(K,m,n,w)
if n==0
  S = 0;
  return
end
%S = w(1)*sum(pl(K)) + Psiw(getK1(K,m,n),m,n-1,w(2:n));
S = w(1)*sum(pl(K)) + Psiw(getK1(K,m,n),m,n-1,suff(w));
return
%
S = 0;
for t=1:n
  S = S + w(t)*sum(K.*(K>0));
  K = getK1(K,m,n-t+1);
end
return

function S = altPsiw(K,m,n,w)
S = 0;
for y=1:m
  Ky = getky(K,y,m,n);
  S = S + Psiw(Ky,m,n-1,pref(w)) + w(n)*pl(sum(Ky));
  %S = S + Psiw(Ky,m,n-1,pref(w)) + w(n)*abs(sum(Ky));
  %%S = S + Psiw(Ky,m,n-1,suff(w)) + w(1)*pl(sum(Ky));
end
return


function S = Psi0(K,m,n)
S = 0;
for t=1:n
  S = S + sum(K.*(K>0));
  K = getK1(K,m,n-t+1);
end
return

function w = pref(w)
n = length(w)-1;
w = w(1:n);
%w = n*w/sum(w);
return

function w = suff(w)
n = length(w)-1;
w = w(2:end);
%w = n*w/sum(w);
return


function [A,B0,B1,LB,UB0,UB1] = alfab(K,m,n,w)
n1 = n-1;
MN  = repmat([m],[1 n ]);
MN1 = repmat([m],[1 n1]);
mn  = m^n;
mn1 = m^n1;

W = sum(w(1:n));

A=zeros(0,mn);
B=zeros(0,1);
rind = 0;
for ind1 = 1:mn
  for ind2 = ind1+1:mn
    x = ind2coord(ind1,MN);
    y = ind2coord(ind2,MN);
    dw = wham(x,y,w);
    d1 = wham(x,y,w*0+1);
    if (sign(K(ind1))==sign(K(ind2)))
      % signs equal -- ordinary Lip. cond
      % (SL2)
      for s = [-1,1]
        rind = rind+1;
        A(rind,ind1) =  s;
        A(rind,ind2) = -s;
        B0(rind,1) = dw;
        B1(rind,1) = dw;
      end
    elseif d1==1
    %else
      ii = find(x~=y);
      wi = w(ii);
      if (K(ind1)>0)
        pos = ind1;
        neg = ind2;
      else
        pos = ind2;
        neg = ind1;
      end
      rind = rind+1;
      A(rind,pos) =  1;
      A(rind,neg) = -1;
      %B(rind,1) = pl(wi-w(1));
      B0(rind,1) = (wi-w(1)); % need like this w/o pl
      B1(rind,1) = (wi-w(1));
      rind = rind+1;
      A(rind,neg) =  1;
      A(rind,pos) = -1;
      B0(rind,1) = wi + w(1);
      B1(rind,1) = wi + w(1);
    end
  end
end

LB = 0*K;
UB0 = LB + W   - (~~pl(K))*w(1);
%UB1 = LB + n+1 - (~~pl(K))*w(1);
%UB1 = LB + W + (~pl(K))*w(1);
if length(w)==n
  UB1 = []; % can't return a valid bet constraint
else
  UB1 = LB + W   - (~~pl(K))*w(1) + w(n+1);
end
%UB1 = UB0;
return


function alfabcheck(K,m,n,w,ksi)
n1 = n-1;
MN  = repmat([m],[1 n ]);
MN1 = repmat([m],[1 n1]);
mn  = m^n;
mn1 = m^n1;

for ind1 = 1:mn
  for ind2 = ind1+1:mn
    x = ind2coord(ind1,MN);
    y = ind2coord(ind2,MN);
    dw = wham(x,y,w);
    d1 = wham(x,y,w*0+1);
    if (sign(K(ind1))==sign(K(ind2)))
      % signs equal -- ordinary Lip. cond
      % (SL2)
      g = abs(ksi(ind1)-ksi(ind2)) < dw + 1e-9;
      if ~g
        'violates SL2'
        x,y
        %keyboard
      end
    elseif d1==1
      ii = find(x~=y);
      wi = w(ii);
      if (K(ind1)>0)
        pos = ind1;
        neg = ind2;
      else
        pos = ind2;
        neg = ind1;
      end
      pn = ksi(pos) - ksi(neg);
      g1 = -wi-w(1) < pn + 1e-9;
      g2 = pn < wi-w(1) + 1e-9;
      g = g1 & g2;
      if ~g
        'violates SL3'
        [ind1 ind2]
        x,y
        %keyboard
      end
    end
  end
end
return


function b = polytopcheck(K,w,alf)
global m n
for y = 1:m
  Ky = getky(K,y,m,n);
  Ay = getky(alf,y,m,n);
  %alfabcheck(Ky,m,n-1,suff(w),Ay)
  [A,B0,B1,LB,UB0,UB1] = alfab(Ky,m,n-1,pref(w));
  b = all(A*Ay' < B1+1e-9);
  if ~b
    'Lemma 4.4 fails'
    keyboard
  end
end
return

function h = wham(x,y,w)
ii = (x~=y);
w = w(1:length(x));
h = sum(ii(:).*w(:));
return


function good= tryfail(lhs,rhs,failstr)
global Ainf Binf mylinopt m n K w ptil phi alf bet 
good = (lhs <= rhs+1e-7);
if ~good
  failstr
  save PSIFAIL m n K w ptil phi alf bet 
  %keyboard
end
return



function [A,B,LB,UB] = Awlt0(Kt,m,n,w,l,t)
ML  = repmat([m],[1 l]);
ml = m^l;

W = sum(w);

A=zeros(0,ml);
B=zeros(0,1);
rind = 0;
for ind1 = 1:ml
  for ind2 = ind1+1:ml
    x = ind2coord(ind1,ML);
    y = ind2coord(ind2,ML);
    dw = wham(x,y,w(1:l));
    d1 = wham(x,y,w(1:l)*0+1);
    if (sign(Kt(ind1))==sign(Kt(ind2)))
      % signs equal -- ordinary Lip. cond
      % (SL2)
      for s = [-1,1]
        rind = rind+1;
        A(rind,ind1) =  s;
        A(rind,ind2) = -s;
        B(rind,1) = dw;
      end
    elseif d1==1
      ii = find(x~=y);
      wi = w(ii);
      if (Kt(ind1)>0)
        pos = ind1;
        neg = ind2;
      else
        pos = ind2;
        neg = ind1;
      end
      rind = rind+1;
      A(rind,pos) =  1;
      A(rind,neg) = -1;
      B(rind,1) = (wi-w(1));
      rind = rind+1;
      A(rind,neg) =  1;
      A(rind,pos) = -1;
      B(rind,1) = wi + w(1);
    end
  end
end

LB = 0*Kt;
UB = LB + W   - (~~pl(Kt))*w(1);
return



function [A,B,LB,UB] = Awlt1(Kt,m,n,w,l,t)
ML  = repmat([m],[1 l]);
ml = m^l;

W = sum(w);

A=zeros(0,ml);
B=zeros(0,1);
rind = 0;
for ind1 = 1:ml
  for ind2 = ind1+1:ml
    x = ind2coord(ind1,ML);
    y = ind2coord(ind2,ML);
    dw = wham(x,y,w(1:l));
    d1 = wham(x,y,w(1:l)*0+1);
    if (Kt(ind1)>0)
      pos = ind1;
      neg = ind2;
    else
      pos = ind2;
      neg = ind1;
    end
    rind = rind+1;
    A(rind,pos) =  1;
    A(rind,neg) = -1;
    B(rind,1) = (dw-w(1));
  end
end

LB = 0*Kt;
UB = LB + W   - (~~pl(Kt))*w(1);
return