function [mval,y] = optmask0(A0)
global A mask
A = A0;
[m,m,T] = size(A);
n = T+1;

global nstate
nstate = m;

y = zeros(1,n);
dims = repmat([m],[1 n]);

FV = [];

for yind = 1:m^n
  mask = ind2coord(yind,dims);
  FV(end+1) = pheval;
end
[mval,yind] = max(FV);
y = ind2coord(yind,dims);
return



function r = F(x)
global mask
r = rho(x,mask);
return
global nstate
lenseq = length(x);
global FN
FF = FN{lenseq};
ind = coord2ind(x,(repmat(nstate,1,lenseq)));
r = FF(ind);
return



function b = myeq(x,y)
b = abs(x-y) < 1e-9;
return

function r = rho(x,y)
r = sum(x~=y);
return



function y = pheval
global A
[m,m,T] = size(A);
n = T;

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

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

y = abs(s);
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