function ksipsi(m0,n0)
global m n
randomize
m = m0;
n = n0;

global pstr
pstr = mfilename
global lstr rstr


global Ainf Binf mylinopt
mylinopt = optimset('Diagnostics','off');
mylinopt = optimset('Display','off');
[Ainf,Binf] = lipab(m,n);  

Ain = -ones(1,m^n);
Bin = -1e-6;

global maxr mind
%maxr = 0;
mind = 1e10;
while 1
  x = randn(m^n,1); 
  if sum(x)<0
    x = -x;
end

  y = fmincon(@myobj,x,Ain,Bin,[],[],0*x-1e6,0*x+1e6);
  %y = fminunc(@myobj,x);
  
  fprintf('%s: PARAMS: n=%d m=%d \n',pstr,n,m);
  %fprintf('RAT _________________________ maxr = %1.9f \n',maxr);
  fprintf('DIF _________________________ mind = %1.9f \n',mind);
end

return


function y = myobj(x)
global m n
global pstr
if sum(x)<0
    x = -x;
end


lhs = Ksi(x);
%rhs = Psi(-x,m,n);
rhs = npsi(x,m,n);
%rhs = sum(pl)
%rhs = sum(x);

y = lhs - rhs;

%r = rhs/(lhs+eps);
%y = -r;
global maxr mind
%if (r > maxr)
if (y < mind)
  global optx
  optx = x;
  %maxr = r;
  mind = y;
  %if y<-10
  %    'bla'
  %end
  fprintf('%s: PARAMS: n=%d m=%d \n',pstr,n,m);
  %fprintf('RAT _________________________ maxr = %1.9f \n',maxr);
  fprintf('DIF _________________________ mind = %1.9f \n',mind);
end

return


function [y,phi] = Ksi(K)
global Ainf Binf m n mylinopt
try
[phi,y] = linprog(K,Ainf,Binf,[],[],0*K+1,0*K+1e5,0,mylinopt);
catch
[phi,y] = linprog(K,Ainf,Binf,[],[],0*K+1,0*K+1e2,0,mylinopt);
end
%phi = phi';
return

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


%--------------------------------------------------------------
function [Ainf,Binf] = lipab(m,n)
d = 1;
mn = m^n;
L = mn;
Ainf = zeros(L,mn);
Binf = zeros(L,1);
MN = repmat([m],[1 n]);
rind = 0;
for ind1 = 1:mn
  for ind2 = ind1+1:mn
    x = ind2coord(ind1,MN);
    y = ind2coord(ind2,MN);
    if length(find(x~=y))==1  % Hamming dist == 1
      for s = [-1,1]
        rind = rind+1;
        Ainf(rind,ind1) =  s;
        Ainf(rind,ind2) = -s;
        Binf(rind,1) = d;
      end
    end
  end
end
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 [S,SS] = npsi(K,m,n)
S = 0;
for t=1:n
  SS(t) = sum(K.*(K<0));
  S = S + SS(t);
  K = getK1(K,m,n-t+1);
end
%S = S+pl(K);
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 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-----------------------





%--------------------------------------------------------------
n = round(randinseg(10,20));
x = now;
for t=1:n
  x = 1/(x-floor(x));
end
x = round((x-floor(x))*1e9);
rand('seed',x)


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


%--------------------------------------------------------------
function r = randinseg(a,b)
r = rand*(b-a) + a;
return


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