function [phi,fval] = afixed(A0,d)

global A
A = A0;

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

[m,m,T] = size(A);


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

% searching over phi

n = T;
nn = T+1;

global nstate lenseq
nstate = m;
lenseq = n;

adim = 0;
fdim = m^nn;
xdim = adim + fdim;


% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
UB = LB;
UB = UB + nn + d/3;  % no box constraints on phi
%UB = [];

%Aeq = zeros(1,fdim);
%Aeq(1) = 1;
%Beq = 0;
Aeq = []; Beq = [];

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

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

[x,fval] = fmincon(@myobjm,x0,Ainf,Binf,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));

y = myobjm(x);

global FN

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

return


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

adim = 0;
fdim = m^nn;
xdim = adim + fdim;

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'

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

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