function refinestrict(N0,h)
global P Ltv TS FV kap N h0
N = N0
h0 = h

Ltv = totvarL(N);
kap = 1;

[P,fval] = minmaxtvdialate

return

TS = [];
FV = [];
for t=1:1000
  P = pnormdim(rand(N),1)
  
  TS(end+1) = totstrict(P);
  FV(end+1) = -fval;
  
  htil = totstrict(P)/(N-1);
  h = strict(P);
  if htil < h
    'nontriv'
    keyboard
  end
end
  
return

function [P,fval] = minmaxtvdialate
global N
P = eye(N);
P = P + randn(N)/1000;
P = pnormdim(P,1);
LB = zeros(N^2,1);
UB = ones(N^2,1);

Aeq = zeros(N,N^2);
for i=1:N
  for j=1:N
    Aeq(i,N*(i-1)+j) = 1;
  end
end
Beq = ones(N,1);

global myopt
[x,fval] = fmincon(@myobjminmax,P(:),[],[],Aeq,Beq,LB,UB,@mycon,myopt);
P = reshape(x,[N N]);
return


function ymax = myobjminmax(x)
global N h0
P = reshape(x,[N N]);
ymax = 0;
for t=1:10
  [x,y] = maxtvdialate(P);
  if y>ymax
    ymax = y;
  end
end
return


function [C, Ceq] = mycon(x)
global N h0
A = reshape(x,[N N]);
h = strict(A);
C = h0-h;
Ceq = 0;
return

function [C, Ceq] = trivcon(x)
C = 0;
Ceq = 0;
return



function [x,y] = maxtvdialate(P0)
global P Ltv TS FV kap N
P = P0;
x = randn(N,1); x = x-sum(x)/N;
Aeq = ones(1,N);
Beq = 0;
x = x/1000;
global myopt
[x,fval] = fmincon(@myobj,x,Ltv,kap*ones(size(Ltv,1),1),Aeq,Beq,[],[],[],myopt);
y = -fval;
return





function y = myobj(x)
global P Ltv
z = P*x;
y = -max(Ltv*z);
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 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;
return




function h = newstrict(P)
N = size(P,1);
h = 0;
for i=2:2^N
%for i=2^N
  setstr = dec2bin(i-1,N);
  ii1 = find(setstr=='1');
  Q = sum(P(ii1,:),1)
  for j=2:2^N-1
    setstr = dec2bin(j-1,N);
    jj0 = find(setstr=='0');
    jj1 = find(setstr=='1');
    d = abs(sum(Q(jj0)) - sum(Q(jj1)));
    if d > h
      h = d;
    end
  end
end

return