function [sols, mbs, d] = bicriterion3(a, c)

%BICRITERION3 performs bicriterion clustering analysis.
%          PRIMARY OBJECTIVE: Minimize partition diameter using
%                             exchange algorithm
%          SECONDARY OBJECTIVE: Minimize within-cluster sum of dissimilarity
%                               (normalized by dividing by n) using the
%                               exchange algorithm.
% INPUTS:  a = an n x n dissimilarity matrix
%          c = a desired number of clusters
% OUTPUTS: sols =  a matrix of index values
%             col. 1 = of diameter values for each replication
%             col 2 = a vector of within cluster sums (WCS) for each rep.
%             col 3 = a vector of WCS divided by n
%             col 4 = a vector of WCS divided by n(n-1)/2
%          mbs = a matrix containing membership assgnmts for each rep.
%          d = a matrix of cluster diameters.

eps=.0000001;
tic;
n=size(a,1);                   % a is an n x n dissimilarity matrix
nreps = 100;                    % Total number of replications
sols = zeros(nreps,4);         % Matrix of index values for each rep.
mbs = zeros(nreps, n);         % Matrix of group memberships.
d = zeros(nreps, c);           % Matrix of diameters.
nca = zeros(1,n);              % number in cluster

k = 0;
for i = 1:n
  k = k + 1;
  if k > c
    k = 1;
  end
  memb(i)=k;                   % Membership of object i is cluster k
  nca(k)=nca(k)+1;             % Number of objects in cluster k
end

for iii = 1:nreps              % Perform nreps repetitions of SA
  rp = randperm(n);            % Select a random assignment of objects
  p = memb(rp);                % and let p be the initial memberships
  for k = 1:c
    s(k)=0;
    nc(k)=nca(k);              
  end
   
  trig = 1;
  while trig == 1
    trig = 0;
    trig1 = 1;
    while trig1 == 1
      trig1=0;
      for i = 1:n
        k1=p(i);
        if nc(k1) <= 1
          continue
        end
        for k = 1:c
          if k == k1
            continue
          end
          delt1=0;
          delt2=0;
          for j = 1:n
            if j == i
              continue
            end
            k2 = p(j);
            if k2 == k
              if a(i,j) > delt2
                delt2 = a(i,j);
              end
            end
            if k2 == k1
              if a(i,j) > delt1
                delt1 = a(i,j);
              end
            end
          end

          if delt2 < delt1
            nc(k1) = nc(k1) - 1;
            nc(k) = nc(k) + 1;
            p(i)=k;
            trig1 = 1;
            trig = 1;
            k1 = k;
          end
        end
      end
    end

    trig2 = 1;
    while trig2 == 1;
      trig2 = 0;
      for i = 1:n-1
        k1 = p(i);
        for j = i+1:n
          k2=p(j);
          if k1 == k2
            continue
          end
          delt1=0;
          delt2=0;
          for m = 1:n
            if m == i || m == j
              continue
            end
            k3=p(m);
            if k3 == k1
              if a(m,j) > delt2
                delt2 = a(m,j);
              end
              if a(m,i) > delt1
                delt1 = a(m,i);
              end
            end
            if k3 == k2
              if a(m,i) > delt2
                delt2 = a(m,i);
              end
              if a(m,j) > delt1
                delt1 = a(m,j);
              end
            end
          end

          if delt2  < delt1
            p(i)=k2;
            p(j)=k1;
            trig2 = 1;
            trig = 1;
            k1 = k2;
          end
        end
      end
    end

  end


  for i = 1:c
    nc(i)=0;
    s(i)=0;
  end
  for i = 1:n
    k1=p(i);
    nc(k1)=nc(k1)+1;
  end
  diamb=0;
  for i = 1:n-1
    k1 = p(i);
    for j = i+1:n
      k2 = p(j);
      if k1 == k2
        s(k1) = s(k1) + a(i,j);
        if a(i,j) > diamb
          diamb = a(i,j);
        end
      end
    end
  end
  sumc = 0;
  for k = 1:c
    sumc = sumc + s(k) ./ nc(k);
  end

  trig = 1;
  while trig == 1
    trig = 0;
    trig1 = 1;
    while trig1 == 1
      trig1=0;
      for i = 1:n
        k1=p(i);
        if nc(k1) <= 1
          continue
        end
        for k = 1:c
          if k == k1
            continue
          end
          delt1=0;
          delt2=0;
          infeas=0;
          for j = 1:n
            if j == i
              continue
            end
            k2 = p(j);
            if k2 == k
              if a(i,j) > diamb
                infeas = 1;
                break
              end
              delt2 = delt2 + a(i,j);
            end
            if k2 == k1
              delt1 = delt1 + a(i,j);
            end
          end
          if infeas == 1
            continue
          end
          
          delta = (s(k1)-delt1)./(nc(k1)-1) + (s(k)+delt2)./(nc(k)+1) - s(k1)./nc(k1) - s(k)./nc(k);
          if delta < -eps
            sumc = sumc + delta;
            s(k1) = s(k1) - delt1;
            s(k) = s(k) + delt2;
            nc(k1) = nc(k1) - 1;
            nc(k) = nc(k) + 1;
            p(i)=k;
            trig1 = 1;
            trig = 1;
            k1 = k;
%            break
          end
        end
      end
    end

    trig2 = 1;
    while trig2 == 1;
      trig2 = 0;
      for i = 1:n-1
        k1 = p(i);
        for j = i+1:n
          k2=p(j);
          if k1 == k2
            continue
          end
          delt1=0;
          delt2=0;
          infeas = 0;
          for m = 1:n
            if m == i || m == j
              continue
            end
            k3=p(m);
            if k3 == k1
              if a(m,j) > diamb
                infeas = 1;
                break
              end
              delt1 = delt1 + a(m,j) - a(m,i);
            end
            if k3 == k2
              if a(m,i) > diamb
                infeas = 1;
                break
              end
              delt2 = delt2 + a(m,i) - a(m,j);
            end
          end
          if infeas == 1
            continue
          end

          delta = (s(k1)+delt1)./nc(k1) + (s(k2)+delt2)./nc(k2) - s(k1)./nc(k1) - s(k2)./nc(k2);
          if delta < -eps
            sumc = sumc + delta;
            s(k1) = s(k1) + delt1;
            s(k2) = s(k2) + delt2;
            p(i)=k2;
            p(j)=k1;
            trig2 = 1;
            trig = 1;
            k1 = k2;
          end
        end
      end
    end

  end

  diamb=0;
  for i = 1:c
    s(i)=0;
    nc(i)=0;
  end
  for i = 1:n
    k1=p(i);
    nc(k1)=nc(k1)+1;
  end
  for i = 1:n-1
    k1 = p(i);
    for j = i+1:n
      k2 = p(j);
      if k1 == k2
        s(k1) = s(k1) + a(i,j);
        if a(i,j) > diamb
          diamb = a(i,j);
        end
        if a(i,j) > d(iii,k1)
          d(iii,k1) = a(i,j);
        end        
      end
    end
  end
  sumc1 = 0;
  sumc2 = 0;
  sumc3 = 0;
  for k = 1:c
    sumc1 = sumc1 + s(k);
    sumc2 = sumc2 + s(k) ./ nc(k);
    if nc(k) > 1
      sumc3 = sumc3 + s(k) ./ (nc(k).*(nc(k)-1)./2);
    end
  end    
  sols(iii,1)=diamb;
  sols(iii,2)=sumc1;
  sols(iii,3)=sumc2;
  sols(iii,4)=sumc3;
  mbs(iii,:)=p;

end
toc