function [A,phi,fval] = sancheck4(h,d)
% searching over A and phi
% this "numerically proves" the theorem
% for n=1
global nstate 
m = nstate;

% the first 2m components of x are for A;
% the remaining m^2 components are for phi

xdim = 2*m+m^2;

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

% try killing box constraints
%LB = [];
%UB = [];
% nope, really do need them!!

% normalization constraints
Aeq = zeros(2,xdim);
for col=1:2
  for row=1:m
    Aeq(col,m*(col-1)+row) = 1;
  end
end
Beq = ones(2,1);

SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina = zeros(0,2*m);
rind = 0;
for col1 = 1
  for col2 = 2
    for s = 1:size(SIGNS,1)
      ss = SIGNS(s,:);
      rind = rind+1;
      for j=1:m
        ind1 = coord2ind([j col1],[m m]);
        ind2 = coord2ind([j col2],[m m]);
        Aina(rind,ind1) =  ss(j);
        Aina(rind,ind2) = -ss(j);
      end
    end
  end
end
Bina = ones(size(Aina,1),1)*h*2;

Ainf = zeros(0,m^2);
Binf = zeros(0,1);
rind = 0;
for ind1 = 1:m^2
  for ind2 = ind1+1:m^2
    x = ind2coord(ind1,[m m]);
    y = ind2coord(ind2,[m m]);
    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

d1a = size(Aina,1);
d1f = size(Ainf,1);

Ain = zeros(d1a+d1f,xdim);
Ain(1:d1a,1:2*m) = Aina;
Ain(d1a+1:end,2*m+1:end) = Ainf;

Bin = [Bina;Binf];

A0 = pnormdim(ones(m,2),1);
A0 = A0+randn(m,2)*1e-3;
A0 = pnormdim(A0,1);
FN = randFN(m,2,1,'not');
phi0 = FN{2};
x0 = [A0(:); phi0(:)];

%[x,fval] = fmincon(@myobjm,x0,Ain,Bin,Aeq,Beq,LB,UB);
[x,fval] = fmincon(@myobjm1,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));
A = reshape(x(1:2*m),[m 2]);
phi = reshape(x(2*m+1:end),[m^2 1]);
%keyboard
return

function y = myobjm(x)
global FN
global nstate 
m = nstate;
a = reshape(x(1:2*m),[m 2]);
FN{2} = reshape(x(2*m+1:end),[m^2 1]);
S = 0;
for j=1:m
  S = S + a(j,1)*F([1 j]) - a(j,2)*F([2 j]);
end
y = -abs(S)^2;
return


function y = myobjm1(x)
global FN
global nstate 
m = nstate;
a = reshape(x(1:2*m),[m 2]);
FN{2} = reshape(x(2*m+1:end),[m^2 1]);
S = 0;
for j=1:m
  S = S + (a(j,1) - a(j,2))*(F([1 j])+F([2 j]));
end
y = -abs(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