function tvtens
global m
%m = 4;

mrat = 0;

while 1
  
  m = round(randinseg(2,9));
  m = 2;
  
  mp = round(randinseg(2,9));
  mq = round(randinseg(2,9));
  mp = m; mq = m;
  p1 = pnormdim(rand(mp,1),1);
  p2 = pnormdim(rand(mp,1),1);
  q1 = pnormdim(rand(mq,1),1);
  q2 = pnormdim(rand(mq,1),1);
  
  p1q1 = tensprod(p1,q1);
  p2q2 = tensprod(p2,q2);

  lhs = TV(p1q1-p2q2);
  
  rhs = TV(p1-p2) + TV(q1-q2) - TV(p1-p2) * TV(q1-q2);
  %rhs = TV(p1-p2) * TV(q1-q2);  % NO! FALSE!!
  %rhs = min(TV(p1-p2),TV(q1-q2));  % NO! FALSE!!
  
  r = lhs / rhs;
  %if r>1.01
  %  keyboard
  %end
  
  mrat = max(r,mrat);
  fprintf('mrat = %1.9f \n',mrat);
end

return


function tvtens0
global m
%m = 4;

mrat = 0;

while 1
  
  m = round(randinseg(2,9));
  m = 3;
  
  mp = round(randinseg(2,9));
  mq = round(randinseg(2,9));
  mp = m; mq = m;
  p1 = pnormdim(rand(mp,1),1);
  p2 = pnormdim(rand(mp,1),1);
  q1 = pnormdim(rand(mq,1),1);
  q2 = pnormdim(rand(mq,1),1);
  
  A = pnormdim(rand(m,m),1);
  B = pnormdim(rand(m,m),1);
  
  p1q1 = tensprod(p1,q1);
  p2q2 = tensprod(p2,q2);
  AB = ABprod(A,B);

  lhs = TV(AB*(p1q1-p2q2));
  
  tp1 = AB*(p1q1-p2q2);
  tp2 = tensprod(A*(p1-p2),B*(q1-q2));
  %if totdif(tp1,tp2) > 1e-9
  %  'oops'
  %  keyboard
  %end

  thA = getTH(A);
  thB = getTH(B);
  %rhs = (thA + thB) * (TV(p1-p2) + TV(q1-q2));
  rhs = thA * TV(p1-p2) + thB * TV(q1-q2);
  %rhs = TV(tp2);
  
  r = lhs / rhs;
  
  if r>1.01
    keyboard
  end
  
  mrat = max(r,mrat);
  fprintf('mrat = %1.9f \n',mrat);
end

return


function tvtens1(m0)
global m mrat
m = m0;
mrat = 0;

[Aeq, Beq] = normab(m,2*m+4);
xdim = size(Aeq,2);
LB = zeros(xdim,1)+1e-12;
UB = LB*0+1;
while 1
  x = pnormdim(rand(m,2*m+4),1);
  x = x(:);
  
  [x,fv] = fmincon(@myopt,x,[],[],Aeq,Beq,LB,UB);
  
  x = reshape(x,[m 2*m+4]);

  fprintf('mrat = %1.9f \n',mrat);
end

return

function y = myopt(x)
global m mrat
x = reshape(x,[m 2*m+4]);
A = x(:,1:m);
B = x(:,m+1:2*m);
p1 = x(:,2*m+1);
p2 = x(:,2*m+2);
q1 = x(:,2*m+3);
q2 = x(:,2*m+4);

p1q1 = tensprod(p1,q1);
p2q2 = tensprod(p2,q2);
AB = ABprod(A,B);

lhs = TV(AB*(p1q1-p2q2));
  
thA = getTH(A);
thB = getTH(B);
rhs = thA * TV(p1-p2) + thB * TV(q1-q2);
  
r = lhs / rhs;
y = -r;

if r > mrat
  global mx
  mx = x;
  mrat = r;
end
fprintf('mrat = %1.9f \n',mrat);
return





function pq = tensprod(p,q)
mp = length(p);
mq = length(q);
pq = zeros(mp*mq,1);
for x = 1:mp
  for y = 1:mq
    xyind = coord2ind([x y],[mp mq]);
    pq(xyind) = p(x) * q(y);
  end
end
return


function AB = ABprod(A,B)
[mA,nA] = size(A);
[mB,nB] = size(B);
AB = [];
for x=1:nA
  for y = 1:nB
    abind = coord2ind([x y],[nA nB]);
    AB(:,abind) = tensprod(A(:,x),B(:,y));
  end
end
return

function d = TV(x)
d = .5 * sum(abs(x));
return



function h = getTH(A)
[n,n] = size(A);
h = 0;
for i=1:n
  for j=i+1:n
    x = A(:,i);
    y = A(:,j);
    d = TV(x-y);
    h = max(h,d);
  end
end
return