function lipsketch
global myopt
myopt = optimset('Diagnostics','off');
myopt = optimset('Display','off');


%dokeyb
lipL1check
%numlipconst
return

function lipL1check
global m n maxr Ainf Binf Ainf1 Ainf0 Binf1 Binf0 myopt
m = 3; n = 3;
[Ainf1,Binf1] = lipabl1(m,n);
[Ainf0,Binf0] = lipab(m,n);
maxr = 0;

while 1

  f = randn(m^n,1);
  haml1opt(f);

  if 0
    Ainf = Ainf1; Binf = Binf1;
    [g,gv] = maxg(f);
    Ainf = Ainf0; Binf = Binf0;
    [h,hv] = maxphi(f);
    
    %S = max(Psi(f,m,n),Psi(-f,m,n));
    %r = gv/S;
    r = gv/hv;
    if r > maxr
      global mg mf
      mg = g;
      mf = f;
      maxr = r;
    end
  end
  fprintf('                             maxr = %1.6f \n',maxr);
end
return

function y = haml1opt(f)
global m n maxr Ainf Binf Ainf1 Ainf0 Binf1 Binf0
global maxr myopt
Ainf = Ainf1; Binf = Binf1;
[g,gv] = maxg(f);
Ainf = Ainf0; Binf = Binf0;
[h,hv] = maxphi(f);
r = gv/hv;
if r > maxr
  global mg mf mh
  mg = g;
  mf = f;
  mh = h;
  maxr = r;
end
fprintf('                             maxr = %1.6f \n',maxr);
y = -r;
return


function numlipconst
global n A b maxr maxD
n = 2;

maxr = 0;

while 1

  A = randn(n,n);
  b = randn(n,1);
  x = rand(n,1);
  y = rand(n,1);
  
  maxD = maxgradL2(A,b);
  dxy = L2(x-y)+eps;

  Fx = F(A,b,x);
  Fy = F(A,b,y);

  %df = dF(A,b,x);

  r = abs(Fx-Fy) / dxy / maxD;

  if r > maxr
    global mA mb mx my
    mA = A; mb = b; mx = x; my = y;
    maxr = r;
    %fprintf('L1 maxr = %1.6f \n',maxr);
    fprintf('L2 maxr = %1.6f \n',maxr);
  end
  
end
return



function numlipconst0
global n A b maxr maxD
n = 1;

maxr = 0;

while 1

  Adim = n^2;
  bdim = n;
  xdim = n;
  ydim = n;
  fdim = Adim+bdim+xdim+ydim;
  %LB = zeros(fdim,1)-10;
  %UB = zeros(fdim,1)+10;
  %LB(Adim+bdim+1:end) = zeros(2*n,1);
  %UB(Adim+bdim+1:end) = ones(2*n,1);
  LB = zeros(2*n,1);
  UB = ones(2*n,1);


  
  A = randn(n,n);
  b = randn(n,1);
  x = rand(n,1);
  y = rand(n,1);
  
  maxD = maxgradL2(A,b);
  
  %f = [A(:); b; x; y;];
  f = [x;y];
  
  [f,v] = fmincon(@lipobj,f,[],[],[],[],LB,UB);

  fprintf('maxr = %1.6f \n',maxr);
end
return


function v = lipobj(xy)
global n maxr maxD A b
Adim = n^2;
bdim = n;
xdim = n;
ydim = n;

%A = reshape(Abxy(1:Adim),[n n]);
%b = reshape(Abxy(Adim+1:Adim+bdim),[n 1]);
%x = reshape(Abxy(Adim+bdim+1:Adim+bdim+xdim),[n 1]);
%y = reshape(Abxy(Adim+bdim+xdim+1:Adim+bdim+xdim+ydim),[n 1]);

x = xy(1:n);
y = xy(n+1:end);

Fx = F(A,b,x);
Fy = F(A,b,y);

dxy = L2(x-y)+eps;
%df = dF(A,b,x);

r = abs(Fx-Fy) / dxy / maxD;

if r > maxr
  global mA mb mx my
  mA = A; mb = b; mx = x; my = y;
  maxr = r;
  fprintf('maxr = %1.6f \n',maxr);
end

v = -r;
return


function f = F(A,b,x)
f = x'*A*x + b'*x;
return

function df = dF(A,b,x)
df = A*x + A'*x + b;
return


function d = L1(x)
d = sum(abs(x));
return

function d = L2(x,y)
d = sqrt(sum(x.^2));
return


function dokeyb
global m n maxr maxD
global mA mb mx my mf mg

keyboard
A=mA,b=mb,x=mx,y=my

Fx = F(A,b,x)
Fy = F(A,b,y)

dxy = L2(x-y)
%df = dF(A,b,x)

keyboard
return

function D = maxgradL1(A,b)
global n
xx = (list_all_ind(repmat(2,[1 n]))-1)';
D = 0;
for j=1:size(xx,2)
  x = xx(:,j);
  df = dF(A,b,x);
  D = max(D,L1(df));
end
return

function D = maxgradL2(A,b)
global n
xx = (list_all_ind(repmat(2,[1 n]))-1)';
D = 0;
for j=1:size(xx,2)
  x = xx(:,j);
  df = dF(A,b,x);
  D = max(D,L2(df));
end
return



function [g,gv] = maxg(f)
global Ainf Binf n myopt
if isempty(Ainf)
    error('Ainf not initialized')
end
[xp,fvalp] = linprog( f,Ainf,Binf,[],[],0*f,0*f+n,0,myopt);
[xn,fvaln] = linprog(-f,Ainf,Binf,[],[],0*f,0*f+n,0,myopt);
if abs(fvalp)>abs(fvaln)
  gv = fvalp;
  g = xp;
else
  gv = fvaln;
  g = xn;
end
gv = abs(gv);
return