function trymaskbound
global m n Ainf Binf
global A B C n mhid mobs
mhid = 3;
mobs = 3;
n = 3;

m = mobs;

myopt = optimset('MaxIter',1000);
%myopt = optimset('MaxIter',10000000);

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

adim = mhid^2;
bdim = mhid*mobs;
fdim = mobs^n;
%xdim = adim+bdim+fdim;
xdim = adim+bdim;

% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
%LB(adim+1:end) = LB(adim+1:end) - 1e6;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+bdim+1:end) = UB(adim+bdim+1:end) + n;  % no box constraints on phi


% normalization constraints
Aeqa = zeros(mhid,adim);
for i=1:mhid
  for j=1:mhid
    Aeqa(i,mhid*(i-1)+j) = 1;
  end
end
Beqa = ones(mhid,1);

Aeqb = zeros(mhid,bdim);
rind = 0;
for c=1:mhid
  rind = rind+1;
  for r=1:mobs
    ind = coord2ind([r c],[mobs mhid]);
    Aeqb(rind,ind) = 1;
  end
end
Beqb = ones(mhid,1);

Aeq = assembleblock({Aeqa,Aeqb});
Beq = [Beqa;Beqb];

global maxK maxr
maxr = 0;


while 1
  %K = randn(m^n,1);
  A = pnormdim(rand(mhid,mhid),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);

  %global A0 B0 C0
  %A = A0, B = B0, C = C0

  K = hmmphicoef(A,B,C,n,[1]);
  
  while 0
    [xp,fvalp] = linprog(K,Ainf,Binf,[],[],0*K,0*K+n);
    [xn,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n);
    if abs(fvalp)>abs(fvaln)
      fv = fvalp;
      x = xp;
    else
      fv = fvaln;
      x = xn;
    end
    
    if x'*K < 0
      K = -K;
    end
    
    %mask = getmask1(x,m,n);
    mask = optmask(K,m,n);
    phim = fillgmask(mask,m);
  
    r = abs(fv)/abs(phim*K)
    if r > 1,
      %r
      break
    end
  end
  
  x = [A(:); B(:)];
  %[K,fval] = fmincon(@myobj,K,[],[],[],[],0*K-1,0*K+1);
  [x,fval] = fmincon(@myobj1,x,[],[],Aeq,Beq,LB,UB,[],myopt);
  
  r = sqrt(abs(fval));
  
  if r > maxr
    maxK = K;
  end
  
  maxr = max(r,maxr);
  fprintf('                   maxr = %1.5f \n',maxr);
end

return


function y = myobj(K)
global m n Ainf Binf

[xp,fvalp] = linprog( K,Ainf,Binf,[],[],0*K,0*K+n);
[xn,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n);
if abs(fvalp)>abs(fvaln)
  fv = fvalp;
  x = xp;
else
  fv = fvaln;
  x = xn;
end

if x'*K < 0
  K = -K;
end

mask = optmask(K,m,n);
phim = fillgmask(mask,m);

y = -(fv/(phim*K))^2;

  
r = sqrt(abs(y));

global maxK maxr

if r > maxr
  maxr = r;
  maxK = K;
end
  
fprintf('                           maxr = %1.5f \n',maxr);

return

function y = myobj1(x)
global A B C n mhid mobs phi
global m Ainf Binf
global mask

adim = mhid^2;
bdim = mhid*mobs;
fdim = mobs^n;
%xdim = adim+bdim+fdim;
xdim = adim+bdim;

A = reshape(x(     1:adim     ),[mhid mhid]);
B = reshape(x(adim+1:adim+bdim),[mobs mhid]);

K = hmmphicoef(A,B,C,n,[1]);

if 0
  [xp,fvalp] = linprog( K,Ainf,Binf,[],[],0*K,0*K+n);
  [xn,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n);
  if abs(fvalp)>abs(fvaln)
    fv = fvalp;
    x = xp;
  else
    fv = fvaln;
    x = xn;
  end

  if x'*K < 0
    K = -K;
  end
else
  psi1 = Psi( K,m,n);
  psi2 = Psi(-K,m,n);
  if psi1 > psi2
    fv = psi1;
  else
    fv = psi2;
    K = -K;
  end
end

mask = optmask(K,m,n);
phim = fillgmask(mask,m);

y = -(fv/(phim*K))^2;

  
r = sqrt(abs(y));

global maxK maxr maxA maxB maxC

if r > maxr
  maxr = r;
  maxK = K;
  maxA = A;
  maxB = B;
  maxC = C;
  
  save HMMASKDAT maxK maxr maxA maxB maxC
end
  
fprintf('                           maxr = %1.5f \n',maxr);


return



function K=genK(m,n)
K = randn(m^n,1);
return



function mask = optmask(K,m,n)
%mask = {};
mask = repmat({[]},[1 n]);
for i=1:n
  mu = getmu(K,m,n);
  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