Skip to content

Commit

Permalink
Added new flags for binning.
Browse files Browse the repository at this point in the history
  • Loading branch information
@matthewbellis
matthewbellis committed Apr 2, 2012
1 parent 5a342e4 commit 93308ce
Showing 1 changed file with 95 additions and 18 deletions.
113 changes: 95 additions & 18 deletions CUDA_code/Calculate_arc_length_two_datasets.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ using namespace std;
// Kernel to calculate angular distances between galaxies and histogram
// the distances.
////////////////////////////////////////////////////////////////////////
__global__ void distance(float *a0, float *d0, float *a1, float *d1, int xind, int yind, int *dev_hist, int nbins, float* dev_bin_edges, bool two_different_files=1)
__global__ void distance(float *a0, float *d0, float *a1, float *d1, int xind, int yind, int *dev_hist, float hist_min, float hist_max, int nbins, float bin_width, bool log_binning=0, bool two_different_files=1)
{

////////////////////////////////////////////////////////////////////////////
Expand All @@ -31,22 +31,24 @@ __global__ void distance(float *a0, float *d0, float *a1, float *d1, int xind, i
idx += xind;

int i=0;
int j=0;
//int j=0;

////////////////////////////////////////////////////////////////////////////
// Copy over a local version of the bin edges. This will prevent threads
// from accessing the same memory that holds this information.
////////////////////////////////////////////////////////////////////////////
/*
float bin_edges[1024]; // This is hard-coded and will result in bugs if a histogram
// with greater than 1024 bins is used.
// with greater than 1024 bins is used.
int bin_offset = thread_idx*nbins;
for (i=0;i<nbins;i++)
{
bin_edges[i] = dev_bin_edges[i+bin_offset];
}
*/

float hist_min = bin_edges[0];
float hist_max = bin_edges[nbins-1];
//float hist_min = bin_edges[0];
//float hist_max = bin_edges[nbins-1];

float alpha = a0[idx], delta0 = d0[idx];
float cos_d0 = cos(delta0), sin_d0 = sin(delta0), dist;
Expand Down Expand Up @@ -100,6 +102,16 @@ __global__ void distance(float *a0, float *d0, float *a1, float *d1, int xind, i
else if(dist >= hist_max)
bin_index = nbins + 1;
else
if (!log_binning)
{
bin_index = int((dist-hist_min)/bin_width) + 1;
}
else // log binning
{
bin_index = int((log(dist)-log(hist_min))/bin_width) + 1;
}

/*
{
//bin_index = int(((dist - hist_min) * nbins / hist_max) +1);
bin_index = 0;
Expand All @@ -113,6 +125,7 @@ __global__ void distance(float *a0, float *d0, float *a1, float *d1, int xind, i
}
}
}
*/

offset = ((nbins+2)*thread_idx);
bin_index += offset;
Expand Down Expand Up @@ -141,10 +154,33 @@ int main(int argc, char **argv)
char *binning_filename = NULL;
FILE *binning_file = NULL;

while ((c = getopt(argc, argv, "ab:o:")) != -1) {
float hist_lower_range = 0.0000001;
float hist_upper_range = 0;
int nbins = 100;
float hist_bin_width = 0.05;
bool log_binning_flag = 0; // False

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


while ((c = getopt(argc, argv, "ab:o:L:N:lw:")) != -1) {
switch(c) {
case 'a':
printf("a is set\n");
case 'N':
printf("N is set\n");
nbins = atoi(optarg);
break;
case 'L':
printf("L is set\n");
hist_lower_range = atof(optarg);
break;
case 'w':
hist_bin_width = atof(optarg);
printf("Histogram bin width: %f\n",hist_bin_width);
break;
case 'l':
log_binning_flag = 1;
printf("Will use log binning.\n");
break;
case 'b':
binning_filename = optarg;
Expand Down Expand Up @@ -177,6 +213,28 @@ int main(int argc, char **argv)
printf("Output filename is %s\n", outfilename);
}

float temp_lo = hist_lower_range;
if (hist_upper_range == 0)
{
if (log_binning_flag==0)
{
for (int i=0;i<nbins;i++)
{
hist_upper_range = temp_lo + hist_bin_width;
temp_lo = hist_upper_range;
}
}
else
{
for (int i=0;i<nbins;i++)
{
hist_upper_range = exp(log(temp_lo) + hist_bin_width);
temp_lo = hist_upper_range;
}
}
}
printf("hist_upper_range: %f\n",hist_upper_range);

//printf ("%d\n", optind);
//printf ("%d\n", argc);
//printf ("%d\n", optind);
Expand Down Expand Up @@ -299,6 +357,7 @@ int main(int argc, char **argv)
// Read in the file that defines the bin edges.
////////////////////////////////////////////////////////////////////////////

/*
int default_nbins = 27;
float default_bin_edges[DEFAULT_NBINS] = {0.0000,0.001000,0.001585,0.002512,0.003981,0.006310,0.010000,0.015849,0.025119,0.039811,0.063096,0.100000,0.158489,0.251189,0.398107,0.630957,1.000000,1.584893,2.511886,3.981072,6.309573,10.000000,15.848932,25.118864,39.810717,63.095734,100.000000};
Expand Down Expand Up @@ -359,6 +418,7 @@ int main(int argc, char **argv)
}
float hist_min = h_bin_edges[0];
float hist_max = h_bin_edges[nbins-1];
*/

////////////////////////////////////////////////////////////////////////////
// Finished defining the bin edges.
Expand Down Expand Up @@ -407,12 +467,12 @@ int main(int argc, char **argv)
///////////////////////////////////////////////////////////////////////////

int *hist, *dev_hist;
float *dev_bin_edges;
//float *dev_bin_edges;

cudaMalloc((void **) &dev_bin_edges,size_bin_edges_array);
printf("Size of dev_bin_edges: %d bytes\n",size_bin_edges_array);
cudaMemset(dev_bin_edges, 0, nbins*SUBMATRIX_SIZE);
cudaMemcpy(dev_bin_edges, h_bin_edges,size_bin_edges_array, cudaMemcpyHostToDevice );
//cudaMalloc((void **) &dev_bin_edges,size_bin_edges_array);
//printf("Size of dev_bin_edges: %d bytes\n",size_bin_edges_array);
//cudaMemset(dev_bin_edges, 0, nbins*SUBMATRIX_SIZE);
//cudaMemcpy(dev_bin_edges, h_bin_edges,size_bin_edges_array, cudaMemcpyHostToDevice );

int size_hist = SUBMATRIX_SIZE * (nbins+2);
int size_hist_bytes = size_hist*sizeof(int);
Expand Down Expand Up @@ -488,7 +548,8 @@ int main(int argc, char **argv)
// Set the histogram to all zeros each time.
cudaMemset(dev_hist,0,size_hist_bytes);

distance<<<grid,block>>>(d_alpha0, d_delta0,d_alpha1, d_delta1, x, y, dev_hist, nbins, dev_bin_edges, two_different_files);
//__global__ void distance(float *a0, float *d0, float *a1, float *d1, int xind, int yind, int *dev_hist, float hist_min, float hist_max, int nbins, float bin_width, bool log_binning=0, bool two_different_files=1)
distance<<<grid,block>>>(d_alpha0, d_delta0,d_alpha1, d_delta1, x, y, dev_hist, hist_lower_range, hist_upper_range, nbins, hist_bin_width, log_binning_flag, two_different_files);
cudaMemcpy(hist, dev_hist, size_hist_bytes, cudaMemcpyDeviceToHost);

////////////////////////////////////////////////////////////////////
Expand All @@ -503,22 +564,38 @@ int main(int argc, char **argv)
}

unsigned long total = 0;
float bin_width = (hist_max - hist_min) / nbins;
//float bin_width = (hist_upper_range - hist_lower_range) / nbins;
float bins_mid = 0;

fprintf(outfile, "%s %s\n", "Angular Distance(radians)","Number of Entries");
float lo = hist_lower_range;
float hi = 0;
//printf("hist_lower_range: %f\n",hist_lower_range);
//printf("hist_upper_range: %f\n",hist_upper_range);
//printf("hist_bin_width: %f\n",hist_bin_width);
for(int k=0; k<nbins+1; k++)
{
//bins_mid = bin_width*(k - 0.5);

float lo = h_bin_edges[k];
float hi = h_bin_edges[k+1];
//float lo = h_bin_edges[k];
//float hi = h_bin_edges[k+1];
if (!log_binning_flag)
{
hi = lo + hist_bin_width;
}
else
{
//printf("lo: %f\t\tlog(lo): %f\n",lo,log(lo));
hi = exp(log(lo) + hist_bin_width);
}

bins_mid = (hi+lo)/2.0;

fprintf(outfile, "%.3e %s %lu \n", bins_mid, ",", hist_array[k]);
total += hist_array[k];

lo = hi;

}
printf("total: %lu \n", total);

Expand All @@ -537,7 +614,7 @@ int main(int argc, char **argv)
cudaFree(d_alpha1);
cudaFree(d_delta1);
cudaFree(dev_hist);
cudaFree(dev_bin_edges);
//cudaFree(dev_bin_edges);

return 0;
}
Expand Down

0 comments on commit 93308ce

Please sign in to comment.