Skip to content

ddtheta: bin in sep^2 instead of cos(theta) #297

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

ddtheta: bin in sep^2 instead of cos(theta) #297

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 19 additions & 10 deletions mocks/DDtheta_mocks/countpairs_theta_mocks_impl.c.src
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,12 @@
#include "countpairs_theta_mocks_kernels_DOUBLE.c"


DOUBLE chord2_DOUBLE(DOUBLE theta){
// Squared chord length for a given angle on the unit sphere
DOUBLE sinval = SIND(theta/2);
return 4*sinval*sinval;
}


int interrupt_status_wtheta_mocks_DOUBLE=EXIT_SUCCESS;

Expand Down Expand Up @@ -566,12 +572,13 @@ int countpairs_theta_mocks_DOUBLE(const int64_t ND1, DOUBLE *ra1, DOUBLE *dec1,
}


DOUBLE costheta_upp[nthetabin];
DOUBLE sep2_upp[nthetabin];
for(int i=0;i<nthetabin;i++) {
costheta_upp[i] = COSD(theta_upp[i]);
// costheta_upp[i] = COSD(theta_upp[i]);
sep2_upp[i] = chord2_DOUBLE(theta_upp[i]);
}
const DOUBLE costhetamin=costheta_upp[0];
const DOUBLE costhetamax=costheta_upp[nthetabin-1];
const DOUBLE sep2min=sep2_upp[0];
const DOUBLE sep2max=sep2_upp[nthetabin-1];

DOUBLE *X1 = my_malloc(sizeof(*X1),ND1);
DOUBLE *Y1 = my_malloc(sizeof(*Y1),ND1);
Expand Down Expand Up @@ -644,12 +651,13 @@ int countpairs_theta_mocks_DOUBLE(const int64_t ND1, DOUBLE *ra1, DOUBLE *dec1,
}

//this is equivalent to brute force calculating on the entire dataset
assert(1); // TODO
int status = countpairs_theta_mocks_brute_force_DOUBLE(ND1, X1, Y1, Z1,
ND2, X2, Y2, Z2,
numthreads,
costhetamax, costhetamin, nthetabin,
sep2max, sep2min, nthetabin,
theta_upp,
costheta_upp,
sep2_upp,
results,
options,
extra);
Expand Down Expand Up @@ -878,12 +886,13 @@ int countpairs_theta_mocks_DOUBLE(const int64_t ND1, DOUBLE *ra1, DOUBLE *dec1,
}

//this is equivalent to brute force calculating on the entire dataset
assert(1); // TODO
int status = countpairs_theta_mocks_brute_force_DOUBLE(ND1, X1, Y1, Z1,
ND2, X2, Y2, Z2,
numthreads,
costhetamax, costhetamin, nthetabin,
sep2max, sep2min, nthetabin,
theta_upp,
costheta_upp,
sep2_upp,
results,
options,
extra);
Expand Down Expand Up @@ -1019,8 +1028,8 @@ int countpairs_theta_mocks_DOUBLE(const int64_t ND1, DOUBLE *ra1, DOUBLE *dec1,
second->nelements, second->x, second->y, second->z, &(second->weights),
this_cell_pair->same_cell,
options->fast_acos,
costhetamax, costhetamin, nthetabin,
costheta_upp,
sep2max, sep2min, nthetabin,
sep2_upp,
this_cell_pair->min_dx, this_cell_pair->min_dy, this_cell_pair->min_dz,
this_cell_pair->closest_x1, this_cell_pair->closest_y1, this_cell_pair->closest_z1,
this_thetaavg, npairs,
Expand Down
46 changes: 16 additions & 30 deletions mocks/DDtheta_mocks/countpairs_theta_mocks_kernels.c.src
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@ static inline int countpairs_theta_mocks_fallback_DOUBLE(const int64_t N0, DOUBL
const int64_t N1, DOUBLE *x1, DOUBLE *y1,
DOUBLE *z1, const weight_struct_DOUBLE *weights1,
const int same_cell, const int order,
const DOUBLE costhetamax, const DOUBLE costhetamin,
const int nthetabin, const DOUBLE *costheta_upp,
const DOUBLE sep2max, const DOUBLE sep2min,
const int nthetabin, const DOUBLE *sep2_upp,
const DOUBLE min_xdiff, const DOUBLE min_ydiff,
const DOUBLE min_zdiff, const DOUBLE closest_icell_xpos,
const DOUBLE closest_icell_ypos, const DOUBLE closest_icell_zpos,
Expand All @@ -42,16 +42,6 @@ static inline int countpairs_theta_mocks_fallback_DOUBLE(const int64_t N0, DOUBL
return EXIT_FAILURE;
}

/* const DOUBLE max_chord_sep = 2.0*SIND(0.5*thetamax); */
/* C = 2.0 * SIN(thetamax/2)
-> C^2 = 4.0 * SIN^2 (thetamax/2.0)
-> C^2 = 2.0 * (2 * SIN^2(thetamax/2.0))
-> C^2 = 2.0 * (1 - COS(thetamax))

-> COS(theta) = (1 - 0.5*C^2)
*/
const DOUBLE sqr_max_chord_sep = ((DOUBLE) 2.0) * (((DOUBLE) 1.0) - costhetamax);

const int32_t need_rpavg = src_rpavg != NULL;
const int32_t need_weightavg = src_weightavg != NULL;
uint64_t npairs[nthetabin];
Expand Down Expand Up @@ -81,7 +71,7 @@ static inline int countpairs_theta_mocks_fallback_DOUBLE(const int64_t N0, DOUBL
}

const DOUBLE *zstart = z1, *zend = z1 + N1;
const DOUBLE max_all_dz = SQRT(sqr_max_chord_sep - min_xdiff*min_xdiff - min_ydiff*min_ydiff);
const DOUBLE max_all_dz = SQRT(sep2max - min_xdiff*min_xdiff - min_ydiff*min_ydiff);
for(int64_t i=0;i<N0;i++) {
const DOUBLE xpos = *x0++;
const DOUBLE ypos = *y0++;
Expand Down Expand Up @@ -112,10 +102,10 @@ static inline int countpairs_theta_mocks_fallback_DOUBLE(const int64_t N0, DOUBL
const DOUBLE min_dy = min_ydiff > 0 ? min_ydiff + FABS(ypos - closest_icell_ypos):min_ydiff;
const DOUBLE min_dz = min_zdiff > 0 ? (this_dz > 0 ? this_dz:min_zdiff + FABS(zpos - closest_icell_zpos)):min_zdiff;
const DOUBLE sqr_min_sep_this_point = min_dx*min_dx + min_dy*min_dy + min_dz*min_dz;
if(sqr_min_sep_this_point >= sqr_max_chord_sep) {
if(sqr_min_sep_this_point >= sep2max) {
continue;
}
max_dz = SQRT(sqr_max_chord_sep - min_dx*min_dx - min_dy*min_dy);
max_dz = SQRT(sep2max - min_dx*min_dx - min_dy*min_dy);

const DOUBLE target_z = zpos - max_all_dz;
while(z1 != zend && *z1 <= target_z) {
Expand Down Expand Up @@ -175,27 +165,23 @@ static inline int countpairs_theta_mocks_fallback_DOUBLE(const int64_t N0, DOUBL
if(dz >= max_dz) break;

const DOUBLE sqr_chord_sep = dx*dx + dy*dy + dz*dz;
if(sqr_chord_sep >= sqr_max_chord_sep) {
if(sqr_chord_sep >= sep2max) {
continue;
}

DOUBLE theta = ZERO, pairweight = ZERO;
const DOUBLE one = (DOUBLE) 1.0, half = (DOUBLE) 0.5;

DOUBLE costheta = (one - half*sqr_chord_sep);
/* DOUBLE costheta = x2*xpos + y2*ypos + z2*zpos; */
if (costheta < -one) costheta = -one;
if (costheta > one) costheta = one;

if(costheta > costhetamin || costheta <= costhetamax) continue;
if(sqr_chord_sep >= sep2max || sqr_chord_sep < sep2min) continue;

if(need_rpavg) {
if(order) {
theta = INV_PI_OVER_180*FAST_ACOS(costheta);
} else {
theta = INV_PI_OVER_180*ACOS(costheta);
}
}
// if(need_rpavg) {
// // TODO
// if(order) {
// theta = INV_PI_OVER_180*FAST_ACOS(costheta);
// } else {
// theta = INV_PI_OVER_180*ACOS(costheta);
// }
// }
if(need_weightavg){
// These are only used for passing to weights
// Too expensive?
Expand All @@ -214,7 +200,7 @@ static inline int countpairs_theta_mocks_fallback_DOUBLE(const int64_t N0, DOUBL
}

for(int ibin=nthetabin-1;ibin>=1;ibin--) {
if(costheta <= costheta_upp[ibin-1]) {
if(sqr_chord_sep > sep2_upp[ibin-1]) {
npairs[ibin]++;
if(need_rpavg) {
thetaavg[ibin] += theta;
Expand Down