www.pudn.com > Mean-Shift.rar > fams.cpp
/////////////////////////////////////////////////////////////////////////////
// Name: fams.cpp
// Purpose: fast adaptive meanshift implementation
// Author: Ilan Shimshoni
// Modified by: Bogdan Georgescu
// Created: 08/14/2003
// Version: v0.1
/////////////////////////////////////////////////////////////////////////////
#include
#include
#include
#include
#include
#include
#include "fams.h"
FAMS::FAMS(int no_lsh)
{
hasPoints_ = 0;
nsel_ = 0;
npm_ = 0;
hashCoeffs_ = NULL;
noLSH_=no_lsh;
time_t tt1;
time(&tt1);
srand((unsigned)tt1);
// srand(100);
nnres1_ = 0;
nnres2_ = 0;
}
FAMS::~FAMS()
{
CleanData();
}
void FAMS::CleanData()
{
CleanPoints();
CleanSelected();
CleanPrunedModes();
CleanHash();
}
void FAMS::CleanPoints()
{
if (hasPoints_)
{
delete [] points_;
delete [] data_;
delete [] rr_;
hasPoints_ = 0;
}
}
void FAMS::CleanSelected()
{
if (nsel_ > 0)
{
delete [] psel_;
delete [] modes_;
delete [] hmodes_;
nsel_ = 0;
}
}
void FAMS::CleanPrunedModes()
{
if (npm_ > 0)
{
delete [] prunedmodes_;
delete [] nprunedmodes_;
npm_ = 0;
}
}
void FAMS::CleanHash()
{
if (hashCoeffs_ != NULL)
{
delete [] hashCoeffs_;
hashCoeffs_ = NULL;
}
}
int FAMS::LoadPoints(char* filename)
{
bgLog("Load data points from %s...",filename);
CleanPoints();
FILE* fd;
fd = fopen(filename, "r");
if (!fd)
{
bgLog("Error opening %s\n", filename);
return 1;
}
fscanf(fd, "%d %d", &n_, &d_);
if ((n_<1) || (d_<1))
{
bgLog("Error reading %s\n", filename);
fclose(fd);
return 1;
}
// allocate data
float* pttemp;
pttemp = new float[n_*d_];
int i;
for (i=0; (i<(n_*d_)) && (fscanf(fd, "%g", &pttemp[i]) == 1); i++);
fclose(fd);
if (i!= (n_*d_))
{
bgLog("Error reading %s\n", filename);
delete [] pttemp;
return 1;
}
// allocate and convert to integer
for (i=0, minVal_=pttemp[0], maxVal_=pttemp[0]; i<(n_*d_); i++)
{
if (minVal_>pttemp[i])
minVal_ = pttemp[i];
else if (maxVal_ M)
{
bgLog("LSH hash table overflow exiting\n");
exit(-1);
}
if(hs[where] == Bs)
continue;
HT[where][hs[where]].pt_ = &pt;
HT[where][hs[where]].whichCut_ = which;
HT[where][hs[where]].which2_ = which2;
hs[where]++;
break;
}
}
// compute the pilot h_i's for the data points
void FAMS::ComputePilot(block *HT,int *hs, fams_partition *cuts, char* pilot_fn)
{
const int win_j = 10,max_win=7000;
int i,j;
unsigned int nn;
unsigned int wjd = (unsigned int)(win_j*d_);
int num_l[1000];
fams_res_cont res(n_);
if(LoadBandwidths(pilot_fn)==0)
{
bgLog("compute bandwidths...");
for(j=0; jk_)
{
break;
}
}
points_[j].window_=(nn+1)*wjd;
}
SaveBandwidths(pilot_fn);
}
else
bgLog("load bandwidths...");
for(j=0; jk_)
{
break;
}
}
points_[who].window_=(nn+1)*win_j;
}
} else
{
for(j=0; jk_)
break;
}
for(; (num_l[nl]-1) == i; nl++)
scores[nl] += (float) (((nn+1.0)*win_j)/points_[who].window_);
}
}
}
else
{
GetNearestNeighbours(points_[who],HT,hs,cuts,res,0,num_l);
nel = res.nel_;
num_l[L_] = n_+1;
for(i=0; ik_)
break;
}
for(; (num_l[nl]-1) == i; nl++)
scores[nl] += (float) (((nn+1.0)*win_j)/points_[who].window_);
}
}
}
numn=0;
for(nn=0; nnk_)
break;
}
}
for(j=0; j M2)
{
bgLog(" Hash Table2 full\n");
exit(-1);
}
for(uu=0; uudata_, dataSize_);
crtH = &hmodes_[jj];
*crtH = currentpt->window_;
tMode[jj]=(unsigned short*)1;
int iter;
for(iter=0; NotEq(oldMean,crtMean) && (iter=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> FAMS_PRUNE_HDIV;
if (cminDist < hprune)
{
// aready in, just add
for (cd=0; cd2000)
{
for (i=0; i=0; i--)
{
if (stemp[i]>=npmin)
nrel++;
else
break;
}
if (nrel > FAMS_PRUNE_MAXM)
nrel = FAMS_PRUNE_MAXM;
// rearange only relevant modes
mcount2 = new int[nrel];
cmodes2 = new float[d_*nrel];
for (i=0; i> FAMS_PRUNE_HDIV;
if (cminDist < hprune)
{
// aready in, just add
for (cd=0; cd=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; i0)
{
adaptive = 0;
hWidth = (int) (65535.0*(width)/(maxVal_-minVal_));
}
k_=k;
epsilon += 1;
// select points on which test is run
SelectMsPoints(FAMS_FKL_NEL*100.0/n_, 0);
// compute bandwidths for selected points
ComputeRealBandwidths(hWidth);
// start finding the correct l for each k
float scores[FAMS_FKL_TIMES*FAMS_MAX_L];
int Lcrt, Kcrt;
int nBest;
int LBest[FAMS_MAX_K];
int KBest[FAMS_MAX_K];
int ntimes, is;
Lcrt = Lmax;
bgLog(" find valid pairs");
for (Kcrt = Kmax, nBest=0; Kcrt >= Kmin; Kcrt -= Kjump, nBest++)
{
// do iterations for crt K and L = 1...Lcrt
for (ntimes=0; ntimes0)
Lcrt = LBest[nBest]+2;
}
bgLog("done\n");
//start finding the pair with best running time
double run_times[FAMS_FKL_TIMES];
int iBest, i;
double timeBest=-1;
bgLog(" select best pair\n");
for (i=0; irun_times[FAMS_FKL_TIMES/2]))
{
iBest = i;
timeBest = run_times[FAMS_FKL_TIMES/2];
}
bgLog(" K=%d L=%d time: %g\n", KBest[i], LBest[i], run_times[FAMS_FKL_TIMES/2]);
}
K = KBest[iBest];
L = LBest[iBest];
bgLog("done\n");
return 0;
}
double FAMS::DoFindKLIteration(int K,int L, float* scores)
{
K_ = K;
L_ = L;
int i, j;
// Allocate memory for the hash table
M_ = GetPrime(3*n_*L_/(Bs));
block *HT = new block[M_];
int *hs = new int[M_];
InitHash(K_+L_);
memset(hs,0,sizeof(int)*M_);
// Build partitions
fams_partition *cuts = new fams_partition[L_];
for(i=0; i<20; i++)
rand();
int cut_res[FAMS_MAX_K];
MakeCuts(cuts);
//Insert data into partitions
for(j=0; j0)
{
adaptive = 0;
hWidth = (int) (65535.0*(width)/(maxVal_-minVal_));
}
K_=K; L_=L; k_=k;
SelectMsPoints(percent, jump);
// Allocate memory for the hash table
M_ = GetPrime(3*n_*L_/(Bs));
M2_ = GetPrime((int)(nsel_*20*3/Bs2));
block *HT = new block[M_];
int *hs = new int[M_];
block2 *HT2 = new block2[M2_];
int * hs2 = new int[M2_];
memset(hs,0,sizeof(int)*M_);
// Build partitions
fams_partition *cuts = new fams_partition[L];
for(i=0; i<20; i++)
rand();
int cut_res[FAMS_MAX_K];
MakeCuts(cuts);
InitHash(K_+L_);
//Insert data into partitions
for(j=0; j 6)
{
for (i=6; i1)) percent = 0;
break;
case 'h':
i++;
width = (float) atof(argv[i]);
break;
case 'f':
i++;
epsilon = (float) atof(argv[i++]);
Kmin = atoi(argv[i++]);
Kjump = atoi(argv[i]);
Lmax = L; Kmax = K;
find_kl=1;
break;
default:
bgLog("Error in param %s\n", argv[i]);
usage();
exit(1);
break;
}
}
}
// load points
FAMS cfams(no_lsh);
if (cfams.LoadPoints(fdata_file_name))
return 1;
// find K L (if necessary)
if (find_kl)
{
cfams.FindKL(Kmin, Kmax, Kjump, Lmax, k_neigh, width, epsilon, K, L);
bgLog("Found K = %d L = %d (write them down)\n", K, L);
int ch=' ';
do{
bgLog("Do you want to run FAMS with this (K=%d,L=%d) pair? (y/n)",K,L);
ch = getchar();
if ((ch == 'n') || (ch == 'N'))
return 0;
} while((ch != 'y') && (ch != 'Y'));
}
sprintf(fdata_file_name, "%spilot_%d_%s.txt", input_path, k_neigh, data_file_name);
cfams.RunFAMS(K, L, k_neigh, percent, jump, width, fdata_file_name);
// save the data
sprintf(fdata_file_name,"%sout_%s.txt",input_path,data_file_name);
cfams.SaveModes(fdata_file_name);
// save pruned modes modes
sprintf(fdata_file_name,"%smodes_%s.txt",input_path, data_file_name);
cfams.SavePrunedModes(fdata_file_name);
return 0;
}