function [A,phi,fval] = pgamcheck(m,n,h,d)
% partial bound check
% for now, m=2

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

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

% searching over A and gam

%nn = 3;
nn = n+1;

global nstate lenseq
nstate = m;
lenseq = n;

T = nn-1;

adim = m^2*T;
fdim = m^nn;
xdim = adim + fdim;


% 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) - 1000;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + 1000;  % 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);


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

%ii1 = 1;
%ii2 = setdiff([1:m],ii1);

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;

Ainf = zeros(0,m^nn);
Binf = zeros(0,1);
mn = repmat([m],[1 nn]);
rind = 0;
for ind1 = 1:m^nn
  for ind2 = ind1+1:m^nn
    x = ind2coord(ind1,mn);
    y = ind2coord(ind2,mn);
    if length(find(x~=y))==1  % Hamming dist == 1
      for s = [-1,1]
        rind = rind+1;
        Ainf(rind,ind1) =  s;
        Ainf(rind,ind2) = -s;
        Binf(rind,1) = d;
      end
    end
  end
end

Ain = assembleblock({Aina,Ainf});
Bin = [Bina;Binf];

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

global FN
FN = randFN(m,nn,d,'not');
phi0 = FN{nn};
x0 = [A0(:); phi0(:)];

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

y = myobjm1(x);

global A FN

A = reshape(x(1:adim),[m m T]);
phi = reshape(x(adim+1:end),[1 m^nn]);
FN{end} = phi;

%gam1 = gamma(n,[1]);
%gam2 = gamma(n,[2]);
%fval = max(abs(gam1-gam2));
%fval = abs(max(gam1)-min(gam1));

return

function y = myobjm(x)
global nstate lenseq
global FN
nn = length(FN);
m = nstate;
n = lenseq;
T = nn-1;
global GG

adim = m^2*T;
fdim = m^nn;
xdim = adim + fdim;

global A
%A = reshape(x,[m m T]);
A = reshape(x(1:adim),[m m T]);
FN{nn} = reshape(x(adim+1:end),[1 m^nn]);

GG = [];
for k=1:m
  gam = gamma(n,[k]);
  GG = [GG; gam'];
end

% GG(k,k')
% summing across k'


aa = A(:,:,T);
%s = sum(sum(a.*GG));

%a0 = aa(1,1);
%a1 = aa(2,1);
%b0 = aa(1,2);
%b1 = aa(2,2);

%s1 = (a0+b0)*(GG(1,1) + GG(1,2)) + (a1+b1)*(GG(2,1) + GG(2,2));

%s2 = (a0-b0)*(GG(1,1) - GG(1,2)) + (a1-b1)*(GG(2,1) - GG(2,2));

%s = (abs(s1)+abs(s2))/2;


i0 = 1;
i1 = 2;

a = aa(:,i0);
b = aa(:,i1);

s2 = 0;
for k=1:m
  s2 = s2 + (a(k)-b(k))*(GG(k,i0)-GG(k,i1));
end

s1 = 0;
for k=1:m
  s1 = s1 + (a(k)+b(k))*(GG(k,i0)+GG(k,i1));
end

ii = 1:m; ii = setdiff(ii,[i0 i1]);
for j=ii
  for k=1:m
    s1 = s1 + 2*aa(k,j) * GG(k,j);
  end
end

s = (abs(s1)+abs(s2))/2;

y = -s^2;
return


function y = myobjm1(x)
global nstate lenseq
global FN
nn = length(FN);
m = nstate;
n = lenseq;
T = nn-1;
global GG

adim = m^2*T;
fdim = m^nn;
xdim = adim + fdim;

global A
%A = reshape(x,[m m T]);
A = reshape(x(1:adim),[m m T]);
FN{nn} = reshape(x(adim+1:end),[1 m^nn]);

GG = [];
for k=1:m
  gam = gamma(n,[k]);
  GG = [GG; gam'];
end

% GG(k,k')
% summing across k'


aa = A(:,:,T);
%s = sum(sum(a.*GG));

%a0 = aa(1,1);
%a1 = aa(2,1);
%b0 = aa(1,2);
%b1 = aa(2,2);

%s1 = (a0+b0)*(GG(1,1) + GG(1,2)) + (a1+b1)*(GG(2,1) + GG(2,2));

%s2 = (a0-b0)*(GG(1,1) - GG(1,2)) + (a1-b1)*(GG(2,1) - GG(2,2));

%s = (abs(s1)+abs(s2))/2;


i0 = 1;
i1 = 2;

a = aa(:,i0);
b = aa(:,i1);

s1 = 0;
for k=1:m
  s1 = s1 + (a(k)-b(k))*(GG(k,i0)-GG(k,i1));
end

s2 = 0;
for k=1:m
  %s2 = s2 + (a(k)+b(k))*(GG(k,i0)+GG(k,i1));
  for i=1:m
    s2 = s2 + (a(k)+b(k))*(GG(k,i));
  end
end

s3 = 0;
ii = 1:m; ii = setdiff(ii,[i0 i1]);
for j=ii
  for k=1:m
    s3 = s3 + (2*aa(k,j)-a(k)-b(k)) * GG(k,j);
  end
end

%s = (abs(s1)+abs(s2))/2;
%s = (abs(s1)+abs(s2)+abs(s3))/2;  % this is BAD
%s = (abs(s1+s3)+abs(s2))/2;  % this is BAD

global GS1 GS2
GS1 = (abs(s1))/2;
GS2 = (abs(s2+s3))/2;
%s = (abs(s1)+abs(s2+s3))/2;
s = GS1;

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