function [A,gam,fval,good] = gamcheck(m,h,c1,c2)
% conjecture:
% Suppose $\g_k[x]$ is such that
%\beq
%\sum_k \g_k[x^{(k)}] &\leq&
%C_1 + C_2\diam\{x^{(k)}\}.
%\eeq
%Then
%\beq
%\sum_k \sum_{k'} a_{kk'} \g_{k'}[x^{(k)}] &\leq&
%C_1 + C_2\h(1+\diam\{x^{(k)}\}).
%\eeq


% searching over A and gam
global XX XXX

global nstate lenseq
nstate = m;
n = 2;
lenseq = n;

if 0
  ncopies = 5;
  XXX={};
  for c=1:ncopies
    XX = cell(m,1);
    for i=1:m
      XX{i} = ceil(rand(1,n)*m);
    end
  %XX{:}
  XXX{c}=XX;
  end
end
T = 1;


% constraints on gamma
Aing = zeros(0,m^(n+1));
xinds = 1:m^n;  % all possible values of x
xdims = repmat([m],[1 n]);
mtupx = list_all_ind(repmat([m^n],[1 m]));

rind = 0;
for mtupind = 1:size(mtupx,1)
  xx = mtupx(mtupind,:);
  diamx = indiam(xx);
  for s = [-1 1];
    rind = rind+1;
    for k=1:m
      x = xx(k);
      cind = coord2ind([k x],[m m^n]);
      Aing(rind,cind) = s;
      Bing(rind,1) = c1+c2*diamx;
    end
  end
end

%gdim = size(Aing,2);
gdim = m^n*m;
adim = m^2;
xdim = gdim + adim;

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

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

% 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 i1 = 2
%  for i2 = [1 3]
for i1 = 1:m
  for i2 = i1+1:m
    for t=1:T
      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;
%Bina = Bina + 1;

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


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

g0 = randn(gdim,1)*1e-3;

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

[x,fval] = fmincon(@myobjm,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));
gam = x(1:gdim);
A = reshape(x(gdim+1:end),[m m T]);
%h = strict(A);
%fval
good = 1;
if 1
  for c=1:length(XXX)
    XX = XXX{c};
    y = myobjmX(x,XX);
    fval = sqrt(abs(y))
    diamX = diam(XX);
    Bound = c1 + c2*h*diamX
    c
    good = good & (fval <= Bound + 1e-8)
  end
end

if 0
  mtupx = list_all_ind(repmat([m^n],[1 m]));
  for mtupind = 1:size(mtupx,1)
    xx = mtupx(mtupind,:);
    y = y + myobjmxx(x,xx);
    fval = sqrt(abs(y))
    diamX = indiam(xx);
    Bound = c1 + c2*h*diamX
    good = good & (fval <= Bound + 1e-8);
  end
end
return


function y = myobjm(x)
global nstate lenseq
n = lenseq;
m = nstate;
global XXX
y=0;
%mtupx = list_all_ind(repmat([m^n],[1 m]));
%for mtupind = 1:size(mtupx,1)
%  xx = mtupx(mtupind,:);
%  y = y + myobjmxx(x,xx);
%end

for c=1:length(XXX)
  y = y + myobjmX(x,XXX{c});
end
return

function y = myobjmxx(x,xx)
global nstate lenseq
m = nstate;
n = lenseq;
T = 1;
gdim = m^n*m;

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

s = 0;
for k=1:m
  for k1=1:m
    %s = s + A(k,k1) * geval(gam,k1,XX{k});
    gind = coord2ind([k xx(k)],[m m^n]);
    s = s + A(k,k1) * gam(gind);
  end
end

y = -s^2;
return


function y = myobjmX(x,XX)
global nstate lenseq
m = nstate;
n = lenseq;
T = 1;
gdim = m^n*m;
%global XXX
%XX = XXX{c};


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

s = 0;
for k=1:m
  for k1=1:m
    s = s + A(k,k1) * geval(gam,k1,XX{k});
  end
end

y = -s^2;
return


function g = geval(gam,k,x)
global nstate lenseq
m = nstate;
n = lenseq;
T = 1;
gdim = m^n*m;

xdims = repmat([m],[1 n]);
gind = coord2ind([k x],[m xdims]);
g = gam(gind);
return


function r = F(x)
%r = length(find(x~=1));
%return
global nstate
lenseq = length(x);
global FN
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
r = FF(ind);
return


function p = p_0(x)
global C
p = C(x);
return

function p = p_n(i,j,n)
% primitive cond prob
global A
if i
  p = A(i,j,n);
else
  p = 1;
end
return


function d = indiam(xx)
global nstate lenseq
m = nstate;
n = lenseq;

d=0;
xdims = repmat([m^n],[1 n]);
for i=1:length(xx)
  for j=i+1:length(xx)
    xi = ind2coord(xx(i),xdims);
    xj = ind2coord(xx(j),xdims);
    d = max(d,rho(xi,xj));
  end
end

return