function [A,gam,fval] = gamrm(m,h)
% exploring |Ag|
global nstate
nstate=m;

C = 1+h;

T = 1;

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


% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
LB(adim+1:end) = LB(adim+1:end) - Inf;  % no box constraints on gam
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + Inf;  % no box constraints on gam


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

% a quick'n'dirty hack to fix A
global A0
for i=1:m
  for j=1:m
    rind = rind+1;
    cind = coord2ind([i j],[m m]);
    Aeq(rind,cind) = 1;
    Beq(rind,1) = A0(i,j);
  end
end


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

% stricture constraints
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*2;

Aing = zeros(0,gdim);
Bing = zeros(0,1);
rind = 0;

% 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;
      end
    end
  end
end

% gam_ij sum constraint
rind = rind+1;
Aing(rind,:) = ones(1,gdim);
Bing(rind,1) = C*h;

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

A0 = pnormdim(ones(m,m,T),1);
A0 = A0+randn(size(A0))*1e-2;
A0 = pnormdim(A0,1);

gam0 = randn(m,m)*.1/C;
x0 = [A0(:); gam0(:)];

[x,fval] = fmincon(@myobjm,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));
A = reshape(x(1:adim),[m m T]);
gam = reshape(x(adim+1:end),[m m]);

return


function y = myobjm(x)
global nstate
m = nstate;
T = 1;
adim = m^2;
gdim = m^2;
xdim = adim + gdim;


A = reshape(x(1:adim),[m m T]);
gam = reshape(x(adim+1:end),[m m]);

%global A0 GG0
%A = A0;
%gam(:,1:m-1) = GG0(:,1:m-1);

s = sum(sum(A.*gam));
  
y = -s^2;
return