fams.cpp

mean-shift算法的实现算例

      if(hs2[where] == Bs2)
         continue;
      HT2[where][hs2[where]].dp_ = solution;
      HT2[where][hs2[where]++].whichCut_ = which;
      break;
   }
}

//perform an LSH query

void FAMS::GetNearestNeighbours(fams_point& who,block *HT, int *hs,
     fams_partition *cuts, fams_res_cont& res, int print,int num_l[])
{
   int i;
   for(i=0; i<L_; i++)
      EvalCutRes(who,cuts[i],t_cut_res_[i]);
   if(CompareCutRes(t_cut_res_,t_old_cut_res_)==0)
   {
      return;
   }
   memcpy(t_old_cut_res_,t_cut_res_,sizeof(t_old_cut_res_));
   res.clear();
   nnres1_++;
   for(i=0; i<L_; i++)
   {
      int hjump;
      int m = HashFunction(t_cut_res_[i],i,K_,M_,&hjump);
      int m2 = HashFunction(&t_cut_res_[i][1],i,K_-1);
      AddDataToRes(HT,hs,res,m,Bs,M_,i,nnres1_,m2,hjump);
      num_l[i] = res.nel_;
   }
}

// perform an LSH query using in addition the second hash table

unsigned short* FAMS::GetNearestNeighbours2H(unsigned short *who, block*HT, int *hs,
       fams_partition* cuts, fams_res_cont& res,
       unsigned short **solution, block2 *HT2, int *hs2)
{
   int i;
   for(i=0; i<L_; i++)
   {
      EvalCutRes(who,cuts[i],t_cut_res_[i]);
      t_m_[i] = HashFunction(t_cut_res_[i],i,K_,M_,&t_hjump_[i]);
      t_m2_[i] = HashFunction(&t_cut_res_[i][1],i,K_-1);
   }
   if(FAMS_DO_SPEEDUP)
   {
      int hjump2;
      int hf = HashFunction(t_m_,0,L_,M2_,&hjump2);
      int hf2 = HashFunction(&t_m_[L_/2-1],0,L_/2);
      unsigned short *old_sol =  FindInHash(HT2,hs2,hf,hf2,M2_,hjump2);
      if(old_sol!= NULL && old_sol != (unsigned short*)1)
         return old_sol;
      
      if(old_sol == NULL)
         InsertIntoHash(HT2,hs2,hf,hf2,solution,M2_,hjump2);
   }
   if(memcmp(t_m_,t_old_m_,sizeof(int)*L_)==0)
   {
      return NULL;
   }
   memcpy(t_old_m_,t_m_,sizeof(int)*L_);
   res.clear();
   nnres2_++;
   for(i=0; i<L_; i++)
      AddDataToRes(HT,hs,res,t_m_[i],Bs,M_,i,nnres2_,t_m2_[i],t_hjump_[i]);

   return NULL;
}

inline double SQ(double x) {return x*x;}

// perform an FAMS iteration

unsigned int FAMS::DoMeanShiftAdaptiveIteration(fams_res_cont& res,
              unsigned short *old, unsigned short *ret)
{
   double total_weight =0;
   int i,j;
   double dist;
   for(i=0; i<d_; i++) rr_[i]=0;
   int nel;
   if(noLSH_)
      nel = n_;
   else{
      nel = res.nel_;
   }
   fams_pointp ptp;
   unsigned int crtH;
   double hmdist=1e100;
   for(i=0; i<nel; i++)
   {
      ptp = noLSH_ ? &points_[i] : res.vec_[i];
      if(DistL1Data(old, *ptp, (*ptp).window_, dist))
      {
         double w;
         w = (*(ptp)).weightdp2_*SQ(1.0-(dist/(*ptp).window_));
         total_weight+=w;
         for(j=0; j<d_; j++)
            rr_[j] += (*(ptp)).data_[j]*w;
         if (dist<hmdist)
         {
            hmdist = dist;
            crtH=(*ptp).window_;
         }
      }
   }	
   if(total_weight==0)
   {
      return 0;
   }
   for(i=0; i<d_; i++)
      ret[i] = (unsigned short)(rr_[i]/total_weight);
   return crtH;
}

// perform a query on the second hash table

unsigned short* FAMS::FindInHash(block2* HT2,int *hs2,int where,int which,int M2,int hjump)
{
   int uu;
   int nw=0;
   for(;;where = (where+hjump)%M2)
   {
      nw++;
      if(nw > M2)
      {
         bgLog(" Hash Table2 full\n");
         exit(-1);
      }
      for(uu=0; uu<hs2[where]; uu++)
         if(HT2[where][uu].whichCut_ == which)
            return *(HT2[where][uu].dp_);
         if(hs2[where] < Bs2)
            break;
   }
   return NULL;
}


// Perform a query to one partition and retreive all the points in the cell 

void FAMS::AddDataToRes(block HT[],int hs[],fams_res_cont& res,int where,
        int Bs,int M,int which, unsigned short nnres, int which2, int hjump)
{
   int uu;
   for(;;where = (where+hjump)%M)
   {
      for(uu=0; uu<hs[where]; uu++)
      {
         if(HT[where][uu].whichCut_ == which && HT[where][uu].which2_ == which2)
         {
            if((*HT[where][uu].pt_).usedFlag_!=nnres)
            {
               res.push_back(HT[where][uu].pt_);
               (*HT[where][uu].pt_).usedFlag_=nnres;
            }
         }
      }
      if(hs[where] < Bs)
         break;
   }
}

// perform FAMS starting from a subset of the data points.
void FAMS::DoFAMS(block *HT,int *hs, fams_partition *cuts,
     block2* HT2, int *hs2)
{
   int jj;
   fams_res_cont res(n_);
   memset(HT2,0,sizeof(block2)*M2_);
   memset(hs2,0,sizeof(int)*M2_);
   unsigned short *oldMean;
   unsigned short *crtMean;
   oldMean = new unsigned short[d_];
   crtMean = new unsigned short[d_];
   fams_pointp currentpt;
   unsigned short *sol;
   //unsigned short *crtMode;
   int who;
   unsigned short **tMode;
   unsigned int newH;
   unsigned int *crtH;
   tMode = new unsigned short*[n_];
   bgLog(" Start MS iterations");
   int myPt = nsel_/10;
   for(jj=0; jj<nsel_; jj++)
   {
      if((jj%myPt)==0)
         bgLog(".");
      who = psel_[jj];
      currentpt = &points_[who];
      memcpy(crtMean, currentpt->data_, dataSize_);
      crtH = &hmodes_[jj];
      *crtH = currentpt->window_;
      tMode[jj]=(unsigned short*)1;
      int iter;
      for(iter=0; NotEq(oldMean,crtMean) && (iter<FAMS_MAXITER); iter++)
      {
         if(!noLSH_)
         {
            if((sol=GetNearestNeighbours2H(crtMean,HT,hs,cuts,
               res,&tMode[jj],
               HT2,hs2)))
            {
               if(sol == (unsigned short*)1)
               {
                  tMode[jj] = &modes_[jj*d_];
                  memcpy(tMode[jj], crtMean, dataSize_);
               }
               else
               {
                  tMode[jj] = &modes_[jj*d_];
                  memcpy(tMode[jj], sol, dataSize_);
                  break;
               }
            }
         }
         memcpy(oldMean, crtMean, dataSize_);
         if(!(newH=DoMeanShiftAdaptiveIteration(res,oldMean,crtMean)))
         {
            memcpy(crtMean, oldMean, dataSize_);
            break;
         }
         *crtH=newH;
      }
      if(tMode[jj]==(unsigned short*)1)
      {
         tMode[jj] = &modes_[jj*d_];
         memcpy(tMode[jj], crtMean, dataSize_);
      }
   }
   delete [] oldMean;
   delete [] crtMean;
   delete [] tMode;
   bgLog("done.\n");
}

void FAMS::SaveModes(char* fn)
{
   if (nsel_ < 1)
      return;
   bgLog("Save convergence points in %s ...", fn);
   FILE* fd;
   fd = fopen(fn, "wb");
   int i,j,idx;
   idx = 0;
   for (i=0; i<nsel_; i++)
   {
      for (j=0; j<d_; j++)
      {
         fprintf(fd, "%g ", modes_[idx++]*(maxVal_-minVal_)/65535.0 + minVal_);
      }
      fprintf(fd, "\n");
   }
   fclose(fd);
   bgLog("done\n");
}

void FAMS::SavePrunedModes(char* fn)
{
   if (npm_ < 1)
      return;
   bgLog("Save joined convergence points in %s ...", fn);
   FILE* fd;
   fd = fopen(fn, "wb");
   int i,j,idx;
   idx = 0;
   for (i=0; i<npm_; i++)
   {
      fprintf(fd, "%d  ", nprunedmodes_[i]);
      for (j=0; j<d_; j++)
      {
         fprintf(fd, "%g ", prunedmodes_[idx++]*(maxVal_-minVal_)/65535.0 + minVal_);
      }
      fprintf(fd, "\n");
   }
   fclose(fd);
   bgLog("done\n");
}
/*
int FAMS::PruneModes(int hprune, int npmin)
{
   bgLog(" Join Modes with h=%g and min pt=%d...", hprune*(maxVal_-minVal_)/65535.0, npmin);
   if (nsel_ < 1)
      return 1;
   hprune *= d_;

   int *mcount;
   float *cmodes, *ctmodes;
   unsigned short *pmodes;
   double cminDist, cdist;
   int iminDist, cref;
   mcount = new int[nsel_];
   cmodes = new float[d_*nsel_];

   int i, cd, cm, maxm;

   memset(mcount, 0, nsel_*sizeof(int));

   // set first mode
   for (cd = 0; cd<d_; cd++)
   {
      cmodes[cd] = modes_[cd];
   }
   mcount[0] = 1;
   maxm = 1;

   for (cm = 1; cm<nsel_; cm++)
   {
      bgLog("cm=%d, nm=%d\n",cm, maxm);

      pmodes = modes_+cm*d_;

      // compute closest mode
      cminDist = d_*1e7;
      iminDist = -1;
      for (cref = 0; cref<maxm; cref++)
      {
         cdist = 0;
         ctmodes = cmodes+cref*d_;
         for (cd=0; cd<d_; cd++)
            cdist += fabs(ctmodes[cd]/mcount[cref] - pmodes[cd]);
         if (cdist<cminDist)
         {
            cminDist = cdist;
            iminDist = cref;
         }
      }
      // join
      if (cminDist < hprune)
      {
         // aready in, just add
         for (cd=0; cd<d_; cd++)
         {
            cmodes[iminDist*d_+cd] += pmodes[cd];
         }
         mcount[iminDist] += 1;
      } else
      {
         // new mode, create
         for (cd=0; cd<d_; cd++)
         {
            cmodes[maxm*d_+cd] = pmodes[cd];
         }
         mcount[maxm] = 1;
         maxm += 1;
      }
   }

   // put the modes in the order of importance (count)
   int* stemp;
   int* istemp;
   stemp = new int[maxm];
   istemp = new int[maxm];
   for (i=0; i<maxm; i++)
   {
      stemp[i] = mcount[i];
      istemp[i] = i;
   }
   bgISort(stemp, maxm, istemp); // increasing

   // find number of relevant modes
   int nrel=1;
   for (i=maxm-2; i>=0; i--)
   {
      if (stemp[i]>=npmin)
         nrel++;
      else
         break;
   }

   CleanPrunedModes();
   prunedmodes_ = new unsigned short[d_*nrel];
   nprunedmodes_ = new int[nrel];
   unsigned short* cpm;
   npm_ = nrel;

   cpm = prunedmodes_;
   for (i=0; i<npm_; i++)
   {
      nprunedmodes_[i] = stemp[maxm-i-1];
      cm = istemp[maxm-i-1];
      for (cd=0; cd<d_; cd++)
      {
         *(cpm++) = (unsigned short) (cmodes[cm*d_+cd]/mcount[cm]);
      }
   }


   delete [] istemp;
   delete [] stemp;


   delete [] cmodes;
   delete [] mcount;

   bgLog("done\n");
   return 1;
}
*/

int FAMS::PruneModes(int hprune, int npmin)
{
   // compute jump
   int jm = ceil(((double) nsel_)/FAMS_PRUNE_MAXP);

   bgLog(" Join Modes with adaptive h/%d, min pt=%d, jump=%d\n", (int) pow(2,FAMS_PRUNE_HDIV), npmin, jm);
   bgLog("            pass 1");
   if (nsel_ < 1)
      return 1;
   hprune *= d_;

   int *mcount, *mcount2;
   float *cmodes, *ctmodes, *cmodes2;
   unsigned short *pmodes;
   double cminDist, cdist;
   int iminDist, cref;
   unsigned char *invalidm;
   invalidm = new unsigned char[nsel_];
   mcount = new int[nsel_];
   cmodes = new float[d_*nsel_];

   int i, cd, cm, maxm;

   memset(mcount, 0, nsel_*sizeof(int));
   memset(invalidm, 0, nsel_*sizeof(unsigned char));

   // set first mode
   for (cd = 0; cd<d_; cd++)
   {
      cmodes[cd] = modes_[cd];
   }
   mcount[0] = 1;
   maxm = 1;

   int myPt = FAMS_PRUNE_MAXP/10;
   for (cm = 1; cm<nsel_; cm+=jm)
   {
      if((cm%myPt)==0)
         bgLog(".");

      pmodes = modes_+cm*d_;

      //bgLog("cm=%d, nm=%d, kk=%d, %d %d\n",cm, maxm, hmodes_[cm], pmodes[0], pmodes[d_-1]);

      // compute closest mode
      cminDist = d_*1e7;
      iminDist = -1;
      for (cref = 0; cref<maxm; cref++)
      {
         if (invalidm[cref])
            continue;
         cdist = 0;
         ctmodes = cmodes+cref*d_;
         for (cd=0; cd<d_; cd++)
            cdist += fabs(ctmodes[cd]/mcount[cref] - pmodes[cd]);
         if (cdist<cminDist)
         {
            cminDist = cdist;
            iminDist = cref;
         }
      }
      // join
      hprune = hmodes_[cm] >> FAMS_PRUNE_HDIV;
      if (cminDist < hprune)
      {
         // aready in, just add
         for (cd=0; cd<d_; cd++)
         {
            cmodes[iminDist*d_+cd] += pmodes[cd];
         }
         mcount[iminDist] += 1;
      } else
      {
         // new mode, create
         for (cd=0; cd<d_; cd++)
         {
            cmodes[maxm*d_+cd] = pmodes[cd];
         }
         mcount[maxm] = 1;
         maxm += 1;
      }
      // check for valid modes
      if (maxm>2000)
      {
         for (i=0; i<maxm; i++)
         {
            if (mcount[i] < 3)
               invalidm[i] = 1;
         }
      }
   }

   bgLog("done\n");
   bgLog("            pass 2");
   
   // put the modes in the order of importance (count)
   int* stemp;
   int* istemp;
   stemp = new int[maxm];
   istemp = new int[maxm];
   for (i=0; i<maxm; i++)
   {
      stemp[i] = mcount[i];
      istemp[i] = i;
   }
   bgISort(stemp, maxm, istemp); // increasing

   // find number of relevant modes