function [gam,fval] = maxgfeat(A)
% exploring |Ag|
global nstate
nstate=size(A,1);
m = nstate;
h = strict(A);
C = 1+h;
C1 = C+h^2;

adim = 0;
gdim = m^2;
xdim = gdim;


rind=0;
Aeq = zeros(0,xdim);
Beq = zeros(0,1);
if 0
  % a quick'n'dirty hack to fix gam
  global GG0
  for i=1:m
    for j=1:m
      %if (i>=1) | (j>=1)
      %if (j>=3) & (i>=3)
      %if j>=2;
      %if (j>=2) | (i>=2)
      if 1
        rind = rind+1;
        cind = adim + coord2ind([i j],[m m]);
        Aeq(rind,cind) = 1;
        Beq(rind,1) = GG0(i,j);
      end
    end
  end
end

% gamma sum constraint
for i=1:m
  for s=[-1 1]
    rind = rind+1;
    for j=1:m
      cind = coord2ind([i j],[m m]);
      Aing(rind,cind) = s;
    end
    Bing(rind,1) = C;
  end
end

% gamma Lip constraint:
d=1;
for j=1:m
  for i1=1:m
    for i2=i1+1:m
      ind1 = coord2ind([i1 j],[m m]);
      ind2 = coord2ind([i2 j],[m m]);
      for s=[-1 1]
        rind = rind+1;
        Aing(rind,ind1) =  s;
        Aing(rind,ind2) = -s;
        Bing(rind,1) = 2*d;
        %Bing(rind,1) = .22;
        %Bing(rind,1) = d*h;
        %0.1124    0.1000    0.2257
        %switch j
        % case 1
        %  Bing(rind,1) = 0.1124;
        % case 2
        %  Bing(rind,1) = 0.1000;
        % case 3
        %  Bing(rind,1) = 0.2257;
        %end
      end
    end
  end
end

rind = rind+1;
Aing(rind,:) = A(:)';
Bing(rind,1) = C1;
rind = rind+1;
Aing(rind,:) = -A(:)';
Bing(rind,1) = C1;


Ain = assembleblock({Aing});
Bin = [Bing];

g0 = randn(m)*1e-3;

[x,fval] = fmincon(@myobjm,g0(:),Ain,Bin,Aeq,Beq);
fval = sqrt(abs(fval));

gam = reshape(x,[m m]);

return



function y = myobjm(x)
y = -gfeat(x)^2;
return