FuzzyKMeanAlg.cxx

Go to the documentation of this file.

///
/// \file    FuzzyKMeanAlg.cxx
/// \brief   Fuzzy k-means vertex clustering algorithm
/// \author  eninerr@indiana.edu
/// \version $Id: FuzzyKMeanAlg.cxx,v 1.5 2012-10-01 00:29:02 edniner Exp $
///
#include "FuzzyKVertex/FuzzyKMeanAlg.h"

#include "RecoBase/Cluster.h"
#include "Utilities/func/MathUtil.h"

#include <iostream>
#include <cmath>


namespace fuzz{

  FuzzyKMeanAlg::FuzzyKMeanAlg(const FuzzyKMeanParams& params)
    : fParams(params)
  {
  }
  
  //......................................................................
  
  FuzzyKMeanAlg::~FuzzyKMeanAlg()
  {
  }

  //......................................................................

  ///Function to load in a cluster and build the hitlist
  void FuzzyKMeanAlg::LoadCluster(rb::Cluster slice)
  {
    fHL.clear();
    fHL = rb::HitList(slice, fParams.WeightByQ()); //initialize hitlist
  }
  
  //......................................................................
  
  ///Function performs clustering using a fuzzy possibilistic clustering
  ///approach.  Takes arguements for the minimum number of clusters to start
  ///with and the maximum number of clusters to allow. Unused hits will be
  ///left unclustered.  Other arguments are a vertex position and geometry view
  void FuzzyKMeanAlg::DoClustering(unsigned int minClust, unsigned int maxClust, 
                                double xy, double z, geo::View_t view)
  {
    this->MakeAngles(xy, z, view);  //turn hits into angles around vertex
    fA.resize(minClust);
    fSeed.clear();
    this->SeedA(minClust); //produce intitial cluster seeds

    //loop to increase number of clusters
    while (fA.size() <=maxClust){
      //optimize cluster centers and assign membership for set cluster number
      this->MakeClusters(fA.size());

      //now check that all clusters centers produced are unique
      //remove duplicates and recluster
      if (this->CheckUnique() == true) continue;

      //check to see if there are cells without cluster membership
      //stop clustering if no unassigned hits are found
      if (this->FindUnclustered() == false) break;
    }//end while
    
    //remove a cluster if it excedes the limit
    if (fA.size() > maxClust) fA.pop_back();

    //now check to see if any clusters should be merged
    this->MergeClusters();
    
    //print statements for debugging
    if (fParams.PrintLevel() >=2){
      std::cout<<"****************************************"<<std::endl;
      std::cout<<"Made: "<<fA.size()<<" clusters!"<<std::endl;
      for(unsigned int j=0;j<fA.size();j++){
        std::cout<<"Cluster "<<j<<": "<<fA[j]*180.0/M_PI<<std::endl;
      }
      for (unsigned int j=0; j<fX.size(); ++j) {
        for (unsigned int k=0; k<fA.size(); ++k) {
          std::cout << "U[" << k << "," << j << "]=" 
          << fU[k][j] << " | ";
        }
        std::cout<<std::endl;
      }
    }//end debugging

  }

  //......................................................................
  
  /// Function to make n clusters.
  /// Function loops until optimum number of cluster centers are found
  void FuzzyKMeanAlg::MakeClusters(unsigned int n) 
  {
    unsigned int i;
  
    fD.resize(n); //resize vector holding distance to cluster centers
    fU.resize(n); //resize vector holding membership probabilities
    for (i=0; i<n; ++i) {
      fD[i].resize(fX.size());
      fU[i].resize(fX.size());
    }

    if (fParams.PrintLevel() >= 2){
      std::cout<<"Doing Clustering for "<<n<<" seeds: "<<std::endl;
      for(unsigned int i=0; i<n; i++){
        std::cout<<"seed["<<i<<"]: "<<fA[i]*180.0/M_PI<<std::endl;
      }
    }
  
    this->ComputeDij();//compute distances to cluster centers
    this->ComputeUij();//compute membership probabilities

    //iterate cluster optimization until cluster centers stop changing
    for (i=0; i<fParams.MaxItr(); ++i) {
      bool done = this->ComputeA(); //calculate new clusters centers
      this->ComputeDij(); //compute distances
      this->ComputeUij(); //compute membership
      if (done) break;
    }
  }

  //......................................................................
  
  /// Function to turn the hits in the hitlist from points into angles 
  /// (in radians) with respect to a vertex candidate.  This is done for
  /// a specified view.
  void FuzzyKMeanAlg::MakeAngles(double vXY, double vZ, geo::View_t view)
  {
    fX.clear();
    fdX.clear();
    fW.clear();
    fDist.clear();
    
    //loop over hitlist calculating angles
    for(unsigned int i=0;i<fHL.size();i++){
      //make sure hit is in the right view
      if (fHL[i].fView == view){
        //calculate angle
        fX.push_back(atan2(fHL[i].fXY-vXY,fHL[i].fZ-vZ));
        //add weight (cell energy)
        fW.push_back(fHL[i].fW);
        //distance from vertex
        double dist = util::pythag(fHL[i].fXY-vXY,fHL[i].fZ-vZ);
        //if cell center is the same as the vertex, alter distance slightly to not divide by 0
        if (dist == 0.0) dist = 0.0001;
        fDist.push_back(dist);

        //uncertainty in angle, fit came from muon simulation
        //could be adjusted for different degrees of fuzziness
        //fit is of the form dx = [p1]/x + [p2] + [p3]*x
        if (dist<fParams.Cutoff()) fdX.push_back(fParams.Par1()/dist + fParams.Par2() + fParams.Par3()*dist);
        //beyond the cutoff distance keep uncertainty fixed
        else fdX.push_back(fParams.Par1()/dist + fParams.Par2() + fParams.Par3()*fParams.Cutoff());
      }//end if
    }//end for loop

  }
  
  //......................................................................

  /// Function to produce n seeds for the candidate clusters
  void FuzzyKMeanAlg::SeedA(unsigned int n) 
  {
    //produce intial seeds evenly distributed about circle
    //this is done in case there are not enough dense regions
    //to provide the desired number of seeds
    double step = M_PI/n;
    double angle = 0.0;
    for (unsigned int i=0; i<n; ++i) {
      fA[i]  = angle;
      angle += step;
      if (angle > M_PI) angle -= 2*M_PI;
    }
    
    //calculate the angular density in 360 bins from -M_PI to M_PI
    std::vector<double> w(fParams.Bins());
    double diff;
    step = (double)(2.0*M_PI)/fParams.Bins();
    for (unsigned int k=0; k<w.size(); k++){
      w[k] = 0.0;
      for (unsigned int j=0; j<fX.size(); j++){
        diff = -M_PI + k*step - fX[j];
        while (diff <-M_PI) diff += 2*M_PI;
        while (diff > M_PI) diff -= 2*M_PI;
        w[k] += exp(-diff*diff/(fdX[j]*fdX[j]));
      }
    }
    
    //now find local maxima in weight distribution
    fLMax.clear();
    for (unsigned int k=0; k<w.size(); k++) {
      if (k==0){
        if ((w[k] >= w.back()) && (w[k] > w[k+1])) {
          fLMax.push_back(std::make_pair(w[k],-M_PI+k*step));
        }
      }
      else if (k==w.size()-1) {
        if ((w[k] >= w[k-1]) && (w[k] > w.front())) {
          fLMax.push_back(std::make_pair(w[k],-M_PI+k*step));
        }
      }
      else {
        if ((w[k] >= w[k-1]) && (w[k] > w[k+1])) {
          fLMax.push_back(std::make_pair(w[k],-M_PI+k*step));
        }
      }
    }
    // get sorted density list
    std::sort(fLMax.begin(),fLMax.end());
    std::reverse(fLMax.begin(),fLMax.end());
    
    //populate cluster seeds starting with the maximum density
    for (unsigned int k=0; k<n; k++){
      if (k<fLMax.size()) fA[k] = fLMax[k].second;
      fSeed.push_back(fA[k]);      
    }
        
  }
  
  //......................................................................
  
  /// Function to calculate distance from an angle to cluster center
  double FuzzyKMeanAlg::Distance(double mu, double x, double dx) 
  {
    double d = x-mu;
    while (d<-M_PI) d += 2*M_PI;
    while (d> M_PI) d -= 2*M_PI;
    return d/dx; 
  }
  
  //......................................................................
  
  /// Function to calculate distance metric for all i clusters and j angles, with uncertainty
  void FuzzyKMeanAlg::ComputeDij() 
  {
    unsigned int i, j;
    double d;
    for (i=0; i<fA.size(); ++i) {
      for (j=0; j<fX.size(); ++j) {
        d = this->Distance(fA[i], fX[j], fdX[j]);
        fD[i][j] = d*d;
        if (fParams.PrintLevel() >= 2){
          std::cout << "D[" << i << "," << j << "]=(" 
                    << fX[j] << "-" << fA[i] << ")/" 
                    << fdX[j] << "=" << fD[i][j]
                    << std::endl;
        }
      }//end loop over j
    }//end loop over it
  }
  
  //......................................................................
  
  /// Function to compute membership of cells to clusters
  void FuzzyKMeanAlg::ComputeUij() 
  {
    unsigned int i, j;
    for (j=0; j<fX.size(); ++j) {
      for (i=0; i<fA.size(); ++i) {
        fU[i][j] = exp(-(fParams.Fuzz()*sqrt((double)fA.size())*fD[i][j])/fParams.Beta());
        if (fU[i][j] != fU[i][j]) fU[i][j] = 1e-50; //clean up NaN problem
        if (fU[i][j] < 1e-50) fU[i][j] = 1e-50; //set a membership floor
        if (fParams.PrintLevel() >= 2){
          std::cout << "U[" << i << "," << j << "]=" 
                    << fU[i][j] << " | ";
        }
      } //end loop over i
      if (fParams.PrintLevel() >= 2) std::cout<<std::endl;
    } //end loop over j
  }
  
  //......................................................................
  
  /// Function to update cluster centers.
  /// Returns true if all the centers have stabalized within tolerance,
  /// false if another iteration is needed.
  bool FuzzyKMeanAlg::ComputeA() 
  {
    unsigned int i, j;
    double s;
    double w;
    double musave;
    bool   done = true;
    for (i=0; i<fA.size(); ++i) {
      musave = fA[i];
      s = 0.0;
      w = 0.0;

      for (j=0; j<fX.size(); ++j) {
        double dtheta = fX[j]-musave;
        while (dtheta<-M_PI) dtheta += 2*M_PI;
        while (dtheta> M_PI) dtheta -= 2*M_PI;
        s += pow(fU[i][j],fParams.Fuzz())*dtheta/(fdX[j]*fdX[j]);
        w += pow(fU[i][j],fParams.Fuzz())/(fdX[j]*fdX[j]);
      }

      double thetaNew = musave + s/w;
      while (thetaNew<-M_PI) thetaNew += 2*M_PI;
      while (thetaNew> M_PI) thetaNew -= 2*M_PI;
      fA[i] = thetaNew;
      if (fabs(fA[i]-musave) > fParams.Tol()) done = false;
      if (fParams.PrintLevel() >= 2) std::cout << "COMP MU["<<i<<"]="<<fA[i]<<std::endl;
    }
    return done;
  }
  
  //......................................................................
  
  /// Function to calculate metric for clustering quality
  double FuzzyKMeanAlg::Jpfuz() 
  {
    unsigned int i, j;
    double sum = 0.0;
    for (i=0; i<fA.size(); ++i) {
      for (j=0; j<fX.size(); ++j) {
        sum += pow(fU[i][j],fParams.Fuzz());
      }
    }
    sum = -sum*fParams.Beta()/(fParams.Fuzz()*fParams.Fuzz()*sqrt((double)fA.size()));
    return sum;
  }
  
  //......................................................................
  
  /// Function to find cells left out of clusters, use them to pick new cluster seed
  /// Return true if a new seed candidate is found and clustering should continue,
  /// return false if all hits are in clusters, or there are no new seeds to try.
  bool FuzzyKMeanAlg::FindUnclustered()
  {
    bool fContinue = false; //flag to continue clustering
    fUnused.clear();
    bool unwanted;
    int count=0;
    const double step = (double)(2.0*M_PI)/fParams.Bins();

    //loop over hits check to see if there are leftovers
    for (unsigned int j=0; j<fX.size(); j++){
      unwanted = true;
      for (unsigned int i=0; i<fA.size(); i++){
        //minimum membership set in FHiCL, break if a hit has membership somewhere
        if (fU[i][j] >= fParams.Memb()){ 
          unwanted = false; 
          break;
        }
      } //end loop over i
      if (unwanted){
        fUnused.push_back(std::make_pair(fX[j],fdX[j]));
        count++;
      }
    } //end loop over j
    
    //if there are unclustered hits, find the most appropriate seed to add
    if (count > 0){
      std::vector<double> w(fParams.Bins());
      double diff;
      for (unsigned int k=0; k<w.size(); k++){
        w[k] = 0.0;
        for (unsigned int j=0; j<fUnused.size(); j++){
          diff = -M_PI + k*step - fUnused[j].first;
          while (diff < -M_PI) diff += 2*M_PI;
          while (diff >  M_PI) diff -= 2*M_PI;
          w[k] += exp(-diff*diff/(fUnused[j].second*fUnused[j].second));
        }
      }
    
      //now find local maxima in weight distribution
      fLMax.clear();
      for (unsigned int k=0; k<w.size(); k++) {
        if (k==0){
          if ((w[k] >= w.back()) && (w[k] > w[k+1])) {
            fLMax.push_back(std::make_pair(w[k],-M_PI+k*step));
          }
        }
        else if (k==w.size()-1) {
          if ((w[k] >= w[k-1]) && (w[k] > w.front())) {
            fLMax.push_back(std::make_pair(w[k],-M_PI+k*step));
          }
        }
        else {
          if ((w[k] >= w[k-1]) && (w[k] > w[k+1])) {
            fLMax.push_back(std::make_pair(w[k],-M_PI+k*step));
          }
        }
      }
      std::sort(fLMax.begin(),fLMax.end());
      std::reverse(fLMax.begin(),fLMax.end());
      
      //loop over peaks in the unclustered distribution,
      //try to add highest peak to the list of seed candidates
      //In some cases an unclustered hit is at an angle where even
      //if it has a seed the cluster will converge on an already found
      //cluster center.  So check that the seed to be added has not already
      //been used
      for (unsigned int i=0; i<fLMax.size(); i++){
        bool addSeed = true;
        //loop over seeds to see if it exists
        for (unsigned int j=0; j<fSeed.size(); j++){
          if (fSeed[j] == fLMax[i].second){ 
            addSeed = false; 
            break;
          }
        }

        if (addSeed){
          fA.push_back(fLMax[i].second);
          fSeed.push_back(fLMax[i].second);
          fContinue = true;
          if (fParams.PrintLevel() >= 2){
            std::cout<<"ADDING NEW CLUSTER AT: "
                     <<fLMax[i].second*180.0/M_PI<<std::endl;
          }
          break;
        }

      }//end loop over i

    }//end if statement
    return fContinue;
  }
  
  //......................................................................
  
  /// Function to check the final cluster centers and see if any are duplicates
  /// If a duplicate is found a replacement seed will be substituted or the duplicate
  /// seed will be removed if no replacement is available.  Returns true if the clustering
  /// algorithm needs to run again with the updated seeds.  Returns false if all clusters found
  /// are unique.
  bool FuzzyKMeanAlg::CheckUnique()
  {
    
    bool recluster = false;
    for (unsigned int i=0; i<fA.size(); i++){
      for (unsigned int j=i+1; j<fA.size(); j++){
        //check cluster centers for duplicates
        if (fabs(fA[i]-fA[j]) < fParams.Duplicate()){
          recluster = true;

          //if there is a duplicate see if a new seed on the list can be found
          bool replace = true;
          for (unsigned int k=0; k<fLMax.size(); k++){
            replace = true;
            for(unsigned int m=0; m<fSeed.size(); m++){
              if (fSeed[m] == fLMax[k].second){ 
                replace = false; 
                break;
              }
            }
            if (replace){
              fA[j] = fLMax[k].second;
              fSeed.push_back(fLMax[k].second);
              break;
            }
          }

          //if a replacement seed could not be found, remove the duplicate cluster
          if (!replace){
            fA.erase(fA.begin()+j);
            j--;
          }

        }//end if statement
      }// end loop over j
    }//end loop over i
    return recluster;
  }

  //......................................................................
  
  /// Function to check clusters and determine if they share enough hits to be merged into one cluster
  void FuzzyKMeanAlg::MergeClusters()
  {
    double diff;
    unsigned int i=0;
    unsigned int size = fA.size();

    while(i<size){
      //loop over clusters and search for pairs within FHiCL specified cutoff (default of 30 degrees)
      for (unsigned int j=i+1; j<fA.size(); j++){
        diff = fabs(fA[i] - fA[j]);
        while (diff >  M_PI) diff -= 2*M_PI;

        if (diff < fParams.MergeAng()){
          //now check out overlap between two clusters
          std::vector<int> iClustCells; //cells with membership in i cluster
          std::vector<int> jClustCells; //cells with membership in j cluster
          int overlap = 0; //number of cells appearing in both clusters
          for (unsigned int k=0; k<fX.size(); k++){
            if (fU[i][k] >= fParams.Memb()) iClustCells.push_back(k);
            if (fU[j][k] >= fParams.Memb()) jClustCells.push_back(k);
            if ((fU[i][k] >= fParams.Memb())&&(fU[j][k] >= fParams.Memb())) ++overlap;
          }

          //now computer fraction of smaller cluster that belongs to bigger cluster
          double frac = 0.0;
          if (iClustCells.size() <= jClustCells.size()){
            frac = (double)overlap/iClustCells.size();

            //if small cluster is above minimum membership in large cluster, merge
            if (frac >= fParams.MergeMemb()){
              //add cells to larger cluster with minimum membership
              for (unsigned int k=0; k<fX.size(); k++){
                if ((fU[i][k] >= fParams.Memb())&&(fU[j][k] < fParams.Memb())) fU[j][k] = fParams.Memb();
              }
              //recompute mean of larger cluster
              double s = 0.0;
              double w = 0.0;
              for (unsigned int k=0; k<fX.size(); ++k) {
                s += pow(fU[j][k],fParams.Fuzz())*fX[k]/(fdX[k]*fdX[k]);
                w += pow(fU[j][k],fParams.Fuzz())/(fdX[k]*fdX[k]);
              }
              fA[j] = s/w;
              //now remove the small cluster from the list
              fA.erase(fA.begin()+i);
              fU.erase(fU.begin()+i);
              fD.erase(fD.begin()+i);
              break;
            }

          }
          // else if the j cluster is the smaller one
          else{ 
            frac = (double)overlap/jClustCells.size();

            //if small cluster is above minimum membership in large cluster, merge
            if (frac >= fParams.MergeMemb()){
              //add cells to larger cluster with minimum membership
              for (unsigned int k=0; k<fX.size(); k++){
                if ((fU[j][k] >= fParams.Memb())&&(fU[i][k] < fParams.Memb())) fU[i][k] = fParams.Memb();
              }
              //recompute mean of larger cluster
              double s = 0.0;
              double w = 0.0;
              for (unsigned int k=0; k<fX.size(); ++k) {
                s += pow(fU[i][k],fParams.Fuzz())*fX[k]/(fdX[k]*fdX[k]);
                w += pow(fU[i][k],fParams.Fuzz())/(fdX[k]*fdX[k]);
              }
              fA[i] = s/w;
              //now remove the small cluster from the list
              fA.erase(fA.begin()+j);
              fU.erase(fU.begin()+j);
              fD.erase(fD.begin()+j);
              continue;
            }

          } //end else
        }//end if clusters are close in angle
      }//end loop on j

      //if a merge was made, continue to see if other merges are possible with that cluster
      if (size != fA.size()){
        size = fA.size();
        continue;
      }

      ++i; //if no merge has been made increase i iterator to next cluster
    }//end while loop on i

  }
  
  
}//end fuzz namespace

////////////////////////////////////////////////////////////////////////