function [A,gam,fval] = agcheck(m,n,h,d)
% exploring |Ag|
global nstate lenseq loch locd
nstate=m;
lenseq=n;
loch = h;
locd = d;

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

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

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


% 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);
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==1) | (j==2))
      if ((j==3))
      %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



% 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);
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:
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) = .06;
        %Bing(rind,1) = d*h^(n-1);
      end
    end
  end
end




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);
g0 = rand(m,m)/1000;


x0 = [A0(:); g0(:)];

[x,fval] = fmincon(@myobjm,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));


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

%C = d*sum(h.^(0:n));
%good = (fval <= C + 1e-10)
return


function y = myobjm(x)
global nstate lenseq
m = nstate;

adim = m^2;
gdim = m^2;

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

s = sum(sum(A.*gam));

y = -s^2;
return


function y = myobjm1(x)
global nstate lenseq loch locd
m=nstate;
n=lenseq;
h=loch;
d=locd;

C = d*sum(h.^(0:n));
C = C+.01;

adim = m^2;
gdim = m^2;

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

s = sum(sum(A.*gam));

%y = (s^2-C^2)^2;
y = C-s^2;

%y = -s^2;
return