20#include <gsl/gsl_randist.h>
21#include <gsl/gsl_cdf.h>
47 UINT4 numberofmaxima = fbins * rows - fbins - (
UINT4 )( ( rows * rows - rows ) / 2 );
54 ihsmaxima->
rows = rows;
110 REAL8 dailyharmonic =
params->Tobs / ( 24.0 * 3600.0 );
111 REAL8 siderealharmonic =
params->Tobs / 86164.0905;
112 REAL8 dailyharmonic2 = dailyharmonic * 2.0, dailyharmonic3 = dailyharmonic * 3.0, dailyharmonic4 = dailyharmonic * 4.0;
113 REAL8 siderealharmonic2 = siderealharmonic * 2.0, siderealharmonic3 = siderealharmonic * 3.0, siderealharmonic4 = siderealharmonic * 4.0;
116 memset( markedharmonics->
data, 0,
sizeof(
INT4 )*markedharmonics->
length );
118 if ( !
params->noNotchHarmonics ) {
119 for (
UINT4 ii = 0; ii < markedharmonics->
length; ii++ ) {
120 if ( fabs( dailyharmonic - (
REAL8 )ii ) <= 1.0 || fabs( dailyharmonic2 - (
REAL8 )ii ) <= 1.0 || fabs( dailyharmonic3 - (
REAL8 )ii ) <= 1.0 || fabs( dailyharmonic4 - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic2 - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic3 - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic4 - (
REAL8 )ii ) <= 1.0 ) {
121 markedharmonics->
data[ii] = 1;
127 for (
UINT4 ii = 0; ii < ihss->length; ii++ ) {
130 memcpy(
row->data, &( input->
ffdata->
data[ii * numfprbins] ),
sizeof(
REAL4 )*numfprbins );
131 if ( !
params->noNotchHarmonics )
for (
UINT4 jj = 0; jj <
row->length; jj++ )
if ( markedharmonics->
data[jj] == 1 ) {
136 if ( !
params->weightedIHS ) {
143 memcpy( ihsvectorarray->
data[ii]->
data, ihsvector->data,
sizeof(
REAL4 )*ihsvector->length );
207 for (
UINT4 ii = 0; ii < tempvect->
length; ii++ ) {
236 for (
INT4 ii = 5; ii < highval; ii++ ) {
237 output->data[ii - 5] = 0.0;
238 for (
UINT4 jj = 1; jj <= ihsfactor; jj++ ) {
267 for (
INT4 ii = 5; ii < highval; ii++ ) {
268 output->data[ii - 5] = 0.0;
269 for (
UINT4 jj = 1; jj <= ihsfactor; jj++ ) {
270 REAL4 weight = 1.0 / ( aveNoise->
data[ii * jj] * aveNoise->
data[ii * jj] );
271 output->data[ii - 5] += input->
data[ii * jj] * weight;
318 if ( ihsfarstruct ) {
345 fprintf( stderr,
"Determining IHS FAR values... " );
346 fprintf(
LOG,
"Determining IHS FAR values... " );
350 UINT4 trials = 5 * rows;
351 if ( trials < 1000 ) {
354 if ( trials > 5000 ) {
355 fprintf( stderr,
"Warning: number of trials may be insufficient given the number of rows to sum\n" );
371 REAL8 dailyharmonic = Tobs / ( 24.0 * 3600.0 );
372 REAL8 siderealharmonic = Tobs / 86164.0905;
373 REAL8 dailyharmonic2 = dailyharmonic * 2.0, dailyharmonic3 = dailyharmonic * 3.0, dailyharmonic4 = dailyharmonic * 4.0;
374 REAL8 siderealharmonic2 = siderealharmonic * 2.0, siderealharmonic3 = siderealharmonic * 3.0, siderealharmonic4 = siderealharmonic * 4.0;
377 memset( markedharmonics->
data, 0,
sizeof(
INT4 )*markedharmonics->
length );
378 if ( !
params->noNotchHarmonics ) {
379 for (
UINT4 ii = 0; ii < markedharmonics->
length; ii++ ) {
380 if ( fabs( dailyharmonic - (
REAL8 )ii ) <= 1.0 || fabs( dailyharmonic2 - (
REAL8 )ii ) <= 1.0 || fabs( dailyharmonic3 - (
REAL8 )ii ) <= 1.0 || fabs( dailyharmonic4 - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic2 - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic3 - (
REAL8 )ii ) <= 1.0 || fabs( siderealharmonic4 - (
REAL8 )ii ) <= 1.0 ) {
381 markedharmonics->
data[ii] = 1;
387 memcpy( noise->data, aveNoise->
data,
sizeof(
REAL4 )*aveNoise->
length );
388 for (
UINT4 ii = 0; ii < aveNoise->
length; ii++ )
if ( markedharmonics->
data[ii] == 1 ) {
389 noise->data[ii] = 0.0;
391 if ( !
params->weightedIHS ) {
398 for (
UINT4 ii = 0; ii < trials; ii++ ) {
400 for (
UINT4 jj = 0; jj < aveNoise->
length; jj++ ) {
401 if ( markedharmonics->
data[jj] == 0 ) {
402 noise->data[jj] = (
REAL4 )( gsl_ran_exponential( rng, aveNoise->
data[jj] ) );
404 noise->data[jj] = 0.0;
411 randval = 1.0 + gsl_ran_gaussian( rng, 0.2 );
412 while ( randval <= 0.0 || randval >= 2.0 ) {
413 randval = 1.0 + gsl_ran_gaussian( rng, 0.2 );
420 if ( !
params->weightedIHS ) {
438 for (
UINT4 ii = 0; ii < trials; ii++ ) {
439 FbinMean->
data[ii] = 1.0;
498 for (
UINT4 ii = 0; ii < ihsvalues->
length; ii++ ) {
506 ihsvalues->
data[ii] = ihsvectorarray->
data[ii]->
data[ihslocations->
data[ii] - 5];
516 for (
UINT4 ii = 2; ii <= rows; ii++ ) {
523 if (
params->vectorMath != 0 ) {
532 for (
UINT4 jj = 0; jj < ihsvectorarray->
length - ( ii - 1 ); jj++ ) {
534 if (
params->vectorMath == 0 ) {
535 if ( ii > 2 )
for (
UINT4 kk = 0; kk < tworows->
data[jj]->
length; kk++ ) {
543 memcpy( rowarraylocs->
data, &( ihslocations->
data[jj] ),
sizeof(
INT4 )*ii );
552 if ( ( ihsvectorarray->
length - ( ii - 1 ) ) * ihsvectorarray->
data[0]->
length > 10000 ) {
562 for (
UINT4 jj = 0; jj < ( ihsvectorarray->
length - ( ii - 1 ) ); jj++ ) {
572 if (
params->ihsfar != 1.0 ) {
574 for (
UINT4 jj = 0; jj < sampledtempihsvals->
length; jj++ ) {
576 if ( !
params->fastchisqinv && sampledtempihsvals->
data[jj] != 0.0 ) {
578 farave += 0.5 * gsl_cdf_chisq_Qinv(
params->ihsfar, 2.0 * sampledtempihsvals->
data[jj] );
579 }
else if (
params->fastchisqinv && sampledtempihsvals->
data[jj] != 0.0 ) {
597 if (
params->ihsfomfar != 1.0 &&
params->ihsfom == 0.0 ) {
603 }
else if (
params->ihsfom != 0.0 ) {
653 REAL4VectorAligned *ihsvalues = NULL, *excessabovenoise = NULL, *scaledExpectedIHSVectorValues = NULL;
665 for (
UINT4 ii = 0; ii < ihsvalues->
length; ii++ ) {
672 for (
INT4 jj = 0; jj <
params->harmonicNumToSearch; jj++ ) {
675 ihsvalues->
data[ii] = ihsvectorarray->
data[ii]->
data[ihslocations->
data[ii] - 5];
677 INT4 newIHSlocation =
max_index_in_range( excessabovenoise, ( jj + 1 ) * minIndexForIHS, ( jj + 1 ) * maxIndexForIHS ) + 5;
678 REAL4 newIHSvalue = ihsvectorarray->
data[ii]->
data[newIHSlocation - 5];
679 if ( newIHSvalue > ihsvalues->
data[ii] ) {
680 ihslocations->
data[ii] = newIHSlocation;
681 ihsvalues->
data[ii] = newIHSvalue;
691 for (
UINT4 ii = 1; ii <= rows; ii++ ) {
695 memcpy(
output->locationsForEachFbin->data, ihslocations->
data,
sizeof(
INT4 )*ihslocations->
length );
703 REAL4 sumofnoise = 0.0;
704 INT4 endloc = ( ( ii - 1 ) * ( ii - 1 ) - ( ii - 1 ) ) / 2;
707 if (
params->vectorMath != 0 ) {
716 for (
UINT4 jj = 0; jj < ihsvectorarray->
length - ( ii - 1 ); jj++ ) {
718 if (
params->vectorMath == 0 ) {
719 if ( ii > 2 )
for (
UINT4 kk = 0; kk < tworows->
data[jj]->
length; kk++ ) {
727 if ( jj == 0 )
for (
UINT4 kk = 0; kk < ii; kk++ ) {
728 sumofnoise += FbinMean->
data[kk];
730 sumofnoise -= FbinMean->
data[jj - 1];
731 sumofnoise += FbinMean->
data[jj + ( ii - 1 )];
740 for (
INT4 kk = 0; kk <
params->harmonicNumToSearch; kk++ ) {
743 output->maxima->data[( ii - 2 )*ihsvalues->
length - endloc + jj] = tworows->
data[jj]->
data[(
output->locations->data[( ii - 2 ) * ihsvalues->
length - endloc + jj] - 5 )];
745 INT4 newIHSlocation =
max_index_in_range( excessabovenoise, ( kk + 1 ) * minIndexForIHS, ( kk + 1 ) * maxIndexForIHS ) + 5;
746 REAL4 newIHSvalue = tworows->
data[ii]->
data[newIHSlocation - 5];
747 if ( newIHSvalue >
output->maxima->data[( ii - 2 ) * ihsvalues->
length - endloc + jj] ) {
748 output->locations->data[( ii - 2 )*ihsvalues->
length - endloc + jj] = newIHSlocation;
749 output->maxima->data[( ii - 2 )*ihsvalues->
length - endloc + jj] = newIHSvalue;
754 memcpy( rowarraylocs->
data, &( ihslocations->
data[jj] ),
sizeof(
INT4 )*ii );
815 INT4 numberofIHSvalsChecked = 0, numberofIHSvalsExceededThresh = 0, numberPassingBoth = 0, linesinterferewithnum = 0, skipped = 0, notskipped = 0;
830 REAL8 highestval = 0.0, highestsignificance = 0.0;
831 INT4 highestvalloc = -1, jjloc = 0;
832 for (
UINT4 jj = 0; jj < numfbins - ( ii - 1 ); jj++ ) {
835 memcpy( avgsinrange->data, &( fbinavgs->
data[jj] ),
sizeof(
REAL4 )*ii );
838 numberofIHSvalsChecked++;
840 INT4 locationinmaximastruct = ( ii - 2 ) * numfbins - ( ( ii - 1 ) * ( ii - 1 ) - ( ii - 1 ) ) / 2 + jj;
843 if ( ihsmaxima->
maxima->
data[locationinmaximastruct] > ihsfarstruct->
ihsfar->
data[ii - 2]*meanNoise ) {
845 numberofIHSvalsExceededThresh++;
853 fsig =
params->fmin -
params->dfmax + ( ( 0.5 * ( ii - 1 ) + jj ) - 6.0 ) /
params->Tsft;
854 B = 0.5 * ( ii - 1 ) /
params->Tsft;
857 INT4 nolinesinterfering = 1;
858 if ( trackedlines != NULL ) {
860 while ( kk < trackedlines->length && nolinesinterfering == 1 ) {
861 if ( 2.0 *
B >= ( trackedlines->
data[kk * 3 + 2] - trackedlines->
data[kk * 3 + 1] ) ) {
862 if ( ( trackedlines->
data[kk * 3 + 2] >= (
REAL4 )( fsig -
B ) && trackedlines->
data[kk * 3 + 2] <= (
REAL4 )( fsig +
B ) ) ||
863 ( trackedlines->
data[kk * 3 + 1] >= (
REAL4 )( fsig -
B ) && trackedlines->
data[kk * 3 + 1] <= (
REAL4 )( fsig +
B ) ) ) {
864 nolinesinterfering = 0;
868 if ( ( (
REAL4 )( fsig +
B ) >= trackedlines->
data[kk * 3 + 1] && (
REAL4 )( fsig +
B ) <= trackedlines->
data[kk * 3 + 2] ) ||
869 ( (
REAL4 )( fsig -
B ) >= trackedlines->
data[kk * 3 + 1] && (
REAL4 )( fsig -
B ) <= trackedlines->
data[kk * 3 + 2] ) ) {
870 nolinesinterfering = 0;
877 if ( !nolinesinterfering ) {
878 linesinterferewithnum++;
881 REAL8 totalnoise = meanNoise * noise;
883 REAL8 significance = gsl_cdf_chisq_Q( 2.0 * ihsmaxima->
maxima->
data[locationinmaximastruct], 2.0 * totalnoise );
884 if ( significance == 0.0 ) {
885 significance = log10(
LAL_E ) * ihsmaxima->
maxima->
data[locationinmaximastruct] - ( totalnoise - 1.0 ) * log10( ihsmaxima->
maxima->
data[locationinmaximastruct] ) + lgamma( totalnoise ) /
log( 10.0 );
887 significance = -log10( significance );
890 if ( significance > highestsignificance && (
params->followUpOutsideULrange || ( !
params->followUpOutsideULrange && fsig >=
params->ULfmin && fsig < (
params->ULfmin +
params->ULfspan ) &&
B >=
params->ULminimumDeltaf && B <= params->ULmaximumDeltaf ) ) ) {
891 highestval = ihsmaxima->
maxima->
data[locationinmaximastruct] - totalnoise;
892 highestvalloc = locationinmaximastruct;
893 highestsignificance = significance;
905 if ( highestvalloc != -1 ) {
907 fsig =
params->fmin -
params->dfmax + ( 0.5 * ( ii - 1 ) + jjloc - 6.0 ) /
params->Tsft;
908 B = 0.5 * ( ii - 1 ) /
params->Tsft;
912 if ( ( *candlist )->numofcandidates == ( *candlist )->length - 1 ) {
915 loadCandidateData( &( ( *candlist )->data[( *candlist )->numofcandidates] ), fsig, per0,
B, 0.0, 0.0, ihsmaxima->
maxima->
data[highestvalloc],
h0, highestsignificance, 0, ffdata->
tfnormalization, -1, 0 );
916 ( *candlist )->numofcandidates++;
1067 fprintf( stderr,
"Number of IHS vals checked = %d, number exceeding IHS threshold = %d, number passing both = %d, but lines interfere with %d, number not skipped = %d and number of skipped candidates is %d\n", numberofIHSvalsChecked, numberofIHSvalsExceededThresh, numberPassingBoth, linesinterferewithnum, notskipped, skipped );
1068 fprintf( stderr,
"Candidates found in IHS step = %d\n", ( *candlist )->numofcandidates );
1069 fprintf(
LOG,
"Number of IHS vals checked = %d, number exceeding IHS threshold = %d, number passing both = %d, but lines interfere with %d, number not skipped = %d and number of skipped candidates is %d\n", numberofIHSvalsChecked, numberofIHSvalsExceededThresh, numberPassingBoth, linesinterferewithnum, notskipped, skipped );
1070 fprintf(
LOG,
"Candidates found in IHS step = %d\n", ( *candlist )->numofcandidates );
1086 if ( ihsval <= 0.0 ) {
1089 REAL8 prefact = 1.0;
1092 return prefact * pow( ihsval / (
params->Tsft *
params->Tobs ), 0.25 );
void destroyihsfarStruct(ihsfarStruct *ihsfarstruct)
Destroy ihsfarStruct struct.
void destroyihsMaxima(ihsMaximaStruct *data)
Destroy vectors and the IHS maxima struct.
ihsVals * createihsVals(void)
Allocate memory for ihsVals struct.
REAL8 ihs2h0(const REAL8 ihsval, const UserInput_t *params)
Convert the IHS statistic to an estimated h0, based on injections.
INT4 incHarmSumVector(REAL4VectorAligned *output, const REAL4VectorAligned *input, const UINT4 ihsfactor)
Compute the IHS vector – does not compute the maximum value.
INT4 runIHS(ihsMaximaStruct *output, const ffdataStruct *input, const ihsfarStruct *ihsfarinput, const UserInput_t *params, const UINT4 rows, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *FbinMean)
Run the IHS algorithm.
void destroyihsVals(ihsVals *ihsvals)
Destroy ihsVals struct.
INT4 incHarmSum(ihsVals *output, const REAL4VectorAligned *input, const UINT4 ihsfactor)
Compute the IHS sum maximum.
ihsfarStruct * createihsfarStruct(const UINT4 rows, const UserInput_t *params)
Allocate memory for ihsfarStruct struct.
INT4 incHarmSumVectorWeighted(REAL4VectorAligned *output, const REAL4VectorAligned *input, const REAL4VectorAligned *aveNoise, const UINT4 ihsfactor)
Compute the noise weighted IHS vector – does not compute the maximum value.
REAL4 ihsFOM(const INT4Vector *locs, const INT4 fomnorm)
Calculate the IHS FOM for a number of rows.
INT4 sumIHSarray(ihsMaximaStruct *output, const ihsfarStruct *inputfar, REAL4VectorAlignedArray *ihsvectorarray, const UINT4 rows, const REAL4VectorAligned *FbinMean, const UserInput_t *params)
Compute the IHS sums for a number of rows.
INT4 genIhsFar(ihsfarStruct *output, const UserInput_t *params, const UINT4 rows, const REAL4VectorAligned *aveNoise, const gsl_rng *rng)
Compute the IHS FAR for a sum of a number of rows.
INT4 sumIHSarrayFAR(ihsfarStruct *outputfar, REAL4VectorAlignedArray *ihsvectorarray, const UINT4 rows, const REAL4VectorAligned *FbinMean, const UserInput_t *params, const gsl_rng *rng)
Compute the IHS sums for a number of rows used for the FAR calculation.
INT4 findIHScandidates(candidateVector **candlist, const ihsfarStruct *ihsfarstruct, const UserInput_t *params, const ffdataStruct *ffdata, const ihsMaximaStruct *ihsmaxima, const REAL4VectorAligned *fbinavgs, const REAL4VectorSequence *trackedlines)
Find IHS candidates above thresholds.
ihsMaximaStruct * createihsMaxima(const UINT4 fbins, const UINT4 rows)
Create vectors for IHS maxima struct.
INT4 sort_float_smallest(REAL4VectorAligned *output, const REAL4VectorAligned *input)
Sort a REAL4VectorAligned, keeping the smallest of the values in the output vector.
REAL4 calcMean_ignoreZeros(const REAL4VectorAligned *vector)
Compute the mean value of a REAL4VectorAligned without accepting values of zero.
INT4 calcStddev_ignoreZeros(REAL4 *sigma, const REAL4VectorAligned *vector)
Compute the standard deviation of a REAL4VectorAligned ignoring zero values.
INT4 calcStddev(REAL4 *sigma, const REAL4VectorAligned *vector)
Compute the standard deviation of a REAL4VectorAligned.
UINT4 max_index_in_range(const REAL4VectorAligned *vector, const UINT4 startlocation, const UINT4 lastlocation)
Determine the index value of the maximum value between elements of a REAL4VectorAligned (inclusive)
REAL4 calcMean(const REAL4VectorAligned *vector)
Compute the mean value of a REAL4VectorAligned, computed via recursion like in GSL.
REAL4VectorAligned * sampleREAL4VectorAlignedArray_nozerosaccepted(const REAL4VectorAlignedArray *input, const UINT4 numberofvectors, const UINT4 sampleSize, const gsl_rng *rng)
Sample a number (sampleSize) of values from an REAL4VectorAlignedArray (input) randomly from vector 0...
REAL8 minPeriod(const REAL8 moddepth, const REAL8 cohtime)
Calculates minimum period allowed, equation 6 of E.
candidateVector * resizecandidateVector(candidateVector *vector, const UINT4 length)
Resize a candidateVector.
void loadCandidateData(candidate *output, const REAL8 fsig, const REAL8 period, const REAL8 moddepth, const REAL4 ra, const REAL4 dec, const REAL8 statval, const REAL8 h0, const REAL8 prob, const INT4 proberrcode, const REAL8 normalization, const INT4 templateVectorIndex, const BOOLEAN lineContamination)
Load candidate data.
REAL8 cdf_chisq_Qinv(REAL8 Q, REAL8 nu)
void * XLALMalloc(size_t n)
INT4Vector * XLALCreateINT4Vector(UINT4 length)
void XLALDestroyINT4Vector(INT4Vector *vector)
int XLALVectorScaleREAL4(REAL4 *out, REAL4 scalar, const REAL4 *in, const UINT4 len)
void XLALDestroyREAL4VectorAligned(REAL4VectorAligned *in)
REAL4VectorAligned * XLALCreateREAL4VectorAligned(const UINT4 length, const UINT4 align)
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
#define XLAL_CHECK_NULL(assertion,...)
REAL4VectorAligned ** data
REAL4VectorAligned * ffdata
INT4Vector * locationsForEachFbin
REAL4VectorAligned * maxima
REAL4VectorAligned * maximaForEachFbin
REAL4VectorAligned * foms
REAL4VectorAligned * ihsfomdistMean
REAL4VectorAligned * ihsdistMean
REAL4VectorAligned * expectedIHSVector
REAL4VectorAligned * ihsfomdistSigma
REAL4VectorAligned * ihsfar
REAL4VectorAligned * ihsdistSigma
REAL4VectorAligned * fomfarthresh
INT4 VectorArraySum(REAL4VectorAlignedArray *output, REAL4VectorAlignedArray *input1, REAL4VectorAlignedArray *input2, INT4 vectorpos1, INT4 vectorpos2, INT4 outputvectorpos, INT4 numvectors)
Sum vectors from REAL4VectorAlignedArrays into an output REAL4VectorAlignedArray using SIMD.
REAL4VectorAlignedArray * createREAL4VectorAlignedArray(const UINT4 length, const UINT4 vectorLength, const size_t align)
void destroyREAL4VectorAlignedArray(REAL4VectorAlignedArray *array)
INT4 VectorSubtractREAL4(REAL4VectorAligned *output, REAL4VectorAligned *input1, REAL4VectorAligned *input2, INT4 vectorMath)