function [v,x] = probdev_qpm(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
% actually, v is the max. variance
% the phi is constrained to be a mask

[m,m,t] = size(A);
if ~exist('C','var')
  C = ones(m,1)/m;
end

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

global P Lmf

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;

V = diag(P) - P*P';
% F' * V * F = Var[F]

if size(Lmf,2) ~= size(V,1)
  Lmf = getLmf(m,n);
end

V = Lmf' * V * Lmf;
% since phi = Lmf * mu

global myopt
[x,fval] = quadprog(-V,zeros(size(LB)),[],[],[],[],LB,UB,[],myopt);
v = -fval*2;
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