function philp = optmasklp(A0)
global A mask
A = A0;
[m,m,T] = size(A);
n = T+1;
global FN
FN{n} = rand(1,m^n);


global nstate lenseq
nstate = m;
lenseq = n;

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

global ALLSEQ ALLSEQ1
i0=1;
Xi = 1;
ALLSEQ =  list_all_ind(repmat(nstate,1,lenseq));
ALLSEQ1 = list_all_ind(repmat(nstate,1,lenseq-i0+1));

d = 1;
cc = phicoeff(Xi);

global Ainf Binf
if size(Ainf,2)~=m^n
  Ainf = [];
end
if isempty(Ainf)
  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;
        %if rind==131
        %  global x0 y0
        %  ind1,ind2,x,y
        %  x0=x,y0=y
        %  return
        %end
        Ainf(rind,ind1) =  s;
        Ainf(rind,ind2) = -s;
        Binf(rind,1) = d;
      end
      end
    end
  end
end
cc = cc(:);

%LB = zeros(size(cc))-1e-9;
%UB = LB + d*n + 1e-9;

[xp,fvalp] = linprog( cc,Ainf,Binf,[],[],[],[]);
[xn,fvaln] = linprog(-cc,Ainf,Binf,[],[],[],[]);
if abs(fvalp)>abs(fvaln)
  x = xp;
  fval = fvalp;
else
  x = xn;
  fval = fvaln;
end

philp = x(:)';
fvalp = abs(fval);

%[mval,y] = optmask(A0);
%if abs(fvalp-mval) > 1e-9
%  'not good'
%  keyboard
%end

return



function r = F(x)
%r = sum(x~=1);
%return
%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 cc = phicoeff(Xi)
% Y1 = E[f(X)|X_i] - E[f(X)|X_{i-1}]
global ALLSEQ1 A nstate lenseq
m = nstate;
n = lenseq;

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

i = length(Xi);
%n = size(ALLSEQ1,2);
cc = zeros(m^n,1);

for j=1:size(ALLSEQ1,1)
  x = ALLSEQ1(j,:);
  x = [x 0];
  Xix = [Xi(1:i) x(2:end)];
  Xi_1x = [Xi(1:i-1) x];
  a = pcp(Xix(i+1),Xix(i),i);
  b = pcp(Xi_1x(i+1),Xi_1x(i),i);
  p = PCOND(x,Xi(1:i-1));
  p = p / (b+realmin);

  % f0
  xind = coord2ind(Xi_1x(1:end-1),dims);
  cc(xind) = cc(xind) - p*b;

  % f1
  xind = coord2ind(Xix(1:end-1),dims);
  cc(xind) = cc(xind) + p*a;
  
  %f0 = F(Xi_1x);
  %f1 = F(Xix);
  %y = y + p*(a*f1-b*f0);
end
return


function Ef = E(Xi)
% E[f(X)|Xi]
global ALLSEQ A
Ef = 0;
i = length(Xi);
n = size(ALLSEQ,2);
for k=1:size(ALLSEQ,1)
  x = ALLSEQ(k,:);
  p = P(x);
  if (i & (i<n))
    pxi = A(x(i+1),x(i),i);
    pXi = A(x(i+1),Xi(i),i);
    p = p / pxi * pXi;
  end
  f = F([Xi x(i+1:end)]);
  Ef = Ef + f*p;
end
return


function y = Y(Xi)
% Y1 = E[f(X)|X_i] - E[f(X)|X_{i-1}]
Ef1 = E(Xi);
Ef0 = E(Xi(1:end-1));
y = Ef1 - Ef0;
return



function p = P(x)
global A C
if isempty(x)
  p = 1;
else
  p = C(x(1));
end
for j=2:length(x)
  p = p * A(x(j),x(j-1),j-1); % markov case
end
return

function p = PCOND(x,Xi,t)
global C
if isempty(Xi)
  p = C(x(1));
else
  p = pcp(x(1),Xi(end),t);
end
for j=2:length(x)
  p = p * pcp(x(j),x(j-1),j-1); % markov case
end
return

function p = pcp(i,j,t)
% primitive cond prob
global A
if i
  p = A(i,j,t);
else
  p = 1;
end
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