function p = probdev_optmask(A,C)
% the probability that the rv. phi(X), with X ~ Markov(A,C) 
% deviates from its mean by at least 1
% finds "worst" such phi
global m n

[m,m,t] = size(A);
n = t+1;
dims = repmat([m],[1 n]);

global P

P = zeros(m^n,1);
F = zeros(m^n,1);
for xind=1:m^n
  x = ind2coord(xind,dims);
  P(xind) = Pr(A,C,x);
  %F(xind) = phi(xind);
end

LB = zeros(m*n,1);
UB = LB + 1;

pp = [];
while 1
  phim = 1+randn(m,n)/100;
  [x,fval] = fmincon(@myobj,phim(:),[],[],[],[],LB,UB);
  p = -fval;
  pp(end+1) = p;
  if length(pp)>40
    break
  end
end
pp
p = max(pp);

return

function y = myobj(x)
global P m n
F = mat2mask(x);
Ef = P'*F;
del = abs(F-Ef);

ii = find(del>=1);
p = sum(P(ii));

y = -p;
return

function F = mat2mask(mask);
global m n
mask = reshape(mask,[m n]);
F = zeros(m^n,1);
dims = repmat([m],[1 n]);
for xind=1:m^n
  x = ind2coord(xind,dims);
  for i=1:n
    F(xind) = F(xind) + mask(x(i),i);
  end
end
return


function p = Pr(A,C,x)
%p = 1;
p = C(x(1));
for j=2:length(x)
  p = p * A(x(j),x(j-1),j-1);
end
return