function [v,phi,P] = probdev_qp(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

[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 Ainf Binf

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

%[Ainf,Binf] = lipab(m,n);
LB = 0*P;
UB = LB + n;

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

global myopt phi01
if phi01
  [x,fval] = quadprog(-V,zeros(size(P)),[],[],[],[],LB,LB+1,[],myopt);
else
  [x,fval] = quadprog(-V,zeros(size(P)),Ainf,Binf,[],[],LB,UB,[],myopt);
end
%[x,fval] = quadprog(-V,zeros(size(P)),[],[],[],[],LB,UB);
v = -fval*2;
phi = round(x(:)');
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