function [permut,cumobfun] = class_scaledp(prox,numbclass,membclass)

%  CLASS_SCALEDP carries out a unidimensional seriation or 
%  scaling of a set of object classes defined for a symmetric proximity
%  matrix using dynamic programming.
%
%  syntax: [permut,cumobfun] = class_scaledp(prox,numbclass,membclass)
%
%  PROX is the ($n \times n$) input proximity matrix (with a zero 
%  main diagonal and a dissimilarity interpretation);
%  NUMBCLASS (= $n_{c}$) is the number of object classes to be
%  sequenced; 
%  MEMBCLASS is an $n \times 1$ vector containing the input class 
%  membership and includes all the integers from 1 to NUMBCLASS and
%  zeros when objects are to be deleted from consideration;
%  PERMUT is the order of the classes 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_{c}$.
%
%  Initializations:  The vectors VALSTORE and IDXSTORE store the
%  results of the recursion for the $(2^n_{c})-1$ nonempty subsets of the
%  set of classes. The integer positions in these vectors correspond to 
%  subsets whose binary number equivalents are equal to those integer positions.

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

member = zeros(n,n);
ncsize = zeros(nc,1);

for i = 1:n
    for j = 1:nc
        
        if(membclass(i) == j)
            
            ncsize(j) = ncsize(j) + 1;
            member(j,ncsize(j)) = i;
            
        end
    end
end

%  Initializing the recursion for subsets of size 1

for i = 1:nc
   rowsum = 0;
   for j = 1:nc
       
       for h2 = 1:ncsize(j)
           for h1 = 1:ncsize(i)
               
               if(i ~= j)
                   
                   rowsum = rowsum + prox(member(i,h1),member(j,h2));
                   
               end
           end
       end
                   
    end
   index = 2^(i-1);
   idxstore(index) = i;
   valstore(index) = (rowsum^2)/ncsize(i);
end

%  Carrying out the recursion for subsets of size $k$ from $k = 1$ 
%  to $n_{c}-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:(nc-1)
   nfirst = 0;
   if(mtc == 0)
      m2 = 0;
      nh = k;
      for j = 1:nh
         sub(k+j-nh) = m2 + j;
      end
      if(sub(1) ~= (nc-k+1))
         mtc = 1;
      end
   end
 
   while (mtc == 1)
      if(nfirst == 1)
         if(m2 < (nc-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) ~= (nc-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(nc,1);
      for i = 1:nc
         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 classes from SUBCOMP one at a time to the end of
 %  SUB to see if better placements can be found  
 
      nk = nc - k;
      for jj = 1:nk
         jone = subcomp(jj);
         sum = 0.0;
         for i = 1:k
             for l1 = 1:ncsize(sub(i))
                 for l3 = 1:ncsize(jone)
            ione = sub(i);
            sum = sum + prox(member(jone,l3),member(ione,l1));
                 end
             end
         end
         for i = 1:nk
             for l2 = 1:ncsize(subcomp(i))
                 for l3 = 1:ncsize(jone)
                     ione = subcomp(i);
            if(ione ~= jone)
               sum = sum - prox(member(jone,l3),member(ione,l2));
            end
                 end
             end
         end
         
 %  Evaluating the placement of the classes from SUBCOMP to
 %  the end of those in SUB 
 
         incre = (sum^2)/ncsize(jone);
         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(nc) = idxstore(nsum);
cumobfun(nc) = valstore(nsum);
index = nsum;
lastint = permut(nc);
for i = 1:(nc-1)
   index = index - (2^(lastint-1));
   lastint = idxstore(index);
   permut(nc-i) = lastint;
   cumobfun(nc-i) = valstore(index);
end

   
      
      
      
      
      
      
      
  


          


