Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
IHS.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2010--2012, 2014, 2015 Evan Goetz
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20#include <gsl/gsl_randist.h>
21#include <gsl/gsl_cdf.h>
22
23#include "IHS.h"
24#include "statistics.h"
25#include "candidates.h"
26#include "cdfdist.h"
27#include "vectormath.h"
28#include "TwoSpect.h"
29
30
31/**
32 * Create vectors for IHS maxima struct
33 * \param [in] fbins Number of frequency bins
34 * \param [in] rows Number of neighboring rows to be summed
35 * \return Pointer to a newly created ihsMaximaStruct
36 */
37ihsMaximaStruct *createihsMaxima( const UINT4 fbins, const UINT4 rows )
38{
39
40 XLAL_CHECK_NULL( fbins > 0 && rows > 0, XLAL_EINVAL );
41
42 ihsMaximaStruct *ihsmaxima = NULL;
43 XLAL_CHECK_NULL( ( ihsmaxima = XLALMalloc( sizeof( *ihsmaxima ) ) ) != NULL, XLAL_ENOMEM );
44
45 //The number of ihs maxima = Sum[fbins - (i-1), {i, 2, rows}]
46 // = fbins*rows - fbins - (rows**2 - rows)/2
47 UINT4 numberofmaxima = fbins * rows - fbins - ( UINT4 )( ( rows * rows - rows ) / 2 );
48
49 XLAL_CHECK_NULL( ( ihsmaxima->maxima = XLALCreateREAL4VectorAligned( numberofmaxima, 32 ) ) != NULL, XLAL_EFUNC );
50 XLAL_CHECK_NULL( ( ihsmaxima->locations = XLALCreateINT4Vector( numberofmaxima ) ) != NULL, XLAL_EFUNC );
51 XLAL_CHECK_NULL( ( ihsmaxima->foms = XLALCreateREAL4VectorAligned( numberofmaxima, 32 ) ) != NULL, XLAL_EFUNC );
52 XLAL_CHECK_NULL( ( ihsmaxima->maximaForEachFbin = XLALCreateREAL4VectorAligned( fbins, 32 ) ) != NULL, XLAL_EFUNC );
53 XLAL_CHECK_NULL( ( ihsmaxima->locationsForEachFbin = XLALCreateINT4Vector( fbins ) ) != NULL, XLAL_EFUNC );
54 ihsmaxima->rows = rows;
55
56 return ihsmaxima;
57
58} /* createihsMaxima() */
59
60
61/**
62 * Destroy vectors and the IHS maxima struct
63 * \param [in] data Pointer to an ihsMaximaStruct to be freed
64 */
66{
67 if ( data ) {
69 XLALDestroyINT4Vector( data->locations );
71 XLALDestroyREAL4VectorAligned( data->maximaForEachFbin );
72 XLALDestroyINT4Vector( data->locationsForEachFbin );
74 }
75} /* destroyihsMaxima() */
76
77
78/**
79 * Run the IHS algorithm
80 * \param [out] output Pointer to the ihsMaximaStruct
81 * \param [in] input Pointer to the ffdataStruct
82 * \param [in] ihsfarinput Pointer to the ihsfarStruct
83 * \param [in] params Pointer to UserInput_t
84 * \param [in] rows Number of neighboring rows to be summed
85 * \param [in] aveNoise Pointer to a REAL4VectorAligned of 2nd FFT background powers
86 * \param [in] FbinMean Pointer to a REAL4VectorAligned of normalized SFT background powers
87 * \return Status value
88 */
89INT4 runIHS( ihsMaximaStruct *output, const ffdataStruct *input, const ihsfarStruct *ihsfarinput, const UserInput_t *params, const UINT4 rows, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *FbinMean )
90{
91
92 XLAL_CHECK( output != NULL && input != NULL && ihsfarinput != NULL && params != NULL && rows > 0 && aveNoise != NULL && FbinMean != NULL, XLAL_EINVAL );
93
94 INT4 numfbins = input->numfbins;
95 INT4 numfprbins = input->numfprbins;
96
97 //Allocate memory for the necessary vectors
98 REAL4VectorAligned *row = NULL, *ihss = NULL, *ihsvector = NULL;
99 INT4Vector *locs = NULL;
100 ihsVals *ihsvals = NULL;
101 REAL4VectorAlignedArray *ihsvectorarray = NULL;
102 XLAL_CHECK( ( row = XLALCreateREAL4VectorAligned( numfprbins, 32 ) ) != NULL, XLAL_EFUNC );
103 XLAL_CHECK( ( ihss = XLALCreateREAL4VectorAligned( numfbins, 32 ) ) != NULL, XLAL_EFUNC );
104 XLAL_CHECK( ( ihsvector = XLALCreateREAL4VectorAligned( ( INT4 )floor( ( 1.0 / ( REAL8 )params->ihsfactor ) * numfprbins ) - 5, 32 ) ) != NULL, XLAL_EFUNC );
105 XLAL_CHECK( ( locs = XLALCreateINT4Vector( numfbins ) ) != NULL, XLAL_EFUNC );
106 XLAL_CHECK( ( ihsvals = createihsVals() ) != NULL, XLAL_EFUNC );
107 XLAL_CHECK( ( ihsvectorarray = createREAL4VectorAlignedArray( numfbins, ihsvector->length, 32 ) ) != NULL, XLAL_EFUNC );
108
109 //We want to ignore daily and sidereal harmonics, so mark the values
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;
114 INT4Vector *markedharmonics = NULL;
115 XLAL_CHECK( ( markedharmonics = XLALCreateINT4Vector( row->length ) ) != NULL, XLAL_EFUNC );
116 memset( markedharmonics->data, 0, sizeof( INT4 )*markedharmonics->length );
117 //If the user has specified not to notch the harmonics, then we skip the step to mark the notched values
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;
122 }
123 }
124 }
125
126 //Loop through the rows, 1 frequency at a time
127 for ( UINT4 ii = 0; ii < ihss->length; ii++ ) {
128
129 //For each row, populate it with the data for that frequency bin, excluding harmonics of antenna pattern modulation
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 ) {
132 row->data[jj] = 0.0;
133 }
134
135 //Run the IHS algorithm on the row
136 if ( !params->weightedIHS ) {
137 XLAL_CHECK( incHarmSumVector( ihsvector, row, params->ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
138 } else {
139 XLAL_CHECK( incHarmSumVectorWeighted( ihsvector, row, aveNoise, params->ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
140 }
141
142 //Copy the result into the ihsvector array
143 memcpy( ihsvectorarray->data[ii]->data, ihsvector->data, sizeof( REAL4 )*ihsvector->length );
144
145 } /* for ii < ihss->length */
146
147 //Now do the summing of the IHS values
148 XLAL_CHECK( sumIHSarray( output, ihsfarinput, ihsvectorarray, rows, FbinMean, params ) == XLAL_SUCCESS, XLAL_EFUNC );
149
150 //Destroy stuff
151 destroyREAL4VectorAlignedArray( ihsvectorarray );
155 XLALDestroyINT4Vector( locs );
156 XLALDestroyINT4Vector( markedharmonics );
157 destroyihsVals( ihsvals );
158
159 return XLAL_SUCCESS;
160
161} /* runIHS() */
162
163
164/**
165 * Allocate memory for ihsVals struct
166 * \return Pointer to a newly created ihsVals structure
167 */
169{
170 ihsVals *ihsvals = NULL;
171 XLAL_CHECK_NULL( ( ihsvals = XLALMalloc( sizeof( *ihsvals ) ) ) != NULL, XLAL_ENOMEM );
172 return ihsvals;
173} /* createihsVals() */
174
175
176/**
177 * Destroy ihsVals struct
178 * \param [in] ihsvals Pointer to an ihsVals structure to be freed
179 */
180void destroyihsVals( ihsVals *ihsvals )
181{
182 if ( ihsvals ) {
183 XLALFree( ( ihsVals * )ihsvals );
184 }
185} /* destroyihsVals() */
186
187
188/**
189 * Compute the IHS sum maximum
190 * \param [out] output Pointer to the ihsVals structure
191 * \param [in] input Pointer to a REAL4VectorAligned
192 * \param [in] ihsfactor Number of folds of the 2nd FFT
193 * \return Status value
194 */
195INT4 incHarmSum( ihsVals *output, const REAL4VectorAligned *input, const UINT4 ihsfactor )
196{
197
198 XLAL_CHECK( output != NULL && input != NULL && ihsfactor > 0, XLAL_EINVAL );
199
200 output->ihs = 0.0;
201
202 REAL4VectorAligned *tempvect = NULL;
203 XLAL_CHECK( ( tempvect = XLALCreateREAL4VectorAligned( ( INT4 )floor( ( 1.0 / ( REAL8 )ihsfactor ) * input->length ) - 5, 32 ) ) != NULL, XLAL_EFUNC );
204
205 XLAL_CHECK( incHarmSumVector( tempvect, input, ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
206
207 for ( UINT4 ii = 0; ii < tempvect->length; ii++ ) {
208 if ( tempvect->data[ii] > output->ihs ) {
209 output->ihs = tempvect->data[ii];
210 output->loc = ii + 5;
211 }
212 } /* for ii < tempvect->length */
213
215
216 return XLAL_SUCCESS;
217
218} /* incHarmSum() */
219
220
221/**
222 * Compute the IHS vector -- does not compute the maximum value
223 * \param [out] output Pointer to a REAL4VectorAligned to contain the folded values
224 * \param [in] input Pointer to a REAL4VectorAligned for the input to the IHS
225 * \param [in] ihsfactor Number of folds of the 2nd FFT
226 * \return Status value
227 */
229{
230
231 XLAL_CHECK( output != NULL && input != NULL && ihsfactor > 0, XLAL_EINVAL );
232
233 INT4 highval = ( INT4 )floor( ( 1.0 / ( REAL8 )ihsfactor ) * input->length );
234
235 //From 5 up to the highval
236 for ( INT4 ii = 5; ii < highval; ii++ ) {
237 output->data[ii - 5] = 0.0; //first set the value to zero
238 for ( UINT4 jj = 1; jj <= ihsfactor; jj++ ) {
239 output->data[ii - 5] += input->data[ii * jj];
240 }
241
242 //check that we're not going outside of the allowed limits
243 XLAL_CHECK( ihsfactor * ( ii - 1 ) <= input->length - 1, XLAL_EBADLEN );
244 } /* for ii=5 --> highval */
245
246 return XLAL_SUCCESS;
247
248} /* incHarmSumVector() */
249
250
251/**
252 * Compute the noise weighted IHS vector -- does not compute the maximum value
253 * \param [out] output Pointer to a REAL4VectorAligned to contain the folded values
254 * \param [in] input Pointer to a REAL4VectorAligned for the input to the IHS
255 * \param [in] aveNoise Pointer to a REAL4VectorAligned of 2nd FFT background powers
256 * \param [in] ihsfactor Number of folds of the 2nd FFT
257 * \return Status value
258 */
260{
261
262 XLAL_CHECK( output != NULL && input != NULL && aveNoise != NULL && ihsfactor > 0, XLAL_EINVAL );
263
264 INT4 highval = ( INT4 )floor( ( 1.0 / ( REAL8 )ihsfactor ) * input->length );
265
266 //From 5 up to the highval
267 for ( INT4 ii = 5; ii < highval; ii++ ) {
268 output->data[ii - 5] = 0.0; //first set the value to zero
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;
272 }
273
274 //check that we're not going outside of the allowed limits
275 XLAL_CHECK( ihsfactor * ( ii - 1 ) <= input->length - 1, XLAL_EBADLEN );
276 } /* for ii=5 --> highval */
277
278 return XLAL_SUCCESS;
279
280} // incHarmSumVectorWeighted()
281
282
283/**
284 * Allocate memory for ihsfarStruct struct
285 * \param [in] rows Number of neighbors to sum
286 * \param [in] params Pointer to UserInput_t
287 * \return Pointer to a newly allocated ihsfarStruct
288 */
290{
291
292 XLAL_CHECK_NULL( rows > 0 && params != NULL, XLAL_EINVAL );
293
294 ihsfarStruct *ihsfarstruct = NULL;
295 XLAL_CHECK_NULL( ( ihsfarstruct = XLALMalloc( sizeof( *ihsfarstruct ) ) ) != NULL, XLAL_ENOMEM );
296
297 XLAL_CHECK_NULL( ( ihsfarstruct->ihsfar = XLALCreateREAL4VectorAligned( rows - 1, 32 ) ) != NULL, XLAL_EFUNC );
298 XLAL_CHECK_NULL( ( ihsfarstruct->ihsdistMean = XLALCreateREAL4VectorAligned( rows - 1, 32 ) ) != NULL, XLAL_EFUNC );
299 XLAL_CHECK_NULL( ( ihsfarstruct->ihsdistSigma = XLALCreateREAL4VectorAligned( rows - 1, 32 ) ) != NULL, XLAL_EFUNC );
300 XLAL_CHECK_NULL( ( ihsfarstruct->fomfarthresh = XLALCreateREAL4VectorAligned( rows - 1, 32 ) ) != NULL, XLAL_EFUNC );
301 XLAL_CHECK_NULL( ( ihsfarstruct->ihsfomdistMean = XLALCreateREAL4VectorAligned( rows - 1, 32 ) ) != NULL, XLAL_EFUNC );
302 XLAL_CHECK_NULL( ( ihsfarstruct->ihsfomdistSigma = XLALCreateREAL4VectorAligned( rows - 1, 32 ) ) != NULL, XLAL_EFUNC );
303 XLAL_CHECK_NULL( ( ihsfarstruct->expectedIHSVector = XLALCreateREAL4VectorAligned( ( INT4 )floor( ( 1.0 / ( REAL8 )params->ihsfactor ) * ( ( INT4 )floor( floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 ) * 0.5 ) + 1 ) ) - 5, 32 ) ) != NULL, XLAL_EFUNC );
304
305 memset( ihsfarstruct->expectedIHSVector->data, 0, sizeof( REAL4 )*ihsfarstruct->expectedIHSVector->length );
306
307 return ihsfarstruct;
308
309} // createihsfarStruct()
310
311
312/**
313 * Destroy ihsfarStruct struct
314 * \param [in] ihsfarstruct Pointer to the ihsfarStruct to be destroyed
315 */
317{
318 if ( ihsfarstruct ) {
319 XLALDestroyREAL4VectorAligned( ihsfarstruct->ihsfar );
326 XLALFree( ( ihsfarStruct * )ihsfarstruct );
327 }
328} /* destroyihsfarStruct() */
329
330
331/**
332 * Compute the IHS FAR for a sum of a number of rows
333 * \param [out] output Pointer to the output ihsfarStruct
334 * \param [in] params Pointer to UserInput_t
335 * \param [in] rows Number of neighbors to sum
336 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
337 * \param [in] rng Pointer to GSL random number generator
338 * \return Status value
339 */
340INT4 genIhsFar( ihsfarStruct *output, const UserInput_t *params, const UINT4 rows, const REAL4VectorAligned *aveNoise, const gsl_rng *rng )
341{
342
343 XLAL_CHECK( output != NULL && params != NULL && rows > 0 && aveNoise != NULL, XLAL_EINVAL );
344
345 fprintf( stderr, "Determining IHS FAR values... " );
346 fprintf( LOG, "Determining IHS FAR values... " );
347
348 REAL8 Tobs = params->Tobs;
349
350 UINT4 trials = 5 * rows;
351 if ( trials < 1000 ) {
352 trials = 1000;
353 }
354 if ( trials > 5000 ) {
355 fprintf( stderr, "Warning: number of trials may be insufficient given the number of rows to sum\n" );
356 trials = 5000;
357 }
358
359 //Allocations for IHS values for the number of trials
360 REAL4VectorAligned *ihsvector = NULL, *ihss = NULL, *noise = NULL;
361 REAL4VectorAlignedArray *ihsvectorarray = NULL;
362 INT4Vector *locs = NULL;
363 XLAL_CHECK( ( noise = XLALCreateREAL4VectorAligned( aveNoise->length, 32 ) ) != NULL, XLAL_EFUNC );
364 XLAL_CHECK( ( ihsvector = XLALCreateREAL4VectorAligned( ( INT4 )floor( ( 1.0 / ( REAL8 )params->ihsfactor ) * aveNoise->length ) - 5, 32 ) ) != NULL, XLAL_EFUNC );
365 XLAL_CHECK( ( ihss = XLALCreateREAL4VectorAligned( trials, 32 ) ) != NULL, XLAL_EFUNC );
366 XLAL_CHECK( ( ihsvectorarray = createREAL4VectorAlignedArray( trials, ihsvector->length, 32 ) ) != NULL, XLAL_EFUNC );
367 XLAL_CHECK( ( locs = XLALCreateINT4Vector( trials ) ) != NULL, XLAL_EFUNC );
368
369 //Determine the locations of the harmonics of the earth's rotation in the IHS vector
370 //Amplitude modulations caused by the varying antenna pattern can sometimes cause excess power, so we ignore these harmonics
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;
375 INT4Vector *markedharmonics = NULL;
376 XLAL_CHECK( ( markedharmonics = XLALCreateINT4Vector( aveNoise->length ) ) != NULL, XLAL_EFUNC );
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;
382 }
383 }
384 }
385
386 //Do the expected IHS values from the expected background here
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;
390 }
391 if ( !params->weightedIHS ) {
392 XLAL_CHECK( incHarmSumVector( output->expectedIHSVector, noise, params->ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
393 } else {
394 XLAL_CHECK( incHarmSumVectorWeighted( output->expectedIHSVector, noise, aveNoise, params->ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
395 }
396
397 //Now do a number of trials
398 for ( UINT4 ii = 0; ii < trials; ii++ ) {
399 //Make exponential noise, removing harmonics of 24 hours to match with the same method as real analysis
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] ) );
403 } else {
404 noise->data[jj] = 0.0;
405 }
406 } /* for jj < aveNoise->length */
407
408 //Make a random number 1 +/- 0.2 to create the variations in the nosie that we typically observe
409 //This number needs to be positive
410 REAL8 randval = 1.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 ); //limit range of variation
414 }
415
416 //scale the exponential noise
417 XLAL_CHECK( XLALVectorScaleREAL4( noise->data, randval, noise->data, noise->length ) == XLAL_SUCCESS, XLAL_EFUNC );
418
419 //Compute IHS value on exponential noise
420 if ( !params->weightedIHS ) {
421 XLAL_CHECK( incHarmSumVector( ihsvector, noise, params->ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
422 } else {
423 XLAL_CHECK( incHarmSumVectorWeighted( ihsvector, noise, aveNoise, params->ihsfactor ) == XLAL_SUCCESS, XLAL_EFUNC );
424 }
425
426 //Copy the result into the IHS vector array
427 memcpy( ihsvectorarray->data[ii]->data, ihsvector->data, sizeof( REAL4 )*ihsvector->length );
428 } /* for ii < trials */
429
430 //Destroy stuff
433 XLALDestroyINT4Vector( markedharmonics );
434
435 //Create a fake vector with the same average value in each bin = 1.0
436 REAL4VectorAligned *FbinMean = NULL;
437 XLAL_CHECK( ( FbinMean = XLALCreateREAL4VectorAligned( trials, 32 ) ) != NULL, XLAL_EFUNC );
438 for ( UINT4 ii = 0; ii < trials; ii++ ) {
439 FbinMean->data[ii] = 1.0;
440 }
441
442 //Calculate the IHS sum values for the IHS trials
443 XLAL_CHECK( sumIHSarrayFAR( output, ihsvectorarray, rows, FbinMean, params, rng ) == XLAL_SUCCESS, XLAL_EFUNC );
444
445 //Destroy stuff
446 destroyREAL4VectorAlignedArray( ihsvectorarray );
449 XLALDestroyINT4Vector( locs );
450
451 fprintf( stderr, "done\n" );
452 fprintf( LOG, "done.\n" );
453
454 return XLAL_SUCCESS;
455
456} /* genIhsFar() */
457
458
459/**
460 * Compute the IHS sums for a number of rows used for the FAR calculation
461 * \param [out] outputfar Pointer to the output ihsfarStruct
462 * \param [in] ihsvectorarray Pointer to REAL4VectorAlignedArray to be summed
463 * \param [in] rows Number of neighbors to sum
464 * \param [in] FbinMean Pointer to REAL4VectorAligned of normalized SFT background powers
465 * \param [in] params Pointer to UserInput_t
466 * \param [in] rng Pointer to random number generator
467 * \return Status value
468 */
469INT4 sumIHSarrayFAR( ihsfarStruct *outputfar, REAL4VectorAlignedArray *ihsvectorarray, const UINT4 rows, const REAL4VectorAligned *FbinMean, const UserInput_t *params, const gsl_rng *rng )
470{
471
472 XLAL_CHECK( outputfar != NULL && ihsvectorarray != NULL && rows > 0 && FbinMean != NULL && params != NULL, XLAL_EINVAL );
473
474 //The minimum and maximum index to search in the IHS vector
475 INT4 maxIndexForIHS = ( INT4 )ceil( fmin( params->Tobs / params->Pmin, fmin( params->Tobs / minPeriod( params->dfmin, params->Tsft ), params->Tobs / ( 4.0 * params->Tsft ) ) ) ) - 5;
476 INT4 minIndexForIHS = ( INT4 )floor( fmax( 5.0, params->Tobs / params->Pmax ) ) - 5;
477
478 //Allocate a vector array that holds the summed values of at least two nearest neighbor rows
479 //On the first iteration this holds the nearest neighbor sums, but on subsequent iterations, this holds nearest 3 neighbors
480 //sums, then nearest 4 neighbor sums, and so on.
481 REAL4VectorAlignedArray *tworows = NULL;
482 XLAL_CHECK( ( tworows = createREAL4VectorAlignedArray( ihsvectorarray->length - 1, ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
483 for ( UINT4 ii = 0; ii < tworows->length; ii++ ) {
484 memset( tworows->data[ii]->data, 0, sizeof( REAL4 )*tworows->data[ii]->length ); //Set everything to 0 at the start
485 }
486
487 //Allocate vectors of the ihs values and locations of the maximums
488 //Vectors for values above the noise and scaling the noise
489 REAL4VectorAligned *ihsvalues = NULL;
490 INT4Vector *ihslocations = NULL;
491 XLAL_CHECK( ( ihsvalues = XLALCreateREAL4VectorAligned( ihsvectorarray->length, 32 ) ) != NULL, XLAL_EFUNC );
492 XLAL_CHECK( ( ihslocations = XLALCreateINT4Vector( ihsvectorarray->length ) ) != NULL, XLAL_EFUNC );
493 REAL4VectorAligned *excessabovenoise = NULL, *scaledExpectedIHSVectorValues = NULL;
494 XLAL_CHECK( ( excessabovenoise = XLALCreateREAL4VectorAligned( ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
495 XLAL_CHECK( ( scaledExpectedIHSVectorValues = XLALCreateREAL4VectorAligned( ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
496
497 //Finding the maximum for each IHS vector and the location
498 for ( UINT4 ii = 0; ii < ihsvalues->length; ii++ ) {
499 //Scale the expected IHS vector by the value in FbinMean (in this function, it is 1.0)
500 XLAL_CHECK( XLALVectorScaleREAL4( scaledExpectedIHSVectorValues->data, FbinMean->data[ii], outputfar->expectedIHSVector->data, outputfar->expectedIHSVector->length ) == XLAL_SUCCESS, XLAL_EFUNC );
501 //subtract the noise from the data
502 XLAL_CHECK( VectorSubtractREAL4( excessabovenoise, ihsvectorarray->data[ii], scaledExpectedIHSVectorValues, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
503 //Find the location of the maximum IHS value
504 ihslocations->data[ii] = max_index_in_range( excessabovenoise, minIndexForIHS, maxIndexForIHS ) + 5;
505 //And get the value
506 ihsvalues->data[ii] = ihsvectorarray->data[ii]->data[ihslocations->data[ii] - 5];
507 //fprintf(stderr, "%d %f\n", ihslocations->data[ii], ihsvalues->data[ii]);
508 } /* for ii=0 --> ihsvalues->length */
509
510 //Some useful variables
511 INT4Vector *rowarraylocs = NULL;
512 REAL4VectorAligned *foms = NULL;
513
514 //Starting from a minimum of 2 rows, start determining the FAR for each nearest neighbor sum, up to the maximum number of
515 //rows to be summed
516 for ( UINT4 ii = 2; ii <= rows; ii++ ) {
517 //First allocate the necessary vectors
518 XLAL_CHECK( ( rowarraylocs = XLALCreateINT4Vector( ii ) ) != NULL, XLAL_EFUNC );
519 XLAL_CHECK( ( foms = XLALCreateREAL4VectorAligned( ihsvectorarray->length - ( ii - 1 ), 32 ) ) != NULL, XLAL_EFUNC );
520
521 //If the user has specified that we should use SSE operations, then do the nearest neighbor summing.
522 //The result is in the tworows variable
523 if ( params->vectorMath != 0 ) {
524 if ( ii > 2 ) {
525 XLAL_CHECK( VectorArraySum( tworows, tworows, ihsvectorarray, 0, ii - 1, 0, ( INT4 )( ihsvectorarray->length - ( ii - 1 ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
526 } else {
527 XLAL_CHECK( VectorArraySum( tworows, ihsvectorarray, ihsvectorarray, 0, ii - 1, 0, ( INT4 )( ihsvectorarray->length - ( ii - 1 ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
528 }
529 }
530
531 //Now we are going to loop through the input ihsvectorarray up to the number of rows-1
532 for ( UINT4 jj = 0; jj < ihsvectorarray->length - ( ii - 1 ); jj++ ) {
533 //A bit tricky here: sum IHS values across SFT frequency bins into the tworows variable only if we didn't do it with SSE above
534 if ( params->vectorMath == 0 ) {
535 if ( ii > 2 ) for ( UINT4 kk = 0; kk < tworows->data[jj]->length; kk++ ) {
536 tworows->data[jj]->data[kk] += ihsvectorarray->data[ii - 1 + jj]->data[kk];
537 } else for ( UINT4 kk = 0; kk < tworows->data[jj]->length; kk++ ) {
538 tworows->data[jj]->data[kk] = ihsvectorarray->data[jj]->data[kk] + ihsvectorarray->data[jj + 1]->data[kk];
539 }
540 }
541
542 //Compute IHS FOM value
543 memcpy( rowarraylocs->data, &( ihslocations->data[jj] ), sizeof( INT4 )*ii ); //First copy the necessary values
544 foms->data[jj] = ihsFOM( rowarraylocs, ( INT4 )outputfar->expectedIHSVector->length ); //then compute the FOM
546 } /* for jj< ihsvectorarray->length - (ii-1) */
547
548 //Sample the IHS values that have been summed to compute mean, standard deviation, and FAR threshold values.
549 //We have an if-else statement for when there are fewer than 10000 entries that will be in the tworows varaible
550 //(only considering the number of rows we have summed together).
551 REAL4VectorAligned *sampledtempihsvals = NULL;
552 if ( ( ihsvectorarray->length - ( ii - 1 ) ) * ihsvectorarray->data[0]->length > 10000 ) {
553 //We sample the tworows array (up to the number of rows-1) without accepting any zeros
554 //because zeros would come from the notched harmonics, which we won't want anyway
555 sampledtempihsvals = sampleREAL4VectorAlignedArray_nozerosaccepted( tworows, ihsvectorarray->length - ( ii - 1 ), 10000, rng );
556 outputfar->ihsdistMean->data[ii - 2] = calcMean( sampledtempihsvals ); //And then calculate the mean value
557 XLAL_CHECK( calcStddev( &( outputfar->ihsdistSigma->data[ii - 2] ), sampledtempihsvals ) == XLAL_SUCCESS, XLAL_EFUNC ); //We also calculate the standard deviation
558 } else {
559 //If there were fewer than 10000 entries, then we will keep all those that are part of the nearest neighbor
560 //sum up to this point
561 XLAL_CHECK( ( sampledtempihsvals = XLALCreateREAL4VectorAligned( ( ihsvectorarray->length - ( ii - 1 ) ) * ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
562 for ( UINT4 jj = 0; jj < ( ihsvectorarray->length - ( ii - 1 ) ); jj++ ) {
563 memcpy( sampledtempihsvals->data, tworows->data[jj]->data, sizeof( REAL4 )*tworows->data[0]->length );
564 }
565 outputfar->ihsdistMean->data[ii - 2] = calcMean_ignoreZeros( sampledtempihsvals ); //Calculate the mean value (remember, don't accept zeros)
566 XLAL_CHECK( calcStddev_ignoreZeros( &( outputfar->ihsdistSigma->data[ii - 2] ), sampledtempihsvals ) == XLAL_SUCCESS, XLAL_EFUNC ); //We also calculate the standard deviation
567 }
568
569 //If the user has specified the IHS FAR == 1.0, then we don't need to compute the threshold (it is = 0.0)
570 INT4 numvals = 0;
571 REAL8 farave = 0.0;
572 if ( params->ihsfar != 1.0 ) {
573 //Looping through the sampled values, we are going to compute the average FAR
574 for ( UINT4 jj = 0; jj < sampledtempihsvals->length; jj++ ) {
575 //When the user has not specified using faster chisq inversion, use the GSL function
576 if ( !params->fastchisqinv && sampledtempihsvals->data[jj] != 0.0 ) {
577 numvals++;
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 ) {
580 numvals++;
581 farave += 0.5 * cdf_chisq_Qinv( params->ihsfar, 2.0 * sampledtempihsvals->data[jj] );
583 } //fastchisqinv?
584 } //for jj=0 --> sampledtempihsval->length
585 outputfar->ihsfar->data[ii - 2] = farave / ( REAL8 )numvals;
586 } //if params->ihsfar != 1.0
587 else {
588 outputfar->ihsfar->data[ii - 2] = 0.0;
589 }
590
591 XLALDestroyREAL4VectorAligned( sampledtempihsvals );
592
593 //FOM part
594 outputfar->ihsfomdistMean->data[ii - 2] = calcMean( foms );
595 XLAL_CHECK( calcStddev( &( outputfar->ihsfomdistSigma->data[ii - 2] ), foms ) == XLAL_SUCCESS, XLAL_EFUNC );
596 //We need to find the smallest values
597 if ( params->ihsfomfar != 1.0 && params->ihsfom == 0.0 ) {
598 REAL4VectorAligned *smallestfomvals = NULL;
599 XLAL_CHECK( ( smallestfomvals = XLALCreateREAL4VectorAligned( ( INT4 )round( ( ihsvalues->length - ii + 1 ) * params->ihsfomfar ) + 1, 32 ) ) != NULL, XLAL_EFUNC );
600 XLAL_CHECK( sort_float_smallest( smallestfomvals, foms ) == XLAL_SUCCESS, XLAL_EFUNC ); //Sort, finding the smallest values
601 outputfar->fomfarthresh->data[ii - 2] = smallestfomvals->data[smallestfomvals->length - 1]; //Pick off the last element
602 XLALDestroyREAL4VectorAligned( smallestfomvals );
603 } else if ( params->ihsfom != 0.0 ) {
604 outputfar->fomfarthresh->data[ii - 2] = params->ihsfom;
605 } else {
606 outputfar->fomfarthresh->data[ii - 2] = -1.0;
607 }
608
609 XLALDestroyINT4Vector( rowarraylocs );
610 rowarraylocs = NULL;
612 foms = NULL;
613 } /* for ii <= rows */
614
617 XLALDestroyREAL4VectorAligned( excessabovenoise );
618 XLALDestroyREAL4VectorAligned( scaledExpectedIHSVectorValues );
619 XLALDestroyINT4Vector( ihslocations );
620
621 return XLAL_SUCCESS;
622
623} /*sumIHSarrayFAR() */
624
625
626/**
627 * \brief Compute the IHS sums for a number of rows
628 *
629 * In the function we will select the the location which is the maximum above the noise
630 *
631 * \param [out] output Pointer to the output ihsMaximaStruct
632 * \param [in] inputfar Pointer to ihsfarStruct
633 * \param [in] ihsvectorarray Pointer to REAL4VectorAlignedArray to be summed
634 * \param [in] rows Number of neighbors to sum
635 * \param [in] FbinMean Pointer to REAL4VectorAligned of normalized SFT background powers
636 * \param [in] params Pointer to UserInput_t
637 * \return Status value
638 */
639INT4 sumIHSarray( ihsMaximaStruct *output, const ihsfarStruct *inputfar, REAL4VectorAlignedArray *ihsvectorarray, const UINT4 rows, const REAL4VectorAligned *FbinMean, const UserInput_t *params )
640{
641
642 XLAL_CHECK( output != NULL && inputfar != NULL && ihsvectorarray != NULL && rows > 0 && FbinMean != NULL && params != NULL, XLAL_EINVAL );
643
644 //Again, we start off by allocating a "tworows" vector array of IHS nearest neighbor sums
645 REAL4VectorAlignedArray *tworows = NULL;
646 XLAL_CHECK( ( tworows = createREAL4VectorAlignedArray( ihsvectorarray->length - 1, ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
647 for ( UINT4 ii = 0; ii < tworows->length; ii++ ) {
648 memset( tworows->data[ii]->data, 0, sizeof( REAL4 )*tworows->data[0]->length ); //Set everything to 0.0
649 }
650
651 //Allocation of ihs values and locations
652 //Vectors for values above the noise and scaling the noise
653 REAL4VectorAligned *ihsvalues = NULL, *excessabovenoise = NULL, *scaledExpectedIHSVectorValues = NULL;
654 INT4Vector *ihslocations = NULL;
655 XLAL_CHECK( ( ihsvalues = XLALCreateREAL4VectorAligned( ihsvectorarray->length, 32 ) ) != NULL, XLAL_EFUNC );
656 XLAL_CHECK( ( excessabovenoise = XLALCreateREAL4VectorAligned( ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
657 XLAL_CHECK( ( scaledExpectedIHSVectorValues = XLALCreateREAL4VectorAligned( ihsvectorarray->data[0]->length, 32 ) ) != NULL, XLAL_EFUNC );
658 XLAL_CHECK( ( ihslocations = XLALCreateINT4Vector( ihsvectorarray->length ) ) != NULL, XLAL_EFUNC );
659
660 //The minimum and maximum index to search in the IHS vector
661 INT4 maxIndexForIHS = ( INT4 )ceil( fmin( params->Tobs / params->Pmin, fmin( params->Tobs / minPeriod( params->dfmin, params->Tsft ), params->Tobs / ( 4.0 * params->Tsft ) ) ) ) - 5;
662 INT4 minIndexForIHS = ( INT4 )floor( fmax( 5.0, params->Tobs / params->Pmax ) ) - 5;
663
664 //Finding the maximum for each IHS vector and the location using SSE or not
665 for ( UINT4 ii = 0; ii < ihsvalues->length; ii++ ) {
666 //Scale the expected IHS vector by the FbinMean data value
667 XLAL_CHECK( XLALVectorScaleREAL4( scaledExpectedIHSVectorValues->data, FbinMean->data[ii], inputfar->expectedIHSVector->data, inputfar->expectedIHSVector->length ) == XLAL_SUCCESS, XLAL_EFUNC );
668 //subtract the noise from the data
669 XLAL_CHECK( VectorSubtractREAL4( excessabovenoise, ihsvectorarray->data[ii], scaledExpectedIHSVectorValues, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
670
671 //search over the range of Pmin-->Pmax and higher harmonics the user has specified
672 for ( INT4 jj = 0; jj < params->harmonicNumToSearch; jj++ ) {
673 if ( jj == 0 ) {
674 ihslocations->data[ii] = max_index_in_range( excessabovenoise, minIndexForIHS, maxIndexForIHS ) + 5;
675 ihsvalues->data[ii] = ihsvectorarray->data[ii]->data[ihslocations->data[ii] - 5];
676 } else {
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;
682 } /* if the new value is better than the previous value */
683 }
684 } /* for jj=0 --> jj<harmonicNumToSearch */
685 }
686
687 //Useful variables
688 INT4Vector *rowarraylocs = NULL;
689
690 //Start with the single IHS vector and march up with nearest neighbor sums up to the total number of row sums
691 for ( UINT4 ii = 1; ii <= rows; ii++ ) {
692 if ( ii == 1 ) {
693 //Copy the data into the output
694 memcpy( output->maximaForEachFbin->data, ihsvalues->data, sizeof( REAL4 )*ihsvalues->length );
695 memcpy( output->locationsForEachFbin->data, ihslocations->data, sizeof( INT4 )*ihslocations->length );
696 } else {
697 //For everything 2 nearest neighbors and higher summed
698 XLAL_CHECK( ( rowarraylocs = XLALCreateINT4Vector( ii ) ) != NULL, XLAL_EFUNC );
699
700 //The maximum index to search in the IHS vector
701 maxIndexForIHS = ( INT4 )ceil( fmin( params->Tobs / params->Pmin, fmin( params->Tobs / minPeriod( 0.5 * ( ii - 1 ) / params->Tsft, params->Tsft ), params->Tobs / ( 4.0 * params->Tsft ) ) ) ) - 5;
702
703 REAL4 sumofnoise = 0.0; //To scale the expected IHS background
704 INT4 endloc = ( ( ii - 1 ) * ( ii - 1 ) - ( ii - 1 ) ) / 2;
705
706 //Sum up the IHS vectors using SSE functions
707 if ( params->vectorMath != 0 ) {
708 if ( ii > 2 ) {
709 XLAL_CHECK( VectorArraySum( tworows, tworows, ihsvectorarray, 0, ii - 1, 0, ( INT4 )( ihsvectorarray->length - ( ii - 1 ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
710 } else {
711 XLAL_CHECK( VectorArraySum( tworows, ihsvectorarray, ihsvectorarray, 0, ii - 1, 0, ( INT4 )( ihsvectorarray->length - ( ii - 1 ) ) ) == XLAL_SUCCESS, XLAL_EFUNC );
712 }
713 }
714
715 //Loop through the IHS vector neighbor sums
716 for ( UINT4 jj = 0; jj < ihsvectorarray->length - ( ii - 1 ); jj++ ) {
717 //If we didn't use SSE to sum the vector array (see lines above)
718 if ( params->vectorMath == 0 ) {
719 if ( ii > 2 ) for ( UINT4 kk = 0; kk < tworows->data[jj]->length; kk++ ) {
720 tworows->data[jj]->data[kk] += ihsvectorarray->data[ii - 1 + jj]->data[kk];
721 } else for ( UINT4 kk = 0; kk < tworows->data[jj]->length; kk++ ) {
722 tworows->data[jj]->data[kk] = ihsvectorarray->data[jj]->data[kk] + ihsvectorarray->data[jj + 1]->data[kk];
723 }
724 }
725
726 //To scale the background efficiently
727 if ( jj == 0 ) for ( UINT4 kk = 0; kk < ii; kk++ ) {
728 sumofnoise += FbinMean->data[kk];
729 } else {
730 sumofnoise -= FbinMean->data[jj - 1];
731 sumofnoise += FbinMean->data[jj + ( ii - 1 )];
732 }
733
734 //If using SSE, scale the expected IHS vector, subtract the noise from the data
735 XLAL_CHECK( XLALVectorScaleREAL4( scaledExpectedIHSVectorValues->data, sumofnoise, inputfar->expectedIHSVector->data, inputfar->expectedIHSVector->length ) == XLAL_SUCCESS, XLAL_EFUNC );
736 XLAL_CHECK( VectorSubtractREAL4( excessabovenoise, tworows->data[jj], scaledExpectedIHSVectorValues, params->vectorMath ) == XLAL_SUCCESS, XLAL_EFUNC );
737
738 //Compute the maximum IHS value in the second FFT frequency direction
739 //search over the range of Pmin-->Pmax and higher harmonics the user has specified
740 for ( INT4 kk = 0; kk < params->harmonicNumToSearch; kk++ ) {
741 if ( kk == 0 ) {
742 output->locations->data[( ii - 2 )*ihsvalues->length - endloc + jj] = max_index_in_range( excessabovenoise, minIndexForIHS, maxIndexForIHS ) + 5;
743 output->maxima->data[( ii - 2 )*ihsvalues->length - endloc + jj] = tworows->data[jj]->data[( output->locations->data[( ii - 2 ) * ihsvalues->length - endloc + jj] - 5 )];
744 } else {
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;
750 } /* if the new value is better than the previous value */
751 }
752 } /* for kk=0 --> kk<harmonicNumToSearch */
753
754 memcpy( rowarraylocs->data, &( ihslocations->data[jj] ), sizeof( INT4 )*ii );
755 output->foms->data[( ii - 2 )*ihsvalues->length - endloc + jj] = ihsFOM( rowarraylocs, ( INT4 )inputfar->expectedIHSVector->length );
757 } /* for jj< ihsvectorarray->length - (ii-1) */
758
759 XLALDestroyINT4Vector( rowarraylocs );
760 rowarraylocs = NULL;
761 }
762
763 } /* for ii <= rows */
764
766 XLALDestroyREAL4VectorAligned( scaledExpectedIHSVectorValues );
767 XLALDestroyREAL4VectorAligned( excessabovenoise );
769 XLALDestroyINT4Vector( ihslocations );
770
771 return XLAL_SUCCESS;
772
773} /*sumIHSarray() */
774
775
776/**
777 * Calculate the IHS FOM for a number of rows
778 * \param [in] locs Pointer to INT4Vector of location values
779 * \param [in] fomnorm Normalization value and is the number of XXX
780 * \return IHS FOM value
781 */
782REAL4 ihsFOM( const INT4Vector *locs, const INT4 fomnorm )
783{
784
785 XLAL_CHECK_REAL4( locs != NULL && fomnorm > 0.0, XLAL_EINVAL );
786
787 REAL4 fom = 0.0;
788
789 for ( INT4 ii = 0; ii < ( INT4 )( locs->length * 0.5 ); ii++ ) {
790 fom += ( REAL4 )( ( locs->data[ii] - locs->data[locs->length - ii - 1] ) * ( locs->data[ii] - locs->data[locs->length - ii - 1] ) ) / ( fomnorm * fomnorm );
791 }
792
793 return fom;
794
795} /* ihsFOM() */
796
797
798/**
799 * Find IHS candidates above thresholds
800 * \param [out] candlist Pointer to a pointer containing the candidate list
801 * \param [in] ihsfarstruct Pointer to ihsfarStruct
802 * \param [in] params Pointer to UserInput_t
803 * \param [in] ffdata Pointer to ffdataStruct
804 * \param [in] ihsmaxima Pointer to ihsMaximaStruct containing the data to be tested above thresholds
805 * \param [in] fbinavgs Pointer to REAL4VectorAligned of normalized SFT background powers
806 * \param [in] trackedlines Pointer to REAL4VectorSequence of lines (allowed to be NULL if no lines)
807 * \return Status value
808 */
809INT4 findIHScandidates( candidateVector **candlist, const ihsfarStruct *ihsfarstruct, const UserInput_t *params, const ffdataStruct *ffdata, const ihsMaximaStruct *ihsmaxima, const REAL4VectorAligned *fbinavgs, const REAL4VectorSequence *trackedlines )
810{
811
812 XLAL_CHECK( ihsfarstruct != NULL && params != NULL && ffdata != NULL && ihsmaxima != NULL && fbinavgs != NULL, XLAL_EINVAL );
813
814 REAL8 fsig, per0, B;
815 INT4 numberofIHSvalsChecked = 0, numberofIHSvalsExceededThresh = 0, numberPassingBoth = 0, linesinterferewithnum = 0, skipped = 0, notskipped = 0;
816
817 INT4 numfbins = ffdata->numfbins; //number of fbins
818 INT4 minrows = ( INT4 )round( 2.0 * params->dfmin * params->Tsft ) + 1; //minimum number of sequential rows to search
819
820 REAL4VectorAligned *ihss = NULL, *avgsinrange = NULL;
821 INT4Vector *locs = NULL;
822
823 //Check the IHS values against the FAR, checking between IHS width values
824 //FILE *IHSVALSOUTPUT = fopen("./output/allihsvalspassthresh.dat","w");
825 for ( UINT4 ii = minrows; ii <= ihsfarstruct->ihsfar->length + 1; ii++ ) {
826 XLAL_CHECK( ( ihss = XLALCreateREAL4VectorAligned( ii, 32 ) ) != NULL, XLAL_EFUNC );
827 XLAL_CHECK( ( locs = XLALCreateINT4Vector( ii ) ) != NULL, XLAL_EFUNC );
828 XLAL_CHECK( ( avgsinrange = XLALCreateREAL4VectorAligned( ii, 32 ) ) != NULL, XLAL_EFUNC );
829
830 REAL8 highestval = 0.0, highestsignificance = 0.0; //highestvalnoise = 0.0
831 INT4 highestvalloc = -1, jjloc = 0;
832 for ( UINT4 jj = 0; jj < numfbins - ( ii - 1 ); jj++ ) {
833
834 //Noise in the range of the rows, mean for IHS
835 memcpy( avgsinrange->data, &( fbinavgs->data[jj] ), sizeof( REAL4 )*ii );
836 REAL4 meanNoise = calcMean( avgsinrange );
837
838 numberofIHSvalsChecked++; //count the number of ihs values checked
839
840 INT4 locationinmaximastruct = ( ii - 2 ) * numfbins - ( ( ii - 1 ) * ( ii - 1 ) - ( ii - 1 ) ) / 2 + jj; //the location in the ihsmaximastruct
841
842 //Check the IHS sum against the FAR (scaling FAR with mean of the noise in the range of rows)
843 if ( ihsmaxima->maxima->data[locationinmaximastruct] > ihsfarstruct->ihsfar->data[ii - 2]*meanNoise ) {
844
845 numberofIHSvalsExceededThresh++; //count the number of values exceeding the IHS value threshold
846
847 if ( ihsfarstruct->fomfarthresh->data[ii - 2] == -1.0 || ihsmaxima->foms->data[locationinmaximastruct] <= ihsfarstruct->fomfarthresh->data[ii - 2] ) {
848
849 numberPassingBoth++; //count number passing both IHS value and FOM value thresholds
850
851 INT4 loc = ihsmaxima->locations->data[locationinmaximastruct];
852 per0 = params->Tobs / loc; //Candidate period
853 fsig = params->fmin - params->dfmax + ( ( 0.5 * ( ii - 1 ) + jj ) - 6.0 ) / params->Tsft; //Candidate frequency
854 B = 0.5 * ( ii - 1 ) / params->Tsft; //Candidate modulation depth
855
856 //Test to see if any tracked lines are overlapping the candidate signal
857 INT4 nolinesinterfering = 1;
858 if ( trackedlines != NULL ) {
859 UINT4 kk = 0;
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;
865 }
866 } // if the band spanned by the line is smaller than the band spanned by the signal
867 else {
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;
871 }
872 } // instead if the band spanned by the line is larger than the band spanned by the signal
873 kk++;
874 } // while kk < trackedlines->length && nolinesinterfering==1
875 } // if trackedlines != NULL
876
877 if ( !nolinesinterfering ) {
878 linesinterferewithnum++;
879 } else {
880 REAL8 noise = ihsfarstruct->ihsdistMean->data[ii - 2];
881 REAL8 totalnoise = meanNoise * noise;
882
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 );
886 } else {
887 significance = -log10( significance );
888 }
889
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;
894 jjloc = jj;
895 notskipped++;
896 } else {
897 skipped++;
898 }
899
900 } // if no lines are interfering
901 } // if fom is below or equal to threshold fom
902 } // if val exceeds threshold
903 } // for jj < numfbins-(ii-1)
904
905 if ( highestvalloc != -1 ) {
906 INT4 loc = ihsmaxima->locations->data[highestvalloc];
907 fsig = params->fmin - params->dfmax + ( 0.5 * ( ii - 1 ) + jjloc - 6.0 ) / params->Tsft; //Candidate frequency
908 B = 0.5 * ( ii - 1 ) / params->Tsft; //Candidate modulation depth
909 per0 = params->Tobs / loc; //Candidate period
910 REAL8 h0 = ihs2h0( 2.0 * highestval, params ); //Candidate h0, need factor of 2 for the degrees of freedom counting
911
912 if ( ( *candlist )->numofcandidates == ( *candlist )->length - 1 ) {
913 XLAL_CHECK( ( *candlist = resizecandidateVector( *candlist, 2 * ( *candlist )->length ) ) != NULL, XLAL_EFUNC );
914 }
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++;
917 }
918
919 //Destroy
921 XLALDestroyREAL4VectorAligned( avgsinrange );
922 XLALDestroyINT4Vector( locs );
923 ihss = NULL;
924 locs = NULL;
925 avgsinrange = NULL;
926
927 } /* for ii < ihsfarstruct->ihsfar->length */
928
929 //The outer loop is over "frequency" (the order in the ihs maxima vector, first 2 row sums, then three, and so on)
930 //for (ii=0; ii<(INT4)numfbins; ii++) {
931 /* for (ii=6+(INT4)round(params->dfmax*params->Tsft); ii<(numfbins-6-(INT4)round(params->dfmax*params->Tsft)); ii++) {
932 REAL8 highestval = 0.0, highestsignificance = 0.0; //highestvalnoise = 0.0
933 INT4 highestvalloc = -1, jjloc = 0;
934
935 //This controls the maximum modulation depth to search which depends on the position in the "frequency" loop
936 //INT4 maxrows = numfbins-ii;
937 INT4 jjmax = (INT4)ihsfarstruct->ihsfar->length+1;
938 //if (maxrows<(INT4)ihsfarstruct->ihsfar->length+1) jjmax = maxrows;
939 //else jjmax = (INT4)ihsfarstruct->ihsfar->length+1;
940
941 //Inner loop over modulation depth
942 for (jj=minrows; jj<=jjmax; jj++) {
943 numberofIHSvalsChecked++;
944
945 ihss = XLALCreateREAL4Vector(jj);
946 locs = XLALCreateINT4Vector(jj);
947 avgsinrange = XLALCreateREAL4Vector(jj);
948 if (ihss==NULL) {
949 fprintf(stderr,"%s: XLALCreateREAL4Vector(%d) failed.\n", __func__, jj);
950 XLAL_ERROR_VOID(XLAL_EFUNC);
951 } else if (locs==NULL) {
952 fprintf(stderr,"%s: XLALCreateINT4Vector(%d) failed.\n", __func__, jj);
953 XLAL_ERROR_VOID(XLAL_EFUNC);
954 } else if (avgsinrange==NULL) {
955 fprintf(stderr,"%s: XLALCreateREAL4Vector(%d) failed.\n", __func__, jj);
956 XLAL_ERROR_VOID(XLAL_EFUNC);
957 }
958
959 //Noise in the range of the rows, mean for IHS
960 memcpy(avgsinrange->data, &(fbinavgs->data[ii]), sizeof(REAL4)*jj);
961 REAL4 meanNoise = calcMean(avgsinrange);
962 if (XLAL_IS_REAL4_FAIL_NAN(meanNoise)) {
963 fprintf(stderr,"%s: calcMean() failed.\n", __func__);
964 XLAL_ERROR_VOID(XLAL_EFUNC);
965 }
966
967 INT4 locationinmaximastruct = (jj-2)*numfbins-((jj-1)*(jj-1)-(jj-1))/2 + ii;
968
969 //Check the IHS sum against the FAR (scaling FAR with mean of the noise in the range of rows)
970 if (ihsmaxima->maxima->data[locationinmaximastruct] > ihsfarstruct->ihsfar->data[jj-2]*meanNoise) {
971 numberofIHSvalsExceededThresh++;
972 if (ihsfarstruct->fomfarthresh->data[jj-2]==-1.0 || ihsmaxima->foms->data[locationinmaximastruct]<=ihsfarstruct->fomfarthresh->data[jj-2]) {
973
974 numberPassingBoth++;
975
976 INT4 loc = ihsmaxima->locations->data[locationinmaximastruct];
977 per0 = params->Tobs/loc; //Candidate period
978 fsig = params->fmin - params->dfmax + (0.5*(jj-1) + ii - 6)/params->Tsft; //Candidate frequency
979 B = 0.5*(jj-1)/params->Tsft; //Candidate modulation depth
980
981 //Test to see if any tracked lines are overlapping the candidate signal
982 INT4 nolinesinterfering = 1;
983 if (trackedlines!=NULL) {
984 kk = 0;
985 while (kk<(INT4)trackedlines->length && nolinesinterfering==1) {
986 if (2.0*B>=(trackedlines->data[kk*3+2]-trackedlines->data[kk*3+1])) {
987 if ((trackedlines->data[kk*3+2]>=fsig-B && trackedlines->data[kk*3+2]<=fsig+B) ||
988 (trackedlines->data[kk*3+1]>=fsig-B && trackedlines->data[kk*3+1]<=fsig+B)) {
989 nolinesinterfering = 0;
990 }
991 } // if the band spanned by the line is smaller than the band spanned by the signal
992 else {
993 if ((fsig+B>=trackedlines->data[kk*3+1] && fsig+B<=trackedlines->data[kk*3+2]) ||
994 (fsig-B>=trackedlines->data[kk*3+1] && fsig-B<=trackedlines->data[kk*3+2])) {
995 nolinesinterfering = 0;
996 }
997 } // instead if the band spanned by the line is larger than the band spanned by the signal
998 kk++;
999 } // while kk < trackedlines->length && nolinesinterfering==1
1000 } // if trackedlines != NULL
1001
1002 if (!nolinesinterfering) {
1003 linesinterferewithnum++;
1004 } else {
1005 REAL8 noise = ihsfarstruct->ihsdistMean->data[jj-2];
1006 REAL8 totalnoise = meanNoise*noise;
1007 REAL8 sigma = calcRms(avgsinrange)*ihsfarstruct->ihsdistSigma->data[jj-2];
1008
1009 REAL8 significance = (2.0*ihsmaxima->maxima->data[locationinmaximastruct] - 2.0*totalnoise)/sqrt(2.0*2.0*sigma);
1010
1011 //if (ihsmaxima->maxima->data[locationinmaximastruct]-totalnoise > highestval) {
1012 //if ( ihsmaxima->maxima->data[locationinmaximastruct]-totalnoise > highestval &&
1013 //(params->followUpOutsideULrange || (!params->followUpOutsideULrange &&
1014 //fsig>=params->ULfmin && fsig<=(params->ULfmin+params->ULfspan) &&
1015 //B>=params->ULmindf && B<=params->ULmaxdf)) ) {
1016 if ( significance > highestsignificance &&
1017 (params->followUpOutsideULrange || (!params->followUpOutsideULrange &&
1018 fsig>=params->ULfmin && fsig<=(params->ULfmin+params->ULfspan) &&
1019 B>=params->ULmindf && B<=params->ULmaxdf)) ) {
1020 //highestval = ihsmaxima->maxima->data[locationinmaximastruct];
1021 highestval = ihsmaxima->maxima->data[locationinmaximastruct]-totalnoise;
1022 //highestvalnoise = totalnoise;
1023 highestvalloc = locationinmaximastruct;
1024 highestsignificance = significance;
1025 jjloc = jj;
1026 notskipped++;
1027 } else {
1028 skipped++;
1029 }
1030 } // if no lines are interfering
1031 } //If exceeding the FOM threshold
1032 } //If exceeding the IHS FAR threshold
1033
1034 XLALDestroyREAL4Vector(ihss);
1035 XLALDestroyREAL4Vector(avgsinrange);
1036 XLALDestroyINT4Vector(locs);
1037 } //loop over modulation depths
1038
1039 if (highestvalloc != -1) {
1040 INT4 loc = ihsmaxima->locations->data[highestvalloc];
1041 //Candidate frequency
1042 fsig = params->fmin - params->dfmax + (0.5*(jjloc-1) + ii - 6.0)/params->Tsft;
1043 //Candidate modulation depth
1044 B = 0.5*(jjloc-1)/params->Tsft;
1045 //Candidate period
1046 per0 = params->Tobs/loc;
1047 //Candidate h0
1048 REAL8 h0 = ihs2h0(2.0*highestval, params); //Need factor of 2 for the degrees of freedom counting
1049 REAL8 significance = highestsignificance;
1050 //fprintf(stderr, "%d %d %f\n", ii, jjloc, significance); //remove this
1051
1052 if (candlist->numofcandidates == candlist->length-1) {
1053 candlist = resizecandidateVector(candlist, 2*(candlist->length));
1054 if (candlist->data==NULL) {
1055 fprintf(stderr,"%s: resizecandidateVector() failed.\n", __func__);
1056 XLAL_ERROR_VOID(XLAL_EFUNC);
1057 }
1058 }
1059 //loadCandidateData(&candlist->data[candlist->numofcandidates], fsig, per0, B, 0.0, 0.0, ihsmaxima->maxima->data[locationinmaximastruct], h0, 0.0, 0, sqrt(ffdata->tfnormalization/2.0*params->Tsft));
1060 loadCandidateData(&candlist->data[candlist->numofcandidates], fsig, per0, B, 0.0, 0.0, ihsmaxima->maxima->data[highestvalloc], h0, significance, 0, ffdata->tfnormalization);
1061 (candlist->numofcandidates)++;
1062 }
1063 } //loop over "frequency" */
1064
1065 //fclose(IHSVALSOUTPUT);
1066
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 );
1071
1072 return XLAL_SUCCESS;
1073
1074} /* findIHScandidates() */
1075
1076
1077/**
1078 * Convert the IHS statistic to an estimated h0, based on injections
1079 * \param [in] ihsval The IHS statistic value
1080 * \param [in] params Pointer to UserInput_t
1081 * \return Estimated h0
1082 */
1083REAL8 ihs2h0( const REAL8 ihsval, const UserInput_t *params )
1084{
1085
1086 if ( ihsval <= 0.0 ) {
1087 return 0.0;
1088 }
1089 REAL8 prefact = 1.0;
1090 //prefact = 7.2; //old value for when IHS --> h0 when in the upper limits, the value of the noise subtracted was incorrectly calculated and for >99.999% CL
1091 prefact = 44.7; //new value for IHS --> h0
1092 return prefact * pow( ihsval / ( params->Tsft * params->Tobs ), 0.25 );
1093
1094}
void destroyihsfarStruct(ihsfarStruct *ihsfarstruct)
Destroy ihsfarStruct struct.
Definition: IHS.c:316
void destroyihsMaxima(ihsMaximaStruct *data)
Destroy vectors and the IHS maxima struct.
Definition: IHS.c:65
ihsVals * createihsVals(void)
Allocate memory for ihsVals struct.
Definition: IHS.c:168
REAL8 ihs2h0(const REAL8 ihsval, const UserInput_t *params)
Convert the IHS statistic to an estimated h0, based on injections.
Definition: IHS.c:1083
INT4 incHarmSumVector(REAL4VectorAligned *output, const REAL4VectorAligned *input, const UINT4 ihsfactor)
Compute the IHS vector – does not compute the maximum value.
Definition: IHS.c:228
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.
Definition: IHS.c:89
void destroyihsVals(ihsVals *ihsvals)
Destroy ihsVals struct.
Definition: IHS.c:180
INT4 incHarmSum(ihsVals *output, const REAL4VectorAligned *input, const UINT4 ihsfactor)
Compute the IHS sum maximum.
Definition: IHS.c:195
ihsfarStruct * createihsfarStruct(const UINT4 rows, const UserInput_t *params)
Allocate memory for ihsfarStruct struct.
Definition: IHS.c:289
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.
Definition: IHS.c:259
REAL4 ihsFOM(const INT4Vector *locs, const INT4 fomnorm)
Calculate the IHS FOM for a number of rows.
Definition: IHS.c:782
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.
Definition: IHS.c:639
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.
Definition: IHS.c:340
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.
Definition: IHS.c:469
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.
Definition: IHS.c:809
ihsMaximaStruct * createihsMaxima(const UINT4 fbins, const UINT4 rows)
Create vectors for IHS maxima struct.
Definition: IHS.c:37
UINT2 B
Definition: SFTnaming.c:47
#define fprintf
FILE * LOG
Definition: TwoSpect.c:55
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.
Definition: candidates.c:1427
candidateVector * resizecandidateVector(candidateVector *vector, const UINT4 length)
Resize a candidateVector.
Definition: candidates.c:61
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.
Definition: candidates.c:122
REAL8 cdf_chisq_Qinv(REAL8 Q, REAL8 nu)
Definition: cdfdist.c:38
#define LAL_E
double REAL8
uint32_t UINT4
int32_t INT4
float REAL4
void * XLALMalloc(size_t n)
void XLALFree(void *p)
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 xlalErrno
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
#define XLAL_CHECK_NULL(assertion,...)
XLAL_EBADLEN
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EINVAL
float data[BLOCKSIZE]
Definition: hwinject.c:360
row
INT4 * data
UINT4 length
REAL4VectorAligned ** data
REAL4VectorAligned * ffdata
REAL8 tfnormalization
INT4Vector * locationsForEachFbin
REAL4VectorAligned * maxima
INT4Vector * locations
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.
Definition: vectormath.c:1292
REAL4VectorAlignedArray * createREAL4VectorAlignedArray(const UINT4 length, const UINT4 vectorLength, const size_t align)
Definition: vectormath.c:99
void destroyREAL4VectorAlignedArray(REAL4VectorAlignedArray *array)
Definition: vectormath.c:110
INT4 VectorSubtractREAL4(REAL4VectorAligned *output, REAL4VectorAligned *input1, REAL4VectorAligned *input2, INT4 vectorMath)
Definition: vectormath.c:243