function [maxp,maxx] = maxprobdev(A0,C0)
% 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(A0);
if ~exist('C0','var')
  C0 = ones(m,1)/m;
end
global A C
A = A0; C = C0;
n = t+1;
global P Ainf Binf
LB = zeros(m^n,1);
UB = LB + n;
%phi = rand(size(LB));
% initialize w/quadprog
global domask phi01
if domask
  [v,phi] = probdev_qpm(A,C);
else
  [v,phi] = probdev_qp(A,C);
end

global myopt
maxp = 0;
maxx = [];
global niter
for i=1:niter
  if domask
    [x,fval] = fmincon(@myobjm,phi,[],[],[],[],zeros(size(phi)),ones(size(phi)),[],myopt);
  elseif phi01
    [x,fval] = fmincon(@myobj,phi,[],[],[],[],LB,LB+1,[],myopt);
  else
    [x,fval] = fmincon(@myobj,phi,Ainf,Binf,[],[],LB,UB,[],myopt);
  end
  
  p = -fval;
  if p > maxp
    maxp = p;
    maxx = x;
  end
end
return

function y = myobj(phi)
global A C
p = probdev(A,C,phi);
y = -p;
return

function y = myobjm(mu)
global A C Lmf 
phi = Lmf*mu;
p = probdev(A,C,phi);
y = -p;
return