function [A,fval] = maskA(mask,m,n,h)
global locm
locm = mask;

global A

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

% searching over A


global nstate lenseq
nstate = m;
lenseq = n;

T = n-1;

adim = m^2*T;
fdim = 0;
xdim = adim + fdim;


% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
LB(adim+1:end) = LB(adim+1:end) - 1e6;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + 1e6;  % 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;

%ii1 = 1;
%ii2 = setdiff([1:m],ii1);

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

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

x0 = [A0(:);];

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

y = myobjm(x);


A = reshape(x(1:adim),[m m T]);

return


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

adim = m^2*T;
fdim = 0;
xdim = adim + fdim;

global A
%A = reshape(x,[m m T]);
A = reshape(x(1:adim),[m m T]);

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

% GG(k,k')
% summing across k'

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

y = -s^2;
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 y=mu(n,k)
global locm
mm = locm{n};
y = ismember(k,mm);
%y = (k==2);
return

function [Ak,Bk] = getAB(n)
global A
[m,m,t] = size(A);

Ak = zeros(m,1);
Bk = zeros(m,1);

if n==2
  for k=1:m
    Ak(k) = p_n(k,1,n-1)*mu(1,1) - p_n(k,2,n-1)*mu(1,2) + mu(2,k)*( p_n(k,1,n-1) - p_n(k,2,n-1) );
    Bk(k) = p_n(k,1,n-1) - p_n(k,2,n-1);
  end
else
  [Ak0,Bk0] = getAB(n-1);
  pkk = A(:,:,n-1);
  mk = zeros(m,1);
  for k=1:m
    mk(k) = mu(n,k);
  end
  
  Bk = pkk * Bk0;
  Ak = pkk * Ak0 + Bk.*mk;  
  
end
return


function gam = gamma(n,x)
global locm
[Ak,Bk] = getAB(n);
gam = Ak + Bk*applymask(locm(n+1:end),x);
return





function r = F(x)
global locm
r = applymask(locm,x);
return

function r = applymask(mask,x)
n = length(x);
r = 0;
for i=1:n
  r = r + ismember(x(i),mask{i});
end
return