function hmmcheckmask
global A B C n mhid mobs
global mask

mhid = 3;
mobs = 3;
n = 3;

[Ainf,Binf] = lipab(mobs,n);

cnt = 0;

while 1
  cnt = cnt+1;

  A = pnormdim(rand(mhid,mhid),1);
  %A = pnormdim(A+100,1);
  C = pnormdim(rand(mhid,1   ),1);
  %C = pnormdim(ones(mhid,1   ),1);
  %C = C*0; C(2) = 1;
  B = pnormdim(rand(mobs,mhid),1);
  %B = eye(mhid);
  
  K = hmmphicoef(A,B,C,n,[1]);
  [x,fval] = linprog(K,Ainf,Binf,[],[],0*K,0*K+n);
  if fval < 0
    K = -K;
  end
  %mask = optmask(K,mobs,n); % nope, this does not yield the best mask
  mask = getmask1(x,mobs,n); % but this does seem to... or not...
  phim = fillgmask(mask,mobs);
  %if phim*K < 0
  %    'bla'
  %end
  
  good = abs(fval) < abs(phim*K) + 1e-9;
  if ~good
    'nope'
    keyboard
  end
  fprintf('cnt = %d, so far so good \n',cnt);
end


return


function mask = optmask(K,m,n)
%mask = {};
mask = repmat({[]},[1 n]);
for i=1:n
  mu = getmu(K,m,n);
  %phim = fillgmask(mu,m);
  mask{i} = mu;
  K = getK1(K,m,n);
  n = n-1;
end
return


function mu = getmu(K,m,n)
maxmu = [];
maxv = 0;
mask0 = repmat({[]},[1 n]);

for i=1:2^m
  setstr = dec2bin(i-1,m);
  ii1 = find(setstr=='1');
  mask = mask0;
  mask{1} = ii1;
  phim = fillgmask(mask,m);
  fv = phim*K;
  if fv > maxv
    maxv = fv;
    maxmu = ii1;
  end
end
mu = maxmu;
return