function [maxH,maxpi] = maxhrand(P,m,n,N,H0)
if ~exist('H0','var')
  H0 = 0;
end
%global maxH maxpi
maxH = 0;
for i=1:N
%while 1
  pi = randperm(n);
  P1 = reorderP(P,m,n,pi);
  [hh,Ht,H] = gethhn(P1,m,n);
  Hpi = max(hh(:));
  if Hpi > maxH
    maxH = Hpi;
    maxpi = pi;
  end
  fprintf('maxH = %1.9f       H0 = %1.9f  \n',maxH,H0)
end
return