Commit 2d472ba8 authored by Karl Wette's avatar Karl Wette
Browse files

Merge branch 'RP-fix-negative-2F-values-for-small-number-of-SFTs' into 'master'

add and use new function XLALComputeAntennaPatternSqrtDeterminant()

See merge request !122
parents f26ea14c 5de17d79
Pipeline #14605 failed with stages
in 221 minutes and 56 seconds
......@@ -897,7 +897,8 @@ XLALComputeFstatFromAtoms ( const MultiFstatAtomVector *multiFstatAtoms, ///<
} // loop through detectors
// compute determinant and final Fstat (not twoF!)
REAL4 Dinv = 1.0 / ( mmatrixA * mmatrixB - SQ(mmatrixC) );
REAL4 D = XLALComputeAntennaPatternSqrtDeterminant ( mmatrixA, mmatrixB, mmatrixC, 0);
REAL4 Dinv = 1.0f / D;
twoF = XLALComputeFstatFromFaFb ( Fa, Fb, mmatrixA, mmatrixB, mmatrixC, 0, Dinv );
......
......@@ -34,6 +34,14 @@
#define SQ(x) ( (x) * (x) )
#ifdef __GNUC__
#define likely(x) __builtin_expect(!!(x), 1)
#define unlikely(x) __builtin_expect(!!(x), 0)
#else
#define likely(x) (x)
#define unlikely(x) (x)
#endif
// ---------- Shared global variables ---------- //
// ---------- Shared struct definitions ---------- //
......@@ -79,11 +87,16 @@ XLALComputeFstatFromFaFb ( COMPLEX8 Fa, COMPLEX8 Fb, REAL4 A, REAL4 B, REAL4 C,
REAL4 Fb_re = creal(Fb);
REAL4 Fb_im = cimag(Fb);
REAL4 F = Dinv * ( B * ( SQ(Fa_re) + SQ(Fa_im) )
+ A * ( SQ(Fb_re) + SQ(Fb_im) )
- 2.0 * C * ( Fa_re * Fb_re + Fa_im * Fb_im )
- 2.0 * E * ( - Fa_re * Fb_im + Fa_im * Fb_re ) // nonzero only in RAA case where Ed!=0
);
return 2.0f*F;
REAL4 twoF = 4; // default fallback = E[2F] in noise when Dinv == 0 due to ill-conditionness of M_munu
if ( likely(Dinv > 0) )
{
twoF = 2.0f * Dinv * ( B * ( SQ(Fa_re) + SQ(Fa_im) )
+ A * ( SQ(Fb_re) + SQ(Fb_im) )
- 2.0 * C * ( Fa_re * Fb_re + Fa_im * Fb_im )
- 2.0 * E * ( - Fa_re * Fb_im + Fa_im * Fb_re ) // nonzero only in RAA case where Ed!=0
);
}
return twoF;
} // ComputeFstatFromFaFb()
......@@ -59,9 +59,13 @@
#define SQ(x) (x) * (x)
static REAL4 AntennaPatternMaxCond = 1e4;
static REAL4 AntennaPatternIllCondDeterminant = INFINITY;
/*---------- internal types ----------*/
/*---------- internal prototypes ----------*/
static inline REAL4 estimateAntennaPatternConditionNumber ( REAL4 A, REAL4 B, REAL4 C, REAL4 E );
/*==================== FUNCTION DEFINITIONS ====================*/
......@@ -251,14 +255,8 @@ XLALWeightMultiAMCoeffs ( MultiAMCoeffs *multiAMcoef, const MultiNoiseWeights *
amcoeX->A = AdX;
amcoeX->B = BdX;
amcoeX->C = CdX;
amcoeX->D = AdX * BdX - CdX * CdX - EdX * EdX;
amcoeX->D = XLALComputeAntennaPatternSqrtDeterminant ( AdX, BdX, CdX, EdX );
// in the unlikely event of a degenerate M-matrix with D = sqrt ( det M ) <= 0,
// we set D->inf, in order to set the corresponding F-value to zero rather than >>1
// By setting 'D=inf', we also allow upstream catching/filtering on such singular cases
if ( amcoeX->D <= 0 ) {
amcoeX->D = INFINITY;
}
/* compute multi-IFO antenna-pattern coefficients A,B,C,E by summing over IFOs X */
Ad += AdX;
Bd += BdX;
......@@ -269,15 +267,7 @@ XLALWeightMultiAMCoeffs ( MultiAMCoeffs *multiAMcoef, const MultiNoiseWeights *
multiAMcoef->Mmunu.Ad = Ad;
multiAMcoef->Mmunu.Bd = Bd;
multiAMcoef->Mmunu.Cd = Cd;
multiAMcoef->Mmunu.Dd = Ad * Bd - Cd * Cd - Ed * Ed;
// in the unlikely event of a degenerate M-matrix with D = sqrt(det M_munu) <= 0,
// we set D->inf, in order to set the corresponding F-value to zero rather than >>1
// By setting 'D=inf', we also allow upstream catching/filtering on such singular cases
if ( multiAMcoef->Mmunu.Dd <= 0 ) {
multiAMcoef->Mmunu.Dd = INFINITY;
}
multiAMcoef->Mmunu.Dd = XLALComputeAntennaPatternSqrtDeterminant ( Ad, Bd, Cd, Ed );
if ( multiWeights ) {
multiAMcoef->Mmunu.Sinv_Tsft = multiWeights->Sinv_Tsft;
......@@ -528,3 +518,69 @@ XLALDestroyAMCoeffs ( AMCoeffs *amcoef )
} /* XLALDestroyAMCoeffs() */
/// estimate condition number for given antenna-pattern matrix
static inline REAL4
estimateAntennaPatternConditionNumber ( REAL4 A, REAL4 B, REAL4 C, REAL4 E )
{
REAL4 sumAB = A + B;
REAL4 diffAB = A - B;
REAL4 disc = sqrt ( SQ(diffAB) + 4.0 * ( SQ(C) + SQ(E) ) );
REAL4 denom = sumAB - disc;
REAL4 cond = (denom > 0) ? ((sumAB + disc) / denom) : INFINITY;
return cond;
}
/// Compute (sqrt of) determinant of the antenna-pattern matrix
/// Mmunu = [ A, C, 0, -E;
/// C B E, 0;
/// 0 E A C;
/// -Ed 0 C B; ]
/// which is det(Mmunu) = ( A B - C^2 - E^2 )^2;
///
/// in addition this function checks if the condition number exceeds
/// a module-local tolerance value 'AntennaPatternMaxCond' and outputs a warning
/// and returns NAN if it does.
/// Use XLALSetAntennaPatternMaxCond() to change the acceptable tolerance value.
///
REAL4
XLALComputeAntennaPatternSqrtDeterminant ( REAL4 A, REAL4 B, REAL4 C, REAL4 E )
{
XLAL_CHECK_REAL4 ( (A >= 0) && (B>=0), XLAL_EINVAL );
REAL4 cond = estimateAntennaPatternConditionNumber ( A, B, C, E );
REAL4 det;
if ( cond > AntennaPatternMaxCond )
{
XLALPrintWarning ( "WARNING: Antenna-pattern matrix condition number exceeds maximum allowed value:\n" );
XLALPrintWarning ( "cond{A=%.16g, B=%.16g, C=%.16g, E=%.16g} = %.2e > %.2e ==> setting derminant = %.16g\n",
A, B, C, E, cond, AntennaPatternMaxCond, AntennaPatternIllCondDeterminant );
det = AntennaPatternIllCondDeterminant;
}
else
{
det = A * B - SQ(C) - SQ(E);
}
return det;
}
/// Set a new module-local maximal acceptable condition number of computing
/// antenna-pattern matrix determinant
void
XLALSetAntennaPatternMaxCond ( REAL4 max_cond )
{
AntennaPatternMaxCond = max_cond;
}
/// Set the 'fallback' determinant to use for ill-conditioned antenna-pattern
/// matrix. Use 'INFINITY' (=default) to allow 2F-calculation to continue by
/// 'turning off' the 2F value, or 'NAN' to force an error-condition to terminate
/// the code.
void
XLALSetAntennaPatternIllCondDeterminant ( REAL4 illCondDeterminant )
{
AntennaPatternIllCondDeterminant = illCondDeterminant;
}
......@@ -156,6 +156,9 @@ MultiAMCoeffs *XLALComputeMultiAMCoeffs ( const MultiDetectorStateSeries *multiD
AMCoeffs *XLALCreateAMCoeffs ( UINT4 numSteps );
void XLALDestroyMultiAMCoeffs ( MultiAMCoeffs *multiAMcoef );
void XLALDestroyAMCoeffs ( AMCoeffs *amcoef );
REAL4 XLALComputeAntennaPatternSqrtDeterminant ( REAL4 A, REAL4 B, REAL4 C, REAL4 E );
void XLALSetAntennaPatternMaxCond ( REAL4 max_cond );
void XLALSetAntennaPatternIllCondDeterminant ( REAL4 illCondDeterminant );
/*@}*/
......
......@@ -821,11 +821,8 @@ XLALComputeTransientFstatMap ( const MultiFstatAtomVector *multiFstatAtoms, /**
/* generic F-stat calculation from A,B,C, Fa, Fb */
REAL4 Dd = ( Ad * Bd - Cd * Cd );
REAL4 DdInv = 0;
if ( Dd > 0 ) { /* safety catch as in XLALWeightMultiAMCoeffs(): make it so that in the end F=0 instead of -nan */
DdInv = 1.0f / Dd;
}
REAL4 Dd = XLALComputeAntennaPatternSqrtDeterminant ( Ad, Bd, Cd, 0 );
REAL4 DdInv = 1.0f / Dd;
REAL4 twoF = XLALComputeFstatFromFaFb ( Fa, Fb, Ad, Bd, Cd, 0, DdInv );
REAL4 F = 0.5 * twoF;
/* keep track of loudest F-stat value encountered over the m x n matrix */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment