function [coord,permut,cumobfun,diff] = uniscaldp(prox)
%  UNISCALDP carries out a unidimensional scaling of a symmetric proximity
%  matrix using dynamic programming.
%  PROX is the input proximity matrix (with a zero main diagonal and a 
%  dissimilarity interpretation);
%  COORD is the set of coordinates of the optimal unidimensional scaling 
%  in ascending order;
%  PERMUT is the order of the objects in the optimal permutation (say,
%  $\rho^{*}$);
%  CUMOBFUN gives the cumulative values of the objective function for
%  the successive placements of the objects in the optimal permutation:
%  $\sum_{i=1}^{k} (t_{i}^{(\rho^{*})})^{2}$ for $k = 1, \ldots, n$;
%  DIFF is the value of the least-squares loss function for the optimal
%  coordinates and object permutation.

%  Initializations.  The vectors VALSTORE and IDXSTORE store the results of 
%  of the recursion for the $(2^n)-1$ nonempty subsets of $S$. The integer
%  positions in these vectors correspond to subsets whose binary
%  number equivalents are equal those integer positions.

n = size(prox,1);
bigneg = -1.0e+20;
nsum = (2^n) - 1;
valstore = ones(nsum,1)*bigneg;
idxstore = zeros(nsum,1);
sub = zeros(n,1);
subcomp = zeros(n,1);
coord = zeros(n,1);
permut = zeros(n,1);
cumobfun = zeros(n,1);

%  Initializing the recursion for subsets of size 1
for i = 1:n
   rowsum = 0;
   for j = 1:n
      rowsum = rowsum + prox(i,j);
   end
   index = 2^(i-1);
   idxstore(index) = i;
   valstore(index) = rowsum^2;
end

%  Carrying out the recursion for subsets of size $k$ from $k = 1$ 
%  to $n-1$ by adapting a Fortran subroutine from Nijenhuis and Wilf
%  (1978, Combinatorial Algorithms, Academic Press, pp. 26--38) on 
%  generating the next k-subset from an n-Set
%  (and placed in the variable SUB)
mtc = 0;
for k = 1:(n-1)
   nfirst = 0;
   if(mtc == 0)
      m2 = 0;
      nh = k;
      for j = 1:nh
         sub(k+j-nh) = m2 + j;
      end
      if(sub(1) ~= (n-k+1))
         mtc = 1;
      end
   end
 
   while (mtc == 1)
      if(nfirst == 1)
         if(m2 < (n-nh))
            nh = 0;            
         end
         nh = nh + 1;
         m2 = sub(k+1-nh);
         for j = 1:nh
            sub(k+j-nh) = m2 + j;
         end
         if(sub(1) ~= (n-k+1))
            mtc = 1;
         else
            mtc = 0;
         end
      end
      
      if(nfirst == 0)
         nfirst = 1;
      end
         index = 0;
      for i = 1:k
         index = index + 2^(sub(i)-1);
      end
      
%  Generating the complement (SUBCOMP) of the subset indicated in SUB      
      jj = 1;
      subcomp = zeros(n,1);
      for i = 1:n
         idxcomp = 1;
         for j = 1:k
            if(sub(j) == i)
               idxcomp = 0;
            end
         end
         if(idxcomp == 1)
            subcomp(jj) = i;
            jj = jj + 1;
         end
      end
      
 %  Trying to add objects from SUBCOMP one at a time to the end of
 %  SUB to see if better placements can be found     
      nk = n - k;
      for jj = 1:nk
         jone = subcomp(jj);
         sum = 0.0;
         for i = 1:k
            ione = sub(i);
            sum = sum + prox(jone,ione);
         end
         for i = 1:nk
            ione = subcomp(i);
            if(ione ~= jone)
               sum = sum - prox(jone,ione);
            end
         end
         
 %  Evaluating the placement of the objects from SUBCOMP to
 %  the end of those in SUB   
         incre = sum^2;
         temp = valstore(index) + incre;
         idxtmp = index + 2^(jone-1);
         comp = valstore(idxtmp);
         if(temp > comp)
            valstore(idxtmp) = temp;
            idxstore(idxtmp) = jone;
         end
      end
   end
end

%  Recursion complete;  working backwards to find the optimal solution
permut(n) = idxstore(nsum);
cumobfun(n) = valstore(nsum);
index = nsum;
lastint = permut(n);
for i = 1:(n-1)
   index = index - (2^(lastint-1));
   lastint = idxstore(index);
   permut(n-i) = lastint;
   cumobfun(n-i) = valstore(index);
end

for i = 1:n
   for j = 1:n
      if(i > j)
         coord(i) = coord(i) + prox(permut(i),permut(j));
      end
      if(i < j) 
         coord(i) = coord(i) - prox(permut(i),permut(j));
      end
   end
   coord(i) = coord(i)/n;
end

diff = 0;
for i = 1:(n-1)
   for j = (i+1):n
      diff = diff + (prox(permut(i),permut(j)) - ...
         abs(coord(i) - coord(j)))^2;
   end
end

        
      
      
      
      
      
      
      
  


          


