function [MU1,phi1,K1,maxs,sphi1] = boundsumm(K,m,n,MU,sum0)

maxs = 0;
dims = repmat([m],[1 n]);

t = length(MU)+1;

K1 = getK1(K,m,n);
[Ainf,Binf] = lipab(m,n-1);
[xp,fvalp] = linprog( K1,Ainf,Binf,[],[],[],[]);
[xn,fvaln] = linprog(-K1,Ainf,Binf,[],[],[],[]);
  
for mind = 0:2^m-1
  mset = dec2bin(mind,m);
  mu = find(mset=='1');
  
  MU1 = MU;
  MU1{end+1} = mu;
  
  sMU1 = sumkm(K,m,n,MU1,1);
  
  sp = abs(sum0 + sMU1 + fvalp);
  sn = abs(sum0 + sMU1 + fvaln);
  if sp>sn
    phi1 = xp(:)';
    s = sp;
    sphi1 = abs(fvalp);
  else
    phi1 = xn(:)';
    s = sn;
    sphi1 = abs(fvaln);
  end
  
  %s = evalsum(K,m,n,MU1,phi);
  if abs(s) > maxs
    maxs = s;
    maxMU1 = MU1;
    maxphi = phi1;
    maxsphi1 = sphi1;
    global mmind
    mmind = mind; %debugging
  end
end
MU1 = maxMU1;
phi1 = maxphi;
sphi1 = maxsphi1;
return




function s = evalsum(K,m,n,MU,phi)
s = 0;
dims = repmat([m],[1 n]);
t = length(MU);
dims1 = repmat([m],[1 n-t]);

for xind = 1:m^n
  x = ind2coord(xind,dims);
  c = K(xind);

  muval = applymask(MU,x(1:t));
  
  %xind1 = coord2ind(x(t+1):end,dims1);
  %phval = phi(xind1);
  %sx = muval + phval;
  
  sx = phi(xind) - muval;
  
  s = s + c*sx;
end

return


function K1 = getK1(K,m,n)
dims = repmat([m],[1 n]);
dims1 = repmat([m],[1 n-1]);

K1 = zeros(m^(n-1),1);

for xind2n = 1:m^(n-1)
  x2n = ind2coord(xind2n,dims1);
  for x1 = 1:m
    xind = coord2ind([x1 x2n],dims);
    K1(xind2n) = K1(xind2n) + K(xind);
  end
end
return


function r = applymask(mask,x)
n = length(x);
r = 0;
for i=1:n
  r = r + ismember(x(i),mask{i});
end
return


function y=mu(n,k)
global locm
mm = locm{n};
y = ismember(k,mm);
return