diff --git a/lalapps/src/pulsar/EinsteinAtHome/hs_boinc_extras.c b/lalapps/src/pulsar/EinsteinAtHome/hs_boinc_extras.c index 69a9d9ffc083f8e4dcda3a86d056fb6197803bc9..facda5bfa3e253141c3bd270dd59e61478892ca1 100644 --- a/lalapps/src/pulsar/EinsteinAtHome/hs_boinc_extras.c +++ b/lalapps/src/pulsar/EinsteinAtHome/hs_boinc_extras.c @@ -1143,8 +1143,8 @@ static int worker (void) { second_outfile = -1; } - /* record if there will be a second output file */ - else if (!strcmp("--SortToplist=6",argv[arg])) { + /* record if there will be triple output files */ + else if (!strcmp("--SortToplist=6",argv[arg]) || !strcmp("--SortToplist=7",argv[arg])) { rargv[rarg] = argv[arg]; LogPrintf(LOG_DEBUG,"BSGL output files\n"); bsgl_outfiles = -1; diff --git a/lalapps/src/pulsar/GCT/GCTtoplist.h b/lalapps/src/pulsar/GCT/GCTtoplist.h index 49b607a3f986dd10e0c2578bffe88b802646af87..7fce4548577f78b5cf268d988d7412a2ca2bd0ed 100644 --- a/lalapps/src/pulsar/GCT/GCTtoplist.h +++ b/lalapps/src/pulsar/GCT/GCTtoplist.h @@ -77,6 +77,7 @@ typedef enum SORTBY_BSGLtL = 4, //< sort by transient-line robust statistic B_S/GLtL SORTBY_BtSGLtL = 5, //< sort by transient-CW line robust statistic B_tS/GLtL SORTBY_TRIPLE_BStSGLtL = 6, //< triple toplists: one sorted by B_S/GL, one by B_S/GLtL, one by B_tS/GLtL + SORTBY_F_BSGLtL_BtSGLtL = 7,//< triple toplists: one sorted by 2F, one by B_S/GLtL, one by B_tS/GLtL SORTBY_LAST //< end-marker } SortBy_t; diff --git a/lalapps/src/pulsar/GCT/HierarchSearchGCT.c b/lalapps/src/pulsar/GCT/HierarchSearchGCT.c index 8f4aaa7bf1f4f40fe616458c772389e9d38e5184..c14a42c1d4cf873c200e22b708a7e877144b58b2 100644 --- a/lalapps/src/pulsar/GCT/HierarchSearchGCT.c +++ b/lalapps/src/pulsar/GCT/HierarchSearchGCT.c @@ -213,6 +213,7 @@ void PrintStackInfo( LALStatus *status, const SFTCatalogSequence *catalogSeq, FI void UpdateSemiCohToplists ( LALStatus *status, toplist_t *list1, toplist_t *list2, toplist_t *list3, FineGrid *in, REAL8 f1dot_fg, REAL8 f2dot_fg, REAL8 f3dot_fg, UsefulStageVariables *usefulparams, REAL4 NSegmentsInv, REAL4 *NSegmentsInvX, BOOLEAN have_f3dot ); void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, + SortBy_t toplist_sortby, toplist_t *list1, toplist_t *list2, toplist_t *list3, @@ -350,8 +351,8 @@ int MAIN( int argc, char *argv[]) { /* fstat candidate structure for candidate toplist*/ toplist_t *semiCohToplist=NULL; - toplist_t *semiCohToplist2=NULL; // only used for SORTBY_DUAL_F_BSGL or SORTBY_TRIPLE_BStSGLtL - toplist_t *semiCohToplist3=NULL; // only used for SORTBY_TRIPLE_BStSGLtL + toplist_t *semiCohToplist2=NULL; // only used for SORTBY_DUAL_F_BSGL, SORTBY_TRIPLE_BStSGLtL or SORTBY_F_BSGLtL_BtSGLtL + toplist_t *semiCohToplist3=NULL; // only used for SORTBY_TRIPLE_BStSGLtL and SORTBY_F_BSGLtL_BtSGLtL /* template and grid variables */ static DopplerSkyScanInit scanInit; /* init-structure for DopperScanner */ @@ -549,7 +550,7 @@ int MAIN( int argc, char *argv[]) { XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_oLGX, "oLGX", STRINGVector, 0, OPTIONAL, "BSGL: prior per-detector line-vs-Gauss odds 'oLGX' (Defaults to oLGX=1/Ndet)") == XLAL_SUCCESS, XLAL_EFUNC); XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_BSGLlogcorr, "BSGLlogcorr", BOOLEAN, 0, DEVELOPER, "BSGL: include log-correction terms (slower) or not (faster)") == XLAL_SUCCESS, XLAL_EFUNC); XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_getMaxFperSeg, "getMaxFperSeg", BOOLEAN, 0, OPTIONAL, "Compute and output maximum F and FX over segments") == XLAL_SUCCESS, XLAL_EFUNC); - XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_SortToplist, "SortToplist", INT4, 0, OPTIONAL, "Sort toplist by: 0=Fstat, 1=nc, 2=B_S/GL, 3='Fstat + B_S/GL', 4=B_S/GLtL, 5=B_tS/GLtL, 6='B_S/GL + B_S/GLtL + B_tS/GLtL'") == XLAL_SUCCESS, XLAL_EFUNC); + XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_SortToplist, "SortToplist", INT4, 0, OPTIONAL, "Sort toplist by: 0=Fstat, 1=nc, 2=B_S/GL, 3='Fstat + B_S/GL', 4=B_S/GLtL, 5=B_tS/GLtL, 6='B_S/GL + B_S/GLtL + B_tS/GLtL', 7='Fstat + B_S/GLtL + B_tS/GLtL' ") == XLAL_SUCCESS, XLAL_EFUNC); // -------------------------------------------- XLAL_CHECK_MAIN( XLALRegisterNamedUvarAuxData( &uvar_FstatMethod, "FstatMethod", UserEnum, XLALFstatMethodChoices(), 0, OPTIONAL, "F-statistic method to use" ) == XLAL_SUCCESS, XLAL_EFUNC); @@ -637,12 +638,14 @@ int MAIN( int argc, char *argv[]) { } if ( (uvar_SortToplist == SORTBY_BSGL || uvar_SortToplist == SORTBY_DUAL_F_BSGL || uvar_SortToplist == SORTBY_BSGLtL || uvar_SortToplist == SORTBY_BtSGLtL || - uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ) && !uvar_computeBSGL ) { - fprintf(stderr, "Toplist sorting by BSGL only possible if --computeBSGL given.\n"); + uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL || + uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL) && !uvar_computeBSGL ) { + fprintf(stderr, "Toplist sorting by BSGL[tL] only possible if --computeBSGL given.\n"); return( HIERARCHICALSEARCH_EBAD ); } if ( ( uvar_SortToplist == SORTBY_BSGLtL || uvar_SortToplist == SORTBY_BtSGLtL || - uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ) && !uvar_getMaxFperSeg ) { + uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL || + uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL ) && !uvar_getMaxFperSeg ) { fprintf(stderr, "Toplist sorting by B[t]SGLtL only possible if --getMaxFperSeg given.\n"); return( HIERARCHICALSEARCH_EBAD ); } @@ -665,6 +668,15 @@ int MAIN( int argc, char *argv[]) { XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist3, uvar_nCand1, SORTBY_BtSGLtL ), XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BtSGLtL ); } + else if ( uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL )// special treatement of 'triple' toplists: 1st one sorted by 'F', 2nd one by 'B_S/GLtL', 3rd by 'B_tS/GLtL' + { + XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist, uvar_nCand1, SORTBY_F ), + XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_F ); + XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist2, uvar_nCand1, SORTBY_BSGLtL ), + XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGLtL ); + XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist3, uvar_nCand1, SORTBY_BtSGLtL ), + XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BtSGLtL ); + } else // 'normal' single-sorting toplist cases (sortby 'F', 'nc' or 'BSGL') { XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist, uvar_nCand1, uvar_SortToplist), @@ -1761,8 +1773,9 @@ int MAIN( int argc, char *argv[]) { /* this is necessary here, because UpdateSemiCohToplists() might set a checkpoint that needs some information from here */ - if(uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL && finegrid.sumTwoFX && finegrid.maxTwoFXl ) { - LAL_CALL( UpdateSemiCohToplistsOptimTriple (&status, semiCohToplist, semiCohToplist2, semiCohToplist3, &finegrid, f1dot_fg, f2dot_fg, f3dot_fg, &usefulParams, NSegmentsInv, usefulParams.NSegmentsInvX, XLALUserVarWasSet(&uvar_f3dot) ), &status); + if(( uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL || uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL ) && + finegrid.sumTwoFX && finegrid.maxTwoFXl ) { + LAL_CALL( UpdateSemiCohToplistsOptimTriple (&status, uvar_SortToplist, semiCohToplist, semiCohToplist2, semiCohToplist3, &finegrid, f1dot_fg, f2dot_fg, f3dot_fg, &usefulParams, NSegmentsInv, usefulParams.NSegmentsInvX, XLALUserVarWasSet(&uvar_f3dot) ), &status); } else { LAL_CALL( UpdateSemiCohToplists (&status, semiCohToplist, semiCohToplist2, semiCohToplist3, &finegrid, f1dot_fg, f2dot_fg, f3dot_fg, &usefulParams, NSegmentsInv, usefulParams.NSegmentsInvX, XLALUserVarWasSet(&uvar_f3dot) ), &status); } @@ -1941,7 +1954,8 @@ int MAIN( int argc, char *argv[]) { strcpy(t3_suffix,""); break; }; - case SORTBY_TRIPLE_BStSGLtL : { + case SORTBY_TRIPLE_BStSGLtL: + case SORTBY_F_BSGLtL_BtSGLtL : { strcpy(t1_suffix,""); strcpy(t2_suffix,"-BSGLtL"); strcpy(t3_suffix,"-BtSGLtL"); @@ -2674,6 +2688,7 @@ void PrintStackInfo( LALStatus *status, * This function allows for inserting candidates into up to 3 toplists at once, which might be sorted differently! */ void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, + SortBy_t toplists_sortby, toplist_t *list1, toplist_t *list2, //< optional (can be NULL): insert candidate into this 2nd toplist as well toplist_t *list3, //< optional (can be NULL): insert candidate into this 3rd toplist as well @@ -2691,6 +2706,7 @@ void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, REAL8 freq_fg; UINT4 ifreq_fg; GCTtopOutputEntry line; + BOOLEAN delay_compute_BSGL = (toplists_sortby == SORTBY_F_BSGLtL_BtSGLtL); INITSTATUS(status); ATTATCHSTATUSPTR (status); @@ -2699,9 +2715,11 @@ void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, ASSERT ( in != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL ); ASSERT ( usefulparams != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL ); - /* Optimized version for triple toplist case: + /* Optimized version for triple toplist cases: BSGL, BSGLtL,BtSGLtL or + 2F, BSGLtL,BtSGLtL - First compute all three detection metrics by which any of the toplists is sorted. + First compute all detection metrics by which any of the toplists is sorted (even as seconary + comparison criterion). Next test if the candidate is loud enough to make it into at least one of the toplists. If so, build the candidate structure. @@ -2730,7 +2748,15 @@ void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, } xlalErrno = 0; - line.log10BSGL = XLALComputeBSGL ( sumTwoF, sumTwoFX, usefulparams->BSGLsetup ); + /* Note: sorting a toplist by 2F takes BSGL as a secondary sorting criterion, so even if the 1st toplist is + sorted by 2F, we still provide an overestimate of BSGL for a test if it can make it into the toplist in + the first place */ + + if(delay_compute_BSGL) { + line.log10BSGL = LAL_REAL4_MAX; + } else { + line.log10BSGL = XLALComputeBSGL ( sumTwoF, sumTwoFX, usefulparams->BSGLsetup ); + } if ( xlalErrno != 0 ) { XLALPrintError ("%s line %d : XLALComputeBSGL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno ); ABORT ( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL ); @@ -2769,8 +2795,8 @@ void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, /* take F-stat averages over segments */ line.avTwoF = sumTwoF*NSegmentsInv; /* average multi-2F by full number of segments */ -/* now test if this candidate makes it into any of the toplists, if not we don't need - to do any more copying of data from finegrid to toplist entry structure */ + /* now test if this candidate makes it into any of the toplists, if not we don't need + to do any more copying and computing of data from finegrid to toplist entry structure */ int isIncludedToplists = TEST_FSTAT_TOPLIST_INCLUSION( list1, &line) || @@ -2780,6 +2806,12 @@ void UpdateSemiCohToplistsOptimTriple ( LALStatus *status, if(likely(! isIncludedToplists)) continue; + + /* if we delayed the computation of BSGL for a 2F sorted first toplist, now is the time to actually compute it */ + if(delay_compute_BSGL) { + line.log10BSGL = XLALComputeBSGL ( sumTwoF, sumTwoFX, usefulparams->BSGLsetup ); + } + line.numDetectors = in->numDetectors; line.avTwoFrecalc = -1.0; /* initialise this to -1.0, so that it only gets written out by print_gctFstatline_to_str if later overwritten in recalcToplistStats step */ line.log10BSGLrecalc = -LAL_REAL4_MAX; /* for now, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */ diff --git a/lalapps/src/pulsar/GCT/testGCT.sh b/lalapps/src/pulsar/GCT/testGCT.sh index b474e2bf83d4b0040f64a6427610946d1d81aa58..a654d904902b2112aa0065d569f1ae25e7a7bbd1 100755 --- a/lalapps/src/pulsar/GCT/testGCT.sh +++ b/lalapps/src/pulsar/GCT/testGCT.sh @@ -389,8 +389,8 @@ fi echo echo "----------------------------------------------------------------------------------------------------" -echo " STEP 7: run HierarchSearchGCT using Resampling and triple toplist (BSGL, BSGLtL and BtSGLtL), " -echo " compared with separate runs for each of these 3 toplist rankings" +echo " STEP 7: run HierarchSearchGCT using Resampling and triple toplists (BSGL, BSGLtL and BtSGLtL), " +echo " (2F, BSGLtL and BtSGLtL), compared with separate runs for each of these 3 toplist rankings" echo "----------------------------------------------------------------------------------------------------" echo @@ -398,6 +398,8 @@ rm -f checkpoint.cpt # delete checkpoint to start correctly outfile_GCT_RS_triple="${testDir}/GCT_RS_triple.dat" timingsfile_RS_triple="${testDir}/timing_RS_triple.dat" +outfile_GCT_RS_triple2="${testDir}/GCT_RS_triple2.dat" +timingsfile_RS_triple2="${testDir}/timing_RS_triple2.dat" # set plan mode so that results should be deterministic, easier to compare results of different runs this way @@ -419,7 +421,42 @@ if ! eval "$cmdline"; then fi -## re-run, but now create the three toplists one by one, then compare results +# second variant of triple toplists, with 2F sorted first toplist + +cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple2' --outputTiming='$timingsfile_RS_triple2' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=7" +if [ -n "$DEBUG" ]; then + cmdline="$cmdline" +else + cmdline="$cmdline &> /dev/null" +fi + +echo $cmdline +if ! eval "$cmdline"; then + echo "Error.. something failed when running '$gct_code' ..." + exit 1 +fi + + + + +## re-run, but now create the in total four toplists one by one, then compare results + +outfile_GCT_RS_triple="${testDir}/GCT_RS_triple_0.dat" + +cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=0" +if [ -n "$DEBUG" ]; then + cmdline="$cmdline" +else + cmdline="$cmdline &> /dev/null" +fi + +echo $cmdline +if ! eval "$cmdline"; then + echo "Error.. something failed when running '$gct_code' ..." + exit 1 +fi + + outfile_GCT_RS_triple="${testDir}/GCT_RS_triple_1.dat" @@ -469,6 +506,7 @@ fi # filter out comments +egrep -v "^%" ${testDir}/GCT_RS_triple_0.dat > ${testDir}/GCT_RS_triple_0.txt egrep -v "^%" ${testDir}/GCT_RS_triple_1.dat > ${testDir}/GCT_RS_triple_1.txt egrep -v "^%" ${testDir}/GCT_RS_triple_2.dat > ${testDir}/GCT_RS_triple_2.txt egrep -v "^%" ${testDir}/GCT_RS_triple_3.dat > ${testDir}/GCT_RS_triple_3.txt @@ -478,6 +516,11 @@ egrep -v "^%" ${testDir}/GCT_RS_triple.dat > ${testDir}/GCT_RS_triple.txt egrep -v "^%" ${testDir}/GCT_RS_triple.dat-BSGLtL > ${testDir}/GCT_RS_triple-BSGLtL.txt egrep -v "^%" ${testDir}/GCT_RS_triple.dat-BtSGLtL > ${testDir}/GCT_RS_triple-BtSGLtL.txt +egrep -v "^%" ${testDir}/GCT_RS_triple2.dat > ${testDir}/GCT_RS_triple2.txt +egrep -v "^%" ${testDir}/GCT_RS_triple2.dat-BSGLtL > ${testDir}/GCT_RS_triple2-BSGLtL.txt +egrep -v "^%" ${testDir}/GCT_RS_triple2.dat-BtSGLtL > ${testDir}/GCT_RS_triple2-BtSGLtL.txt + + if ! eval "diff ${testDir}/GCT_RS_triple_1.txt ${testDir}/GCT_RS_triple.txt"; then echo "Error: tripple toplists do not match separately generated toplists (1) " exit 1 @@ -493,6 +536,24 @@ if ! eval "diff ${testDir}/GCT_RS_triple_3.txt ${testDir}/GCT_RS_triple-BtSGLtL. exit 1 fi +if ! eval "diff ${testDir}/GCT_RS_triple_0.txt ${testDir}/GCT_RS_triple2.txt"; then + echo "Error: tripple toplists do not match separately generated toplists (4) " + exit 1 +fi + +if ! eval "diff ${testDir}/GCT_RS_triple_2.txt ${testDir}/GCT_RS_triple2-BSGLtL.txt"; then + echo "Error: tripple toplists do not match separately generated toplists (5) " + exit 1 +fi + +if ! eval "diff ${testDir}/GCT_RS_triple_3.txt ${testDir}/GCT_RS_triple2-BtSGLtL.txt"; then + echo "Error: tripple toplists do not match separately generated toplists (6) " + exit 1 +fi + + + + ## ---------- compute relative differences and check against tolerance -------------------- diff --git a/lalpulsar/src/SSBtimes.c b/lalpulsar/src/SSBtimes.c index 1a254bcd6655f6e034020fc20ea24bb59129a41c..2c39070be7e028bd20ce5c3d79d12a12ea6d4b95 100644 --- a/lalpulsar/src/SSBtimes.c +++ b/lalpulsar/src/SSBtimes.c @@ -42,6 +42,7 @@ const UserChoices SSBprecisionChoices = { { SSBPREC_NEWTONIAN, "newtonian" }, { SSBPREC_RELATIVISTIC, "relativistic" }, { SSBPREC_RELATIVISTICOPT, "relativisticopt" }, + { SSBPREC_DMOFF, "DMoff"}, }; /*---------- internal prototypes ----------*/ @@ -637,6 +638,16 @@ XLALGetSSBtimes ( const DetectorStateSeries *DetectorStates, /**< [in] detector- break; + case SSBPREC_DMOFF: /* switch off all demodulation terms */ + + for ( UINT4 i = 0; i < numSteps; i++ ) + { + DetectorState *state = &(DetectorStates->data[i]); + ret->DeltaT->data[i] = XLALGPSGetREAL8 ( &state->tGPS ) - refTimeREAL8; + ret->Tdot->data[i] = 1.0; + } /* for i < numSteps */ + break; + default: XLAL_ERROR_NULL (XLAL_EFAILED, "\n?? Something went wrong.. this should never be called!\n\n" ); break; diff --git a/lalpulsar/src/SSBtimes.h b/lalpulsar/src/SSBtimes.h index 05e51f703f38fa9f2be32f06f0822a8adef333d1..0e49ebb14fafec91433b83d6594ca5f029b0cd9b 100644 --- a/lalpulsar/src/SSBtimes.h +++ b/lalpulsar/src/SSBtimes.h @@ -46,6 +46,7 @@ typedef enum tagSSBprecision { SSBPREC_NEWTONIAN, /**< simple Newtonian: \f$\tau = t + \vec{r}\cdot\vec{n}/c\f$ */ SSBPREC_RELATIVISTIC, /**< detailed relativistic: \f$\tau=\tau(t; \vec{n}, \vec{r})\f$ */ SSBPREC_RELATIVISTICOPT, /**< optimized relativistic, numerically equivalent to #SSBPREC_RELATIVISTIC, but faster */ + SSBPREC_DMOFF, /**< switch off all demodulatoin terms */ SSBPREC_LAST /**< end marker */ } SSBprecision;