function [A,gam,fval] = gam3m(c1,c2,h,m)
global C1 C2 H
C1 = c1;
C2 = c2;
H = h;

global nstate lenseq
nstate = m;
%lenseq = n;

T = 1;

adim = m^2*T;
gdim = 2*m*m;
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) - 1e6;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + 1e6;  % 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);


[Aina,Bina] = constrict(m,h,T);

rind=0;
Aing = zeros(0,gdim);
mm = repmat([2],[1 m]);
minds = list_all_ind(mm);
for k=1:m
  for s=[-1 1]
    for jind = 1:size(minds,1)
      rind = rind+1;
      iseq = minds(jind,:);
      rho = ~all(iseq==iseq(1));
      for j=1:m
        i = iseq(j);
        cind = coord2ind([i j k],[2 m m]);
        Aing(rind,cind) = s;
        Bing(rind,1) = c1+c2*rho;
      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 = randn(2,m,m)/100;
%G0(:,1) = G0(:,1) + 1e6;
%G0(:,2) = G0(:,2) - 1e6;

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 T]);
gam = reshape(x(adim+1:end),[2 m m]);
%g = gam(2,:);
%gg = repmat(g,[m-2,1]);
%gam = [gam; gg];

return


function y = myobjm(x)
global nstate lenseq
global C1 C2 H

global Gkl

m = nstate;
%n = lenseq;
T = 1;

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

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

Gkl = zeros(m);

for k = 1:m
  gamk = gam(:,:,k);
  g = gamk(2,:);
  gg = repmat(g,[m-2 1]);
  gamk = [gamk; gg];
  
  for l = 1:m
    a = A(l,:);
    g = gamk(l,:);
    Gkl(k,l) = g*a';
  end
end

h = H;


s = sum(diag(Gkl));
%s = sum(Gkl(1,:));
k1 = C1 + C2*h;
k2 = C2*h;
%
y = k1 + k2 - abs(s);


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

%g = gam(2,:);
%gg = repmat(g,[m-2,1]);
%gam = [gam; gg];

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

%y = C1 + C2*h - abs(s);


return


function [Aina,Bina] = constrict(m,h,T)
% stricture constraints

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

SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina = zeros(0,m^2*T);
rind = 0;

%for i1=1
for i1=1:m-1
%for i2=2
  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],[m m]);
        ind2 = coord2ind([j i2],[m m]);
        Aina(rind,ind1) =  ss(j);
        Aina(rind,ind2) = -ss(j);
      end
    end
  end
end
Bina = ones(size(Aina,1),1)*h;
Bina = Bina * 2;
%keyboard
return