function [gam,fval] = gamrmlp(A,d,n)
% exploring |Ag|
global nstate
nstate=size(A,1);
m = nstate;
h = strict(A);
%n = 2;

if ~exist('d','var')
  d=1;
end

C = d*sum(h.^(0:n-1));


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

% box constraints
% to avoid affine degeneracy
%LB = -ones(gdim,1)+1e-10;
%LB = zeros(gdim,1)+1e-10;
%UB =  ones(gdim,1)-1e-10;
LB = [];
UB = [];

rind=0;
Aeq = zeros(0,xdim);
Beq = zeros(0,1);
if 1
  % a quick'n'dirty hack to fix gam
  global GG0
  for i=1:m
    for j=1:m
      if ((i==1) | (i==2))
        %if (i>=1) | (j>=1)
        %if (j>=3) & (i>=3)
        %if j>=2;
        %if ((j==1) | (j==3))
        %if ((j==1))
        %if (j>=2) | (i>=2)
        %if 0
        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:
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) = .1;
        %Bing(rind,1) = 0.05;
        %Bing(rind,1) = d*h^(n-1);
        %Bing(rind,1) = 1e-10;
        %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

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

f = A(:);


[xp,fvalp] = linprog( f,Ain,Bin,Aeq,Beq,LB,UB);
[xn,fvaln] = linprog(-f,Ain,Bin,Aeq,Beq,LB,UB);
if abs(fvalp)>abs(fvaln)
  x = xp;
  fval = fvalp;
else
  x = xn;
  fval = fvaln;
end

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

return