function [K,phi,fval] = trynbreak_selfsuff(K,phi,m0,n0,comstr)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt
%warning off
warning on

mylinopt = optimset('Diagnostics','off');
mylinopt = optimset('Display','off');

m = m0; n = n0;

if isempty(K)
  K = randn(m^n,1);
end

if isempty(phi)
  phi = randlip(m,n)';
end

n1 = n-1;
[Ainl,Binl] = lipab(m,n);
[Ain1,Bin1] = lipab(m,n1);

dum = phi(:);

LBK = -ones(size(dum));
UBK =  ones(size(dum));

LBf =   zeros(size(dum));
UBf =  n*ones(size(dum));

LB = [LBf; LBK];
UB = [UBf; UBK];

Ain = assembleblock({Ainl,zeros(0,m^n)});
Bin = [Binl];

x0 = [phi(:); K(:)];


switch comstr
 case 'debug'
  myobjm_debug;
 case 'greedy'
  [x,fval] = fmincon(@myobjm_greedy,x0,Ain,Bin,[],[],LB,UB);
 case 'psi'
  [x,fval] = fmincon(@myobjm_psi,x0,Ain,Bin,[],[],LB,UB);
 case 'phi1'
  [x,fval] = fmincon(@myobjm_phi1,x0,Ain,Bin,[],[],LB,UB);
 case 'phi1counterx'
  UBf =  n1*ones(size(dum));
  UB = [UBf; UBK];
  [x,fval] = fmincon(@myobjm_phi1counterx,x0,Ain,Bin,[],[],LB,UB);
 case 'phi1counterx2'
  [x,fval] = fmincon(@myobjm_phi1counterx2,x0,Ain,Bin,[],[],LB,UB);
 case 'pos0'
  %LBK = zeros(size(dum));
  %LB = [LBf; LBK];
  [x,fval] = fmincon(@myobjm_pos0,x0,Ain,Bin,[],[],LB,UB);
 case 'pos'
  UBf =  n1*ones(size(dum));
  UB = [UBf; UBK];
  LBK = zeros(size(dum));
  LB = [LBf; LBK];
  [x,fval] = fmincon(@myobjm_pos,x0,Ain,Bin,[],[],LB,UB);
 case 'posneg'
  [x,fval] = fmincon(@myobjm_posneg,x0,Ain,Bin,[],[],LB,UB);
 case 'zajim'
  [x,fval] = fmincon(@myobjm_zajim,x0,Ain,Bin,[],[],LB,UB);
 case 'conv'
  [x,fval] = fmincon(@myobjm_conv,x0,Ain,Bin,[],[],LB,UB);
 case 'phir'
  [x,fval] = fmincon(@myobjm_phir,x0,Ain,Bin,[],[],LB,UB);
 case 'phim'
  [x,fval] = fmincon(@myobjm_phim,x0,Ain,Bin,[],[],LB,UB);
 case 'lips'
  [x,fval] = fmincon(@myobjm_lips,x0,Ain,Bin,[],[],LB,UB);
 case 'greedym'
  [x,fval] = fmincon(@myobjm_greedym,x0,Ain,Bin,[],[],LB,UB);
 case 'phikshlp'
  [x,fval] = fmincon(@myobjm_phikshlp,x0,Ain,Bin,[],[],LB,UB);
 otherwise
  error('unknown obj function');
end
phi = x(1:m^n)';
K = x(m^n+1:end);
%if phi*K<0
%  K = -K;
%end
return


function y = myobjm_phikshlp(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
ksig = (K>0);
phiks=pl(phi-ksig');

K1 = getK1(K,m,n);

[x,fval] = linprog(-K1,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
phi1 = x';

Ksrt = sort(K);
Fsrt = sort(phiks);
fprintf('.');

y = phi1*K1 - Fsrt*Ksrt;
return



function y = myobjm_greedym(x)
global Ainl Binl Ain1 Bin1 GG iif m n

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
ksig = (K>0);
phiks=pl(phi-ksig');

K1 = getK1(K,m,n);

phig = greedymultispike(K1,m,n1);
fprintf('.');

y = phig*K1 - phiks*K;
return

function y = myobjm_debug
global Ainl Binl Ain1 Bin1 GG iif m0 n0 K0 phi0 mylinopt
K=K0,phi=phi0,m=m0,n=n0

n1 = n-1;
ksig = (K>0);
phiks=pl(phi-ksig');

K1 = getK1(K,m,n);

keyboard

phig = greedymultispike(K1,m,n1);
%fprintf('.');

y = phig*K1 - phiks*K;
return






function y = myobjm_phi1counterx(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
K1 = getK1(K,m,n);

[x,fval] = linprog(-K1,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
phi1 = x';

y = phi1*K1 - phi*K;
% DEF not true that phi1 bounds general phi:X^n --> {0,...,n-1}
return


function y = myobjm_phi1counterx2(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
K1 = getK1(K,m,n);

phiks = pl(phi-1);

[x,fval] = linprog(-K1,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
phi1 = x';

y = phi1*K1 - phiks*K;
% CONFIRMED YET AGAIN: need real phiks=phi-ksig, not phi-1
return


function y = myobjm_lips(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
K1 = getK1(K,m,n);
ksig = (K>0);
phiks=pl(phi-ksig');

Kf = getK1(phiks'.*K,m,n);

phil = fitlips(K1,Kf,m,n1);
fprintf('.');

if ~islip(phil)
  'not lip'
  keyboard
end

y = phil*K1 - phiks*K;
return



function y = myobjm_phim(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
K1 = getK1(K,m,n);
ksig = (K>0);
phiks=pl(phi-ksig');

phim = meanphi(phi,m,n);
phim = ceil(phim);
phim = min(phim,n1);

if ~islip(phim)
  'not lip'
  keyboard
end

y = phim*K1 - phiks*K;
return



function y = myobjm_phir(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

phi = x(1:m^n)';
K = x(m^n+1:end);

n1 = n-1;
K1 = getK1(K,m,n);
ksig = (K>0);
phiks=pl(phi-ksig');

Kf = getK1(phiks'.*K,m,n);

phir = recphi(K1,Kf,m,n1);
fprintf('.');

if ~islip(phir)
  'not lip'
  keyboard
end

y = phir*K1 - phiks*K;
return


function y = myobjm_conv(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt
n1 = n-1;

phi = x(1:m^n)';
K = x(m^n+1:end);

K1 = getK1(K,m,n);
ksig = (K>0);
phiks=pl(phi-ksig');

phic = convlips(K1,m,n1)';
fprintf('.');

y = phic*K1 - phiks*K;
return


function y = myobjm_zajim(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt
n1 = n-1;

phi = x(1:m^n)';
K = x(m^n+1:end);

fprintf('.');
global phi00 K00
phi00=phi;K00=K;

K1 = getK1(K,m,n);
ksig = (K>0);
phiks=pl(phi-ksig');
phim = meanphi(phiks,m,n);
%phim = round(phim);
phiz = zajim(K1,phim,m,n1);

y = phiz*K1 - phiks*K;
return


function y = myobjm_posneg(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt
n1 = n-1;

phi = x(1:m^n)';
K = x(m^n+1:end);

Kpl = pl(K);
Kng = pl(-K);
ksig = (K>0);
phiks=pl(phi-ksig');

K1 = getK1(K,m,n);

if phiks*K <= 1e-6
  phib = zeros(1,m^n1);
else
  K1p = getK1(Kpl,m,n);
  K1n = getK1(Kng,m,n);
  
  %phibp = ceil(termwisebd2(Kpl,phiks,m,n));
  %phibn = ceil(termwisebd2(Kng,phi  ,m,n));
  %phibp = termwisebd2(Kpl,phiks,m,n);
  %phibn = termwisebd2(Kng,phi  ,m,n);
  %phibn = meanphi(phiks,m,n);

  phibp = termwisebdu(Kpl,phiks,m,n);
  
  if 0
    phibn = termwisebdl(Kng,phi  ,m,n);
    
    
    %[x,fval] = linprog(-K1p,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
    %phibp = x';
    %[x,fval] = linprog(-K1n,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
    %phibn = x';
    
    phib = pl(phibp - phibn);
    phib = min(phib,n1);
    
    %phib=phibp;
    %d=max(abs(Ain1*phib')) + eps;
    %phib = phib/d;
    
    if ~islip(phib)
      %'not lip'
      phil =  lowelip(phib,m,n1);
      phiu = upperlip(phib,m,n1);
      if (phiu-phil)*K1 > 0
        phib = phiu;
      else
        phib = phil;
      end
    end
  end
  phib = phibp;
end
y = phib*K1 - phiks*K;
return


function y = myobjm_pos(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt
n1 = n-1;

% K is positive
% phi:X^n -> {0,...,n-1}

phi = x(1:m^n)';
K = x(m^n+1:end);

K1 = getK1(K,m,n);

%[x,fval] = linprog(-K1,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
%phi1 = x';

phi1 = termwisebd(K,phi,m,n);

y = phi1*K1 - phi*K;
return
    


function y = myobjm_pos0(x)
global Ainl Binl Ain1 Bin1 GG iif m n
n1 = n-1;

phi = x(1:m^n)';
K = x(m^n+1:end);

Kpl = pl(K);
Kng = pl(-K);
ksig = (K>0);
phiks=pl(phi-ksig');

K1 = getK1(K,m,n);

if phiks*K <= 0
  phib = zeros(1,m^n1);
else

  phibp = termwisebd1(Kpl,phiks,m,n);
  phibn = termwisebd1(Kng,phi,m,n);
  
  phib = pl(phibp - phibn);
  phib = min(phib,n1);
  if ~islip(phib)
    'not lip'
    keyboard
  end
end
y = phib*K1 - phiks*K;
return

function y = myobjm_greedy(x)
global Ainl Binl Ain1 Bin1 GG iif m n
n1 = n-1;

phi = x(1:m^n)';
K = x(m^n+1:end);

ksig = (K>0);
phiks=pl(phi-ksig');

K1 = getK1(K,m,n);
  
phim = meanphi(phiks,m,n);
phim = min(phim,n1);
d = max(Ain1*phim');
%if (d>1)
%  phim = min(phim/d,n1);
%end
phim = phim/d;
phim = phim-min(phim);

%phig = ceil(phim);
%phig = phim;

if 1
  [a,b] = sort(abs(K1));
  %[a,b] = sort(K1);
  b = b(end:-1:1);
  phig = greedymaxf(K1,phim,m,n,b');
  phig = greedymaxf(K1,phig,m,n,b');
  %phig = ceil(phig);
end

y = phig*K1 - phiks*K;
return


function y = myobjm_psi(x)
global Ainl Binl Ain1 Bin1 GG iif m n

phi = x(1:m^n)';
K = x(m^n+1:end);

K1 = getK1(K,m,n);

y = -phi*K;
return

function y = myobjm_phi1(x)
global Ainl Binl Ain1 Bin1 GG iif m n mylinopt
n1 = n-1;

phi = x(1:m^n)';
K = x(m^n+1:end);

ksig = (K>0);
phiks=pl(phi-ksig');
%phiks=pl(phi-1); %bad! BAD!!

K1 = getK1(K,m,n);

[x,fval] = linprog(-K1,Ain1,Bin1,[],[],0*K1,0*K1+n1,0,mylinopt);
phi1 = x';

y = phi1*K1 - phiks*K;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function A = assembleblock(AA)
% takes a cell array of matrices
% and assembles them into a block matrix

nb = length(AA);
mm = zeros(nb,1);
nn = zeros(nb,1);
for i=1:nb
  mm(i) = size(AA{i},1);
  nn(i) = size(AA{i},2);  
end

A = zeros(sum(mm),sum(nn));
row = 0;
col = 0;
for i=1:nb
  A(row+1:row+mm(i),col+1:col+nn(i)) = AA{i};
  row = row + mm(i);
  col = col + nn(i);
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = fitlips(K,Kf,m,n)
%w = n*abs(K)/sum(abs(K));
%w = pl(K)/sum(pl(K));
%w = abs(Kf)/(sum(abs(Kf))+eps);
%w = w/mean(w);

phi = zeros(1,m^n);

iip = find(K>0);

FV = [];
for v=0:n
  phi = multispike(m,n,iip,v);
  FV(end+1) = phi*K;
end

[a,b] = max(FV);
phi = multispike(m,n,iip,b-1);

if 0
  for i=1:m^n
    sig = sign(Kf(i));
    phi = phi + sig*w(i)*lipspike(m,n,i)';
  end
  
  phi = pl(phi);
  phi = min(phi,n);
end


if 0
  [a,b] = max(Kf);
  
  phi = lipspike(m,n,b)';
  
  if sum(Kf)<1e-6
    phi = phi*0;
  end
  
  if all(K>0)
    phi = phi*0 + n;
  end
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function phi = lipspike(m,n,ind)
dims = repmat([m],[1 n]);
phi = zeros(m^n,1);

phi(ind) = n;

for v=n:-1:1
  xind = find(phi==v);
  for i=xind(:)'
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badYi = find(v - phi(Yind) > 1);
    phi(Y(badYi)) = 0*phi(Y(badYi)) + v - 1;
  end
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function Y = allrho1(x,m)
% returns all the Y that are exactly 1 away from x
n = length(x);
dims = repmat([m],[1 n]);
X = list_all_ind(dims);
x = repmat(x,[size(X,1) 1]);
rr = sum(X~=x,2);
yi = find(rr==1);
Y = X(yi,:);
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
% list all n-tuples of the dimension vector v

function list = list_all_ind(v)
list = [];
if ~isempty(v)
  d1 = v(1);

  L1 = list_all_ind(v(2:end));
  [m,n] = size(L1);
  m = max(m,1);
  D1 = repmat(1:d1,[m,1]); D1 = D1(:);
  L2 = repmat(L1,[d1,1]);
  list = [D1,L2];

end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
function ind = coord2ind(coords,dims)
% dims is a ROW vector of dimensions
% coords is a matrix of ROW coordinates vectors
% inds is a list of indices
% this is a recent (16/02/2002) change

% compare with matlab's sub2ind routine

if (length(dims)<size(coords,2))
  dims = repmat(dims,[1,size(coords,2)]);
end

cD = cumprod([1,dims(1:end-1)]);
ind = (coords-1) * cD' + 1;

return


function ind = coord2ind_old(coords,dims)
% NOTE: coords and dims are both ROW vectors

cD = cumprod([1,dims(1:end-1)]);
ind = (coords-1) * cD' + 1;

% REALLY OLD VERSION, to illustrate the idea
%ind = 0;
%coords = coords - 1;
%f = 1;
%for ij = 1:length(dims)
%  ind = ind + coords(ij)*f;
%  f = f*dims(ij);
%end
%ind = ind+1;
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function c = ind2coord(ind,dim)
% dims is a ROW vector of dimensions
% ind is a column vector of indices
% this is a recent (16/02/2002) change

% compare with matlab's ind2sub routine

ind = ind-1;
c = zeros(size(ind,1),size(dim,2));
for ij = 1:length(dim)
  c(:,ij) = mod(ind,dim(ij));
  ind = (ind - c(:,ij))/dim(ij);
end
c = round(c) + 1;
return

function c = ind2coord_old(ind,dim)
ind = ind-1;
c = zeros(size(dim));
for ij = 1:length(dim)
  c(ij) = mod(ind,dim(ij));
  ind = (ind - c(ij))/dim(ij);
end
c = round(c) + 1;
return



%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
function phi = multispike(m,n,inds,val)
dims = repmat([m],[1 n]);
phi = zeros(1,m^n);

phi(inds) = phi(inds) + val;

for v=n:-1:1
  xind = find(phi==v);
  for i=xind(:)'
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badYi = find(v - phi(Yind) > 1);
    phi(Yind(badYi)) = 0*phi(Yind(badYi)) + v - 1;
  end
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function y=pl(x)
y = max(0,x);
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------




%--------------------------------------------------------------
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


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = greedymaxf(K,phi,m,n,perm)
global Ainl Binl Ain1 Bin1
dims = repmat([m],[1 n]);
dims1 = repmat([m],[1 n-1]);

for xind2n = perm
  x2n = ind2coord(xind2n,dims1);
  %ks = 0;
  ks = K(xind2n);
  %for x1 = 1:m
  %  xind = coord2ind([x1 x2n],dims);
  %  ks = ks + K(xind);
  %end
  Y = allrho1(x2n,m);  
  Yind = coord2ind(Y,dims1);
  ymax = max(phi(Yind));
  ymin = min(phi(Yind));
  d = ymax-ymin;
  if ks > 0
    phi(xind2n) = min(ymax + 1 - d,n-1);
  else
    phi(xind2n) = max(ymin - 1 + d,0);
  end
end

return


%--------------------------------------------------------------


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = lowelip(phi,m,n)
dims = repmat([m],[1 n]);

phi = floor(phi);

for val = 0:n-1
  xind = find(phi==val);
  for i=xind(:)'
    v = phi(i);
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badYi = find(phi(Yind) - v > 1);
    phi(Y(badYi)) = 0*phi(Y(badYi)) + v + 1;
  end
end

return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi1 = meanphi(phi,m,n)
global Ainl Binl Ain1 Bin1
dims = repmat([m],[1 n]);
dims1 = repmat([m],[1 n-1]);

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

% init phi1 to be plain-vanilla mean
for xind2n = 1:m^(n-1)
  x2n = ind2coord(xind2n,dims1);
  ff = [];
  for x1 = 1:m
    xind = coord2ind([x1 x2n],dims);
    ff(end+1) = phi(xind);
  end
  phi1(xind2n) = mean(ff);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function S = Psi(K,m,n)
S = 0;
for t=1:n
  S = S + sum(K.*(K>0));
  K = getK1(K,m,n-t+1);
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = randlip(m,n)
[A,B] = lipab(m,n);
GG = gengray(n,m);
iif = invperm(GG+1);

while 1
  phi = zeros(m^n,1);
  rp = randperm(m^n);
  diam = ceil(rand*n);
  %diam = n;
  rp = rp(1:diam);
  phi(rp) = ones(diam,1);
  phi = cumsum(phi);
  phi = phi(iif);
  if all(A*phi<=B)
    break
  end
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function d=diam(X)
d=0;
for i=1:length(X)
  for j=i+1:length(X)
    %xi=X{i}
    %xj=X{j}
    %rij = rho(X{i},X{j})
    d = max(d,rho(X{i},X{j}));
  end
end
return


function r=rho(x,y)
r = length(find(x~=y));
return



%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function GG = gengray(N,K)
% generate the N-ary Gray code of length K
% taken from http://nr.stic.gov.tw/ejournal/ProceedingA/v22n6/841-848.pdf
% http://en.wikipedia.org/wiki/Gray_code
% Guan, Dah-Jyu. "Generalized Gray Codes with Applications." Proc. Katl. Sci. Counc. Repub. Of China (A) 22 (1998): 841-848. [5]. 

GG = [];
g = [];
u = [];
n = [];
for i=0:N
    g(i+1) = 0;
    u(i+1) = 1;
    n(i+1) = K;
end

while ~g(N+1)
    %for j=N-1:-1:0
    %    fprintf('%d',g(j+1));
    %end
    str = '';
    for j=N-1:-1:0
        str = [str, num2str(g(j+1))];
    end
    %fprintf('%s\n',str);
    GG(end+1) = base2dec(str,K);
    
    i = 0;
    k = g(1) + u(1);
    while ( (k>=n(i+1)) | (k<0) )
        u(i+1) = -u(i+1);
        i = i+1;
        k = g(i+1) + u(i+1);
    end
    g(i+1) = k;
end

return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function Q=invperm(P)
n = length(P);
Q = zeros(size(P));
for i=1:n
    Q(i) = find(P==i);
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------


%--------------------------------------------------------------
function p = P(x,A,C)
x=x+1;
if isempty(x)
  p = 1;
else
  p = C(x(1));
end
for j=2:length(x)
  p = p * A(x(j),x(j-1)); % markov case
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------





%--------------------------------------------------------------
function phi = recphi(K,Kf,m,n)
if n==1
  % base case
  %if 1
  %  if sum(K)>0
  %    phi = ones(1,m);
  %  else
  %    phi = zeros(1,m);
  %  end
  %else
  if all(K>=0)
    phi = ones(1,m);
  elseif all(K<0)
    phi = zeros(1,m);
  else
    phi = pl(K');
    phi = phi/(max(phi)+eps);
    %  phi = zeros(1,m);
    %  iip = find(K>1e-6);
    %  phi(iip) = ones(size(phi(iip)));
    %  %iin = find(K<0);
    %  %if 
    %  %[a,b] = max(K);
    %  %phi = zeros(1,m);
    %  %phi(b) = 1;
  end
else
  phi = zeros(1,m^n);
  n1 = n-1;
  dims  = repmat([m],[1 n ]);
  dims1 = repmat([m],[1 n1]);
  Y = list_all_ind(dims1);
  %K1 = getK1(K,m,n);
  %phi1 = recphi(K1,m,n1);
  for x1 = 1:m
    x1Y = [repmat([x1],[size(Y,1) 1]) Y];
    xyi = coord2ind(x1Y,dims);
    K1 = K(xyi);
    Kf1 = Kf(xyi);
    phi1 = recphi(K1,Kf1,m,n1);
    phi(xyi) = phi1;
    %if phi1*K1 > 0
    %  phi(xyi) = phi1 + 1;
    %else
    %  phi(xyi) = phi1;
    %end
    %phi(xyi) = phi1;
    %if sum(K(xyi))>0
    %  phi(xyi) = phi(xyi) + 1;
    %else
    %  phi(xyi) = pl(phi(xyi) - 1);
    %end
  end    
  
  
  phil = pl(phi-1);
  if phil*K > phi*K
    phi = phil;
  end
  
  phih = min(phi+1,n);
  if phih*K > phi*K
    phi = phih;
  end
  
  %if sum(K)<0
  %  phi = 
  %else
  %  phi = min(phi+1,n);
  %end
end
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phib = termwisebd(K,phi,m,n)
n1 = n-1;

global Ainl Binl Ain1 Bin1 GG iif m n mylinopt


% assume K>0

K1 = getK1(K,m,n);
Kf = getK1(phi'.*K,m,n);

phib = zeros(m^n1,1);

if sum(K1)>1e-6
  phib = min(Kf ./ (K1 + ~K1),n1);
end

%phib = ceil(phib);

dims = repmat([m],[1 n1]);

for val = n1:-1:1
  %xind = find(phib==val);
  xind = find(ceil(phib)==val);
  for i=xind(:)'
    v = phib(i);
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badY = find(v-phib(Yind) > 1);
    phib(badY) = 0*phib(badY) + v - 1;
  end
end

phib = phib';
return


function phib = termwisebd0(K,phi,m,n)
n1 = n-1;

global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

if isempty(mylinopt)
  mylinopt = optimset('Diagnostics','off');
  mylinopt = optimset('Display','off');
  [Ain1,Bin1] = lipab(m,n1);
end

%ksig = (K>0);
%phiks=pl(phi-ksig');

K1 = getK1(K,m,n);
Kf = getK1(phi'.*K,m,n);

phib = zeros(m^n1,1);

LB = zeros(size(phib));
%LB = LB - Inf;
iin = find(K1<0);
LB(iin) = Kf(iin)./(K1(iin) + ~K1(iin));

UB = n1+zeros(size(phib));
%UB = UB + Inf;
iip = find(K1>0);
UB(iip) = Kf(iip)./(K1(iip) + ~K1(iip));

[x,fval] = linprog(-K1,Ain1,Bin1,[],[],LB,UB,0,mylinopt);
phib = x';

Kb = phib'.*K1;
marg = min(Kb - Kf);
if marg < -1e-4
  'problem'
  keyboard
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phib = termwisebd1(K,phi,m,n)
n1 = n-1;

global Ainl Binl Ain1 Bin1 GG iif m n mylinopt

if isempty(mylinopt)
  mylinopt = optimset('Diagnostics','off');
  mylinopt = optimset('Display','off');
  [Ain1,Bin1] = lipab(m,n1);
end

% assume K>0

K1 = getK1(K,m,n);
Kf = getK1(phi'.*K,m,n);

phib = zeros(m^n1,1);

if sum(K1)<1e-12
  phib = phib';
  return
end

LB = Kf ./ (K1 + ~K1);

UB = n1+zeros(size(phib));

[x,fval] = linprog(K1,Ain1,Bin1,[],[],LB,UB,0,mylinopt);
phib = x';

Kb = phib'.*K1;
marg = min(Kb - Kf);
if marg < -1e-4
  'problem'
  keyboard
end
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phib = termwisebdl(K,phi,m,n)
n1 = n-1;

global Ainl Binl Ain1 Bin1 GG iif m n mylinopt


% assume K>0

K1 = getK1(K,m,n);
Kf = getK1(phi'.*K,m,n);

phib = zeros(m^n1,1);

if sum(K1)>1e-6
  phib = min(Kf ./ (K1 + ~K1),n1);
end

phib = floor(pl(phib));

dims = repmat([m],[1 n1]);

for val = 0:n1-1
  xind = find(phib==val);
  for i=xind(:)'
    v = phib(i);
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badYi = find(phib(Yind) - v > 1);
    phib(Y(badYi)) = 0*phib(Y(badYi)) + v + 1;
  end
end

phib = phib';
return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phib = termwisebdu(K,phi,m,n)
n1 = n-1;

global Ainl Binl Ain1 Bin1 GG iif m n mylinopt


% assume K>0

K1 = getK1(K,m,n);
Kf = getK1(phi'.*K,m,n);

phib = zeros(m^n1,1);

if sum(K1)>1e-6
  phib = min(Kf ./ (K1 + ~K1),n1);
end

phib = ceil(phib);

dims = repmat([m],[1 n1]);

for val = n1:-1:1
  xind = find(phib==val);
  %xind = find(ceil(phib)==val);
  for i=xind(:)'
    v = phib(i);
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badYi = find(v-phib(Yind) > 1);
    phib(Y(badYi)) = 0*phib(Y(badYi)) + v - 1;
  end
end

phib = phib';
return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = upperlip(phi,m,n)
dims = repmat([m],[1 n]);

phi = ceil(phi);

for val = n:-1:1
  xind = find(phi==val);
  %xind = find(ceil(phi)==val);
  for i=xind(:)'
    v = phi(i);
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    badYi = find(v-phi(Yind) > 1);
    phi(Y(badYi)) = 0*phi(Y(badYi)) + v - 1;
  end
end

return


%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------



%--------------------------------------------------------------
function phi = zajim(K,phi,m,n)
dims = repmat([m],[1 n]);

iip = find(K>0);
iin = find(K<0);

phi(iip) = ceil(phi(iip));
phi(iin) = floor(phi(iin));

%phi = zeros(1,m^n);
%phi(iip) = phi(iip) + n;

x = abs(K);
[a,b] = sort(x);
%b = b(end:-1:1);

while ~islip(phi)
  %x = abs(phi'.*K);
  %[a,b] = sort(x);
  %b = b(end:-1:1);
  
  for i=b(:)'
    v = phi(i);
    x = ind2coord(i,dims);
    Y = allrho1(x,m);
    Yind = coord2ind(Y,dims);
    if K(i)>0
      % its K is positive
      % bring it down to its lowest neighbor
      w = min(phi(Yind));
      if v-w > 1
        phi(i) = max(v-1,w);
      end
    else
      % its K is negative
      % bring it up to its highest neighbor
      w = max(phi(Yind));
      if w-v > 1
        phi(i) = min(v+1,w);
      end
    end
  end
end

b = b(end:-1:1);
phi = greedymaxf(K,phi,m,n,b(:)');

return


function phi = greedymaxf(K,phi,m,n,perm)
global Ainl Binl Ain1 Bin1
dims = repmat([m],[1 n]);

for xind2n = perm
  x2n = ind2coord(xind2n,dims);
  ks = K(xind2n);
  Y = allrho1(x2n,m);  
  Yind = coord2ind(Y,dims);
  ymax = max(phi(Yind));
  ymin = min(phi(Yind));
  d = ymax-ymin;
  if ks > 0
    phi(xind2n) = min(ymax + 1 - d,n);
  else
    phi(xind2n) = max(ymin - 1 + d,0);
  end
end

return


function z0
for i=b(:)'
  v = phi(i);
  x = ind2coord(i,dims);
  Y = allrho1(x,m);
  Yind = coord2ind(Y,dims);
  if K(i)>0
    % its K is positive
    % bring its neighbors up to it
    % find where it exceeds its neighbors
      badYi = find(v - phi(Yind) > 1);
      phi(Y(badYi)) = 0*phi(Y(badYi)) + v - 1;
  else
    % its K is negative
    % bring its neighbors down to it
    badYi = find(phi(Yind) - v > 1);
    phi(Y(badYi)) = 0*phi(Y(badYi)) + v + 1;
  end
end

return

%-----------------------FUNCTIONS INCLUDED FOR SELF-SUFFICIENCY-----------------------