function [A,fval] = minmaxprdev(m0,n0,Hh)
global m n H0 h0
T = n-1;
m=m0; n=n0; 
%H0 = Hh; h = [];
h0 = Hh; H = [];

adim = m^2*T;

LB = zeros(adim,1);
UB = ones(adim,1);

% normalization constraints
Aeq = zeros(m*T,adim);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,1);

C = ones(m,1)/m;

if flip
  A = repmat(eye(m),[1 1 T]);
  A = A + rand(m,m,T)/10000;
  A = pnormdim(A,1);
else
  % initialize via minmaxvar
  [A,fval] = minmaxvar_hcon(m,n,Hh);
end

global myopt
[x,fval] = fmincon(@myobj,A(:),[],[],Aeq,Beq,LB,UB,@myconh,myopt);
A = reshape(x,[m m T]);
return


function [C, Ceq] = myconh(x)
global m n h0
T = n-1;
A = reshape(x,[m m T]);
Ceq = 0;
rind=0;
C = [];
for t=1:T
  for i1=1:m
    for i2=i1+1:m
      rind = rind+1;
      p1 = A(:,i1,t);
      p2 = A(:,i2,t);
      C(rind,1) = h0 - sum(abs(p1-p2))/2;
    end
  end
end
return


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


function [C, Ceq] = myconH(x)
global m n H0
T = n-1;
A = reshape(x,[m m T]);
H = getH(A);
C = H0-H;
Ceq = 0;
return


function y = myobj(x)
global m n H0
T = n-1;
A = reshape(x,[m m T]);
C = ones(m,1)/m;

[y,junk] = maxprobdev(A,C);
return