function [A,phi,fval] = sancheck47(m,h,d)
% searching over A and phi
% this "numerically proves" the theorem
global X Y XX

global nstate lenseq
nstate = m;
n = 3;
lenseq = n;

lx = 0;
%X = ceil(rand(1,lx)*nstate)
%Y = ceil(rand(1,lx)*nstate)
XX = cell(m,1);
for i=1:m
  XX{i} = ceil(rand(1,lx)*nstate);
end


T=1;

% normalization constraints for A1
Aeq1 = zeros(m,m^2);
for i=1:m
  for j=1:m
    Aeq1(i,m*(i-1)+j) = 1;
  end
end
Beq1 = ones(size(Aeq1),1);

% normalization constraints for A2
Aeq2 = zeros(m,2*m);
for i=1:2
  for j=1:m
    Aeq2(i,m*(i-1)+j) = 1;
  end
end
Beq2 = ones(size(Aeq2),1);


% stricture constraints
SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina1 = zeros(size(SIGNS,1)*m*(m-1)/2*T,m^2*T);
rind = 0;
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]);
          Aina1(rind,ind1) =  ss(j);
          Aina1(rind,ind2) = -ss(j);
        end
      end
    end
  end
end
Bina1 = ones(size(Aina1,1),1)*h*2;


Aina2 = 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]);
        Aina2(rind,ind1) =  ss(j);
        Aina2(rind,ind2) = -ss(j);
      end
    end
  end
end
Bina2 = ones(size(Aina2,1),1)*h*2;


Ainf = zeros(0,m^n);
Binf = zeros(0,1);
mn = repmat([m],[1 n]);
rind = 0;
for ind1 = 1:m^n
  for ind2 = ind1+1:m^n
    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({Aina1,Aina2,Ainf});
%Aeq = assembleblock({Aeq1,Aeq2});
Ain = assembleblock({Aina2,Aina1,Ainf});
Aeq = assembleblock({Aeq2,Aeq1,zeros(0,m^n)});

Bin = [Bina2;Bina1;Binf];
Beq = [Beq2; Beq1];


A1 = pnormdim(ones(m,m),1);
A1 = A1+randn(size(A1))*1e-3;
A1 = pnormdim(A1,1);
A2 = pnormdim(ones(m,2),1);
A2 = A2+randn(size(A2))*1e-3;
A2 = pnormdim(A2,1);

FN = randFN(m,n,d,'not');
phi0 = FN{n};
%x0 = [A1(:); A2(:); phi0(:)];
x0 = [A2(:); A1(:); phi0(:)];


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


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


function y = myobjm(x)
global FN
global nstate lenseq
m = nstate;
n = lenseq;
T = n-1;
global X Y XX
n = n+length(XX{1});

A2 = reshape(x(1:m^2),[m m]);
A1 = reshape(x(m^2+1:m^2+2*m),[m 2]);
FN{n} = reshape(x(m^2+2*m+1:end),[1 m^n]);


%here, p0(2)=1



s = 0;
for x1=2
  %for x2=1:m
  for x2=[1 2]
    for x3=1:m
    %for x3=[1 3]
      s = s + A2(x3,x2)*( A1(x2,1)*F([1 x2 x3]) - A1(x2,x1)*F([x1 x2 x3]) );
    end
  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