function [A,fval] = sancheck460(m,h,d)
% searching over A -- NOT 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 = lenseq-1;
n = n+lx;

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

xdim = m^2*T;
adim = m^2*T;
fdim = m^n;

% box constraints
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

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


Ain = Aina;
Bin = Bina;

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

x0 = A0(:);

[x,fval] = fmincon(@myobjm1,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));
A = reshape(x(1:m^2*T),[m m T]);
%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});

global A C
A = reshape(x(1:m^2*T),[m m T]);
%mn = repmat([m],[1 n]);
%here, p0(2)=1

s = 0;
for x1=2
  for x2=1:m
    for x3=1:m
    %for x3=[1 3]
    %if ismember(x3,[1 2])
    %  ff = 1;
    %else
    %  ff = .2;
    %end
      %s = s + A(x3,x2)*( A(x2,1,1)*F([1 x2 x3 X]) - A(x2,x1,1)*F([x1 x2 x3 Y]) );
      s = s + A(x3,x2,2)*( A(x2,1,1)*F([1 x2 x3]) - A(x2,x1,1)*F([x1 x2 x3]) );
    end
  end
end

y = -s^2;
return



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

global A C
A = reshape(x(1:m^2*T),[m m T]);
C = zeros(m,1); C(2)=1;
%mn = repmat([m],[1 n]);

%here, p0(2)=1
gam = gamma(3,[]);
s=0;
for k=1:m
  gam = gamma(3,XX{k});
  s = s+ gam(k);
end

y = -s^2;
return


function gam = gamma(n,x)
global nstate
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



function r = F(x)
r = length(find(x~=1));
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