subroutine clusterfuncf (mat,ki,kj,ccenters,parentcluster,
     + totclusters,ninclusters,maxsize)
      implicit none
      integer i,j,ki,kj,ninclusters(ki),
     + parentcluster(ki),totclusters,t,di,maxsize
      double precision ccenters(ki,maxsize),mat(ki,maxsize),
     + dlimit,distance,tempdist,nspace


      nspace = kj/1D0


C     start with first observation
      do 30 j = 1,kj
         ccenters(1,j) = mat(1,j)
 30   continue
      ninclusters(1) = 1
      totclusters = 1
      dlimit = sqrt(nspace)*.1D0
      parentcluster(1) = 1
      
      
      do 40 i = 2,ki
         distance = sqrt(kj/1D0)
         do 50 t = 1,totclusters
            tempdist = 0
            do 55 j = 1,kj
              tempdist = tempdist + (ccenters(t,j) -mat(i,j))*
     +                              (ccenters(t,j) -mat(i,j))
 55         continue
            tempdist = sqrt(tempdist)
            if (tempdist .lt. distance) then
               distance = tempdist
               di = t
            end if
 50      continue
         if (distance .lt. dlimit) then 
            ninclusters(di) = ninclusters(di)+1
            do 60 j =1,kj
               ccenters(di,j) = (ccenters(di,j)*
     +          (ninclusters(di)-1)+mat(i,j))/ninclusters(di)
 60         continue
            parentcluster(i) = di
         else 
            totclusters = totclusters + 1
            ninclusters(totclusters) = 1
            parentcluster(i) = totclusters
            do 70 j = 1,kj
               ccenters(totclusters,j) = mat(i,j)
 70         continue
         end if
 40   continue
         
      return
      end

      

      subroutine fitnessfuncf (mat,ki,kj,ccenters,parentcluster,
     + totclusters,ninclusters,class01s,n0inclusters,purity,purities,
     + maxsize)
      implicit none
      integer totclusters, ki, kj, parentcluster(ki), 
     + ninclusters(ki),n0inclusters(ki),class01s(ki),i,di,maxsize
      double precision ccenters(ki,maxsize), mat(ki,maxsize), 
     + purities(ki), purity


      call clusterfuncf(mat,ki,kj,ccenters,parentcluster,
     + totclusters,ninclusters,maxsize)

      do 10 i = 1,ki
         if (class01s(i) .eq. 0) then
            n0inclusters(parentcluster(i)) = 
     +      n0inclusters(parentcluster(i)) + 1 
         end if
 10   continue
      
      purity = 0
      do 20 di = 1,totclusters
         purities(di) = max(1d0*n0inclusters(di)/ninclusters(di),
     +   1-1d0*n0inclusters(di)/ninclusters(di))
         purity = purity+purities(di)*ninclusters(di)/ki
 20   continue

      return
      end
      
      


      subroutine fitnessmatf (rset1,fitnessmat,chromes,chromesizes,
     +     nchromes,ki,maxj,maxsize,class01s)
      implicit none
      integer i,ki,j,maxsize,parentcluster(ki), ninclusters(ki), 
     + n0inclusters(ki),class01s(ki),totclusters,n,nchromes,kj,
     + chromesizes(nchromes),chromes(nchromes,maxsize),
     + epoints(maxsize),maxj
      double precision holdme(ki,maxsize), ccenters(ki,maxsize),
     + purities(ki), purity,maxrow,minrow,rset1(ki,maxj),
     + fitnessmat(nchromes,2)

      do 10 n = 1,nchromes
         do 3 i=1,ki
            do 5 j = 1,maxsize
               holdme(i,j) = 0
               ccenters(i,j) = 0
    5       continue
            parentcluster(i) = 0
            ninclusters(i) = 0
            n0inclusters(i) = 0
            purities(i) = 0
 3       continue
         totclusters = 0
         purity = 0
         kj = chromesizes(n)
         do 20 j = 1,kj
            epoints(j) = chromes(n,j)
 20      continue
         
C     compute normalization matrix in main loop, not function
C     or subroutine
C     holdme is name of normalized matrix

        do 30 i = 1,ki
           maxrow = 0
           minrow = 1000
           do 40 j = 1,kj
              holdme(i,j) = rset1(i,epoints(j))
              if (holdme(i,j) .gt. maxrow) then
                 maxrow = holdme(i,j)
              end if
              if (holdme(i,j) .lt. minrow) then
                 minrow = holdme(i,j)
              end if
   40      continue
           if (minrow .eq. maxrow) then
              do 50 j = 1,kj
                 holdme(i,j) = .5
   50         continue
           else
              do 60 j = 1,kj
                 holdme(i,j) = (holdme(i,j)-minrow)/(maxrow-minrow)
   60         continue
           endif
 30     continue
         
        call fitnessfuncf(holdme,ki,kj,ccenters,parentcluster,
     +  totclusters,ninclusters,class01s,n0inclusters,purity,
     +  purities,maxsize)


        fitnessmat(n,1) = purity
        fitnessmat(n,2) = totclusters*1D0
 10   continue
      return
      end