function [A,fval,diamX] = gamcheck1(m,n,h)
% 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

global C
C = zeros(m,1); C(2)=1;

% searching over A and gam
global XX


global nstate lenseq
nstate = m;
lenseq = n;

lx = 1;
XX = cell(m,1);
for i=1:m
  XX{i} = ceil(rand(1,lx)*m);
end
%XX{:}
diamX = diam(XX);

global FN
FN = randFN(m,n+lx,1,'inexact');

T = n-1;

adim = m^2*T;

% box constraints
LB = zeros(adim,1)+1e-10;
UB = ones(adim,1);


% normalization constraints
Aeq = zeros(m*T,adim);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,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 t=1:T
  %if (T==1)
  %  ii1 = 1;
  %  ii2 = setdiff([1:m],ii1);
  %else
    ii1 = 2;
    ii2 = setdiff([1:m],ii1);
  %end
  for i1 = ii1
    for i2 = ii2
      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;

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

x0 = A0(:);

[x,fval] = fmincon(@myobjm,x0,Aina,Bina,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));
A = reshape(x,[m m T]);
%h = strict(A);
%fval
%Bound = c1 + c2*h*(1+diam(XX))
%good = (fval <= c1 + c2*h*(diam(XX)) + 1e-8)
return


function y = myobjm(x)
global nstate lenseq
m = nstate;
n = lenseq;
T = n-1;
global XX
global A
A = reshape(x,[m m T]);


s = 0;
for k=1:m
  gamk = gamma(n-1,XX{k});
  for k1=1:m
    s = s + A(k,k1,T) * gamk(k1);
  end
end
  
y = -s^2;
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 gam = gamma(n,x)
%global nstate
global A
[m,m,t] = size(A);
nstate = m;
gam = zeros(nstate,1);
if (n==2)
  for k=1:nstate
    for x1=1:nstate
      gam(k) = gam(k) + p_0(x1)*(p_n(k ,1 ,1)*F([1 k x]) - p_n(k ,x1,1)*F([x1 k x]));
    end
  end
else
  for k=1:nstate
    gamkx = gamma(n-1,[k x]);
    for k1=1:nstate
      gam(k) = gam(k) + p_n(k,k1,n-1) * gamkx(k1);
    end
  end
end  
return