Commit 22780605 authored by Karl Wette's avatar Karl Wette
Browse files

Merge branch 'modularize-ComputePSD' into 'master'

lalapps_ComputePSD overhaul: modularize, fix some minor bugs, add one new feature

Closes CW/software/lalsuite#79 and CW/software/lalsuite#78

See merge request lscsoft/lalsuite!1294
parents 2e8b60e2 96a1d2ad
This diff is collapsed.
......@@ -2,21 +2,36 @@
psd_code="lalapps_ComputePSD"
mfd_code="lalapps_Makefakedata_v4"
tolerance=1e-5
# ---------- fixed parameter of our test SFTs
IFO=H1
retstatus=0
# ---------- fixed parameters of our test SFTs
sqrtSh=3.25e-22
h0=5e-23
cosi=0
startTime=828002611
fmin=50
linefreq="50.05"
TSFT=1800
# ---------- fixed parameters for PSD runs
blocksRngMed=101
outPSD=psd1.dat
## ----- Correct cropping of normalized SFT
# ---------- define comparison functions
awk_absdiff='{printf "%.6f", sqrt(($1-$2)*($1-$2)) }'
awk_reldev='{printf "%.2e", sqrt(($1-$2)*($1-$2))/(0.5*($1+$2)) }'
awk_isgtr='{if($1>$2) {print "1"}}'
awk_nonpos='{if(0.0>=$1) {print "1"}}'
echo "ComputePSD: --outputNormSFT frequency range..."
## ----- test correct cropping of normalized SFT power on a single SFT
echo "----------------------------------------------------------------------"
echo "STEP 1: testing --outputNormSFT frequency range..."
echo "----------------------------------------------------------------------"
echo
outSFT="./testpsd_sft_H1"
linefreq="50.05"
## ----- run MFDv4
cmdline="${mfd_code} --IFO=$IFO --outSingleSFT=1 --outSFTbname=$outSFT --startTime=828002611 --duration=1800 --fmin=50 --Band=0.1 --noiseSqrtSh=3.25e-22 --lineFeature=1 --h0=5e-23 --cosi=0 --Freq=$linefreq --randSeed=1 "
## ----- create SFTs
cmdline="${mfd_code} --IFO=H1 --outSingleSFT=1 --outSFTbname=$outSFT --startTime=$startTime --duration=$TSFT --fmin=$fmin --Band=0.1 --noiseSqrtSh=$sqrtSh --lineFeature=1 --h0=$h0 --cosi=$cosi --Freq=$linefreq --randSeed=1"
echo $cmdline;
if ! eval $cmdline; then
......@@ -45,6 +60,7 @@ if ! eval $cmdline; then
exit 1
fi
# get the result file row at maximum, which should recover the line feature
topline_psd_full=$(sort -nr -k2,2 $outPSD_full | head -1)
toppsd_full=$(echo $topline_psd_full | awk '{print $2}')
toppsdfreq_full=$(echo $topline_psd_full | awk '{print $1}')
......@@ -63,14 +79,8 @@ echo "Loudest bins:"
echo "==> full SFT: PSD=$toppsd_full at $toppsdfreq_full Hz, normPower=$toppower_full at $toppowerfreq_full Hz"
echo "==> freqband: PSD=$toppsd_band at $toppsdfreq_band Hz, normPower=$toppower_band at $toppowerfreq_band Hz"
retstatus=0
tolerance_freq=0.0
tolerance_psd=1e-46
tolerance_power=0.001
tolerance_rel=0.001
awk_absdiff='{printf "%.6f", sqrt(($1-$2)*($1-$2)) }'
awk_reldev='{printf "%.2e", sqrt(($1-$2)*($1-$2))/(0.5*($1+$2)) }'
awk_isgtr='{if($1>$2) {print "1"}}'
freqdiff_full=$(echo $toppowerfreq_full $linefreq | awk "$awk_absdiff")
freqdiff_band=$(echo $toppowerfreq_band $linefreq | awk "$awk_absdiff")
freqdiff_runs=$(echo $toppowerfreq_full $toppowerfreq_band | awk "$awk_absdiff")
......@@ -90,7 +100,204 @@ if [ "$fail1" -o "$fail2" -o "$fail3" -o "$fail4" -o "$fail5" ]; then
else
echo " ==> OK"
echo
echo "========== OK. All ComputePSD tests PASSED. =========="
echo "========== OK. Line feature successfully recovered. =========="
echo
fi
echo "----------------------------------------------------------------------"
echo "STEP 2: testing math ops over SFTs and IFOs..."
echo "----------------------------------------------------------------------"
echo
tolerance_rel=0.000001
## ----- create SFTs
outSFTbname="./testpsd_multisfts"
IFOs=(H1 L1)
for X in 0 1; do
outSFT="${outSFTbname}_${IFOs[$X]}"
cmdline="${mfd_code} --IFO=${IFOs[$X]} --outSingleSFT=1 --outSFTbname=$outSFT --startTime=$startTime --duration=$(echo "2*$TSFT" | bc) --fmin=$fmin --Band=0.1 --noiseSqrtSh=$sqrtSh --lineFeature=1 --h0=$h0 --cosi=$cosi --Freq=$linefreq --randSeed=$X"
echo $cmdline;
if ! eval $cmdline; then
echo "Error.. something failed when running '$mfd_code' ..."
exit 1
fi
done
echo
# function to de-spaghettify the loop over multiple mathops, IFOs and timestamps below
get_psd () { # expected argument order: mthop IFO startTime endTime extrArgs
psdfile="psd_$2_$3-$4.dat"
# here we only use the --IFO constraint for single IFOs;
# currently the code would silently interpret H1L1 as H1 only!
if [ ${#2} == 2 ]; then
IFObit="--IFO=$2"
else
IFObit=""
fi
cmdline="${psd_code} ${IFObit} --inputData=${outSFTbname}* --outputPSD=$psdfile --blocksRngMed=$blocksRngMed --outputNormSFT=1 --startTime=$3 --endTime=$4 --PSDmthopSFTs=$1 --PSDmthopIFOs=$1 $5"
echo $cmdline;
if ! eval $cmdline; then
echo "Error.. something failed when running '$psd_code' ..."
exit 1
fi
firstline_psd=$(grep -o '^[^%]*' $psdfile | head -1)
freq=$(echo $firstline_psd | awk '{print $1}')
psd=$(echo $firstline_psd | awk '{print $2}')
power=$(echo $firstline_psd | awk '{print $3}')
firstline_psd=($freq $psd $power)
}
## ----- test the various mthops - for simplicity, use the same over SFTs and IFOs
# use the same ordering as the MathOpType enum:
mthops_awk=('{printf "%.6e", $1+$2}' # ARITHMETIC_SUM
'{printf "%.6e", 0.5*($1+$2)}' # ARITHMETIC_MEAN
'{printf "%.6e", 0.5*($1+$2)}' # ARITHMETIC_MEDIAN - note for 2 elements this falls back to mean!
'{printf "%.6e", 1./(1./$1+1./$2)}' # HARMONIC_SUM
'{printf "%.6e", 2./(1./$1+1./$2)}' # HARMONIC_MEAN
'{printf "%.6e", 1./sqrt(1./($1*$1)+1./($2*$2))}' # POWERMINUS2_SUM
'{printf "%.6e", 1./sqrt(0.5*(1./($1*$1)+1./($2*$2)))}' # POWERMINUS2_MEAN
'{printf "%.6e", ($1<$2) ? $1 : $2}' # MINIMUM
'{printf "%.6e", ($1>$2) ? $1 : $2}' # MAXIMUM
)
for mthop in 0 1 2 3 4 5 6 7 8; do
startTime2=$(echo "$startTime + $TSFT" | bc)
endTime1=$(echo "$startTime + 1" | bc)
endTime2=$(echo "$startTime2 + 1" | bc)
echo "Running ComputePSD with PSDmthopSFTs=PSDmthopIFOs=$mthop..."
for IFO in "H1" "L1" "H1L1"; do
for minT in $startTime $startTime2; do
for maxT in $endTime1 $endTime2; do
if (( $maxT > $minT )); then
get_psd $mthop $IFO $minT $maxT
# storing the results in dynamically named variables
# so they can be accessed again outside the loop
psdvar="psd_${IFO}_${minT}_${maxT}"
powvar="power_${IFO}_${minT}_${maxT}"
declare $psdvar=${firstline_psd[1]}
declare $powvar=${firstline_psd[2]}
fi
done
done
done
echo "--------- Compare results (PSDmthopSFTs=PSDmthopIFOs=$mthop) ----------------------------------------------------------------------------------------"
echo " [T0,T1] [T1,T2] [T0,T2] awk (reldev) [Tolerance=$tolerance_rel]"
failed=0
# first check the over-SFTs mathop for each IFO (and multi-IFO)
for IFO in "H1" "L1" "H1L1"; do
psd1=$(eval "echo \$"$(echo "psd_${IFO}_${startTime}_${endTime1}")"")
psd2=$(eval "echo \$"$(echo "psd_${IFO}_${startTime2}_${endTime2}")"")
psd12=$(eval "echo \$"$(echo "psd_${IFO}_${startTime}_${endTime2}")"")
psdawk=$(echo $psd1 $psd2 | awk "${mthops_awk[$mthop]}")
psd_reldev=$(echo $psd12 $psdawk | awk "$awk_reldev")
nonpos_psd1=$(echo $psd1 | awk "$awk_nonpos")
nonpos_psd2=$(echo $psd2 | awk "$awk_nonpos")
nonpos_psd12=$(echo $psd12 | awk "$awk_nonpos")
failed_tolerance=$(echo $psd_reldev $tolerance_rel | awk "$awk_isgtr")
if [ "$nonpos_psd1" -o "$nonpos_psd2" -o "$nonpos_psd12" -o "$failed_tolerance" ]; then
failbit="FAILED!"
failed=1
else
failbit="OK!"
fi
echo "==> $IFO: $psd1 $psd2 $psd12 $psdawk ($psd_reldev) $failbit"
done
# then check the over-IFOs mathop (only for the full SFT set)
psdH1=$(eval "echo \$"$(echo "psd_H1_${startTime}_${endTime2}")"")
psdL1=$(eval "echo \$"$(echo "psd_L1_${startTime}_${endTime2}")"")
psdH1L1=$(eval "echo \$"$(echo "psd_H1L1_${startTime}_${endTime2}")"")
psdawk=$(echo $psdH1 $psdL1 | awk "${mthops_awk[$mthop]}")
psd_reldev=$(echo $psdH1L1 $psdawk | awk "$awk_reldev")
nonpos_psdH1L1=$(echo $psdH1L1 | awk "$awk_nonpos")
failed_tolerance=$(echo $psd_reldev $tolerance_rel | awk "$awk_isgtr")
if [ "$nonpos_psdH1L1" -o "$failed_tolerance" ]; then
failbit="FAILED!"
failed=1
else
failbit="OK!"
fi
echo "==> H1L1(awk): $psdH1L1 $psdawk ($psd_reldev) $failbit"
if (( $failed > 0 )); then
echo " ==> something FAILED"
retstatus=1
else
echo "========== OK. All results consistent. =========="
echo
fi
done # loop over mthops
echo "----------------------------------------------------------------------"
echo "STEP 3: testing special case for combined all-IFO-all-SFTs harmmean/powerminus2mean..."
echo "----------------------------------------------------------------------"
echo
## ----- test the special normalizeByTotalNumSFTs mode
mthops_loop=('harmsum' # together with --normalizeByTotalNumSFTs should emulate harmmean
'powerminus2sum' # together with --normalizeByTotalNumSFTs should emulate powerminus2mean
)
mthops_awk=('{printf "%.6e", 4./(1./$1+1./$2+1./$3+1./$4)}' # harmonic mean over 4 inputs
'{printf "%.6e", sqrt(4.)/sqrt((1./($1*$1)+1./($2*$2)+1./($3*$3)+1./($4*$4)))}' # power-minus-2 mean over 4 inputs
)
for mthopidx in 0 1; do
startTime2=$(echo "$startTime + $TSFT" | bc)
endTime1=$(echo "$startTime + 1" | bc)
endTime2=$(echo "$startTime2 + 1" | bc)
echo "Running ComputePSD with PSDmthopSFTs=PSDmthopIFOs=${mthops_loop[$mthopidx]} and --normalizeByTotalNumSFTs..."
for IFO in "H1" "L1" "H1L1"; do
for minT in $startTime $startTime2; do
for maxT in $endTime1 $endTime2; do
if (( $maxT > $minT )); then
get_psd ${mthops_loop[$mthopidx]} $IFO $minT $maxT "--normalizeByTotalNumSFTs"
# storing the results in dynamically named variables
# so they can be accessed again outside the loop
psdvar="psd_${IFO}_${minT}_${maxT}"
powvar="power_${IFO}_${minT}_${maxT}"
declare $psdvar=${firstline_psd[1]}
declare $powvar=${firstline_psd[2]}
fi
done
done
done
echo "--------- Compare results (PSDmthopSFTs=PSDmthopIFOs=${mthops_loop[$mthopidx]} --normalizeByTotalNumSFTs) ----------------------------------------------------------------------------------------"
echo " [T0,T1] [T1,T2] [T0,T2] awk (reldev) [Tolerance=$tolerance_rel]"
for IFO in "H1" "L1" "H1L1"; do
psd1=$(eval "echo \$"$(echo "psd_${IFO}_${startTime}_${endTime1}")"")
psd2=$(eval "echo \$"$(echo "psd_${IFO}_${startTime2}_${endTime2}")"")
echo "==> $IFO: $psd1 $psd2"
done
failed=0
# check the combination over the full all-IFOs-all-SFTs set
psdH11=$(eval "echo \$"$(echo "psd_H1_${startTime}_${endTime1}")"")
psdH12=$(eval "echo \$"$(echo "psd_H1_${startTime2}_${endTime2}")"")
psdL11=$(eval "echo \$"$(echo "psd_L1_${startTime}_${endTime1}")"")
psdL12=$(eval "echo \$"$(echo "psd_L1_${startTime2}_${endTime2}")"")
psdH1L1=$(eval "echo \$"$(echo "psd_H1L1_${startTime}_${endTime2}")"")
psdawk=$(echo $psdH11 $psdH12 $psdL11 $psdL12 | awk "${mthops_awk[$mthopidx]}")
psd_reldev=$(echo $psdH1L1 $psdawk | awk "$awk_reldev")
nonpos_psdH1L1=$(echo $psdH1L1 | awk "$awk_nonpos")
failed_tolerance=$(echo $psd_reldev $tolerance_rel | awk "$awk_isgtr")
if [ "$nonpos_psdH1L1" -o "$failed_tolerance" ]; then
failbit="FAILED!"
failed=1
else
failbit="OK!"
fi
echo "==> H1L1(awk): $psdH1L1 $psdawk ($psd_reldev) $failbit"
if (( $failed > 0 )); then
echo " ==> something FAILED"
retstatus=1
else
echo "========== OK. All results consistent. =========="
echo
fi
done # loop over mthops
exit $retstatus
......@@ -182,7 +182,32 @@ int XLALCWGPSinRange( const LIGOTimeGPS gps, const LIGOTimeGPS* minGPS, const LI
return 1;
}
return 0;
}
} /* XLALCWGPSinRange() */
/**
* Round a REAL8 frequency down to the nearest integer SFT bin number.
*
* This function provides an official rounding convention,
* including a "fudge" factor.
*/
UINT4 XLALRoundFrequencyDownToSFTBin( const REAL8 freq, const REAL8 df )
{
return (UINT4) floor (freq / df * fudge_up);
} /* XLALRoundFrequencyDownToSFTBin() */
/**
* Round a REAL8 frequency up to the nearest integer SFT bin number.
*
* This function provides an official rounding convention,
* including a "fudge" factor.
*/
UINT4 XLALRoundFrequencyUpToSFTBin( const REAL8 freq, const REAL8 df )
{
return (UINT4) ceil (freq / df * fudge_down);
} /* XLALRoundFrequencyUpToSFTBin() */
/**
* Find the list of SFTs matching the \a file_pattern and satisfying the given \a constraints,
......@@ -740,11 +765,11 @@ XLALLoadSFTs (const SFTCatalog *catalog, /**< The 'catalogue' of SFTs to load
if (fMin < 0)
firstbin = minbin;
else
firstbin = (UINT4) floor (fMin / deltaF * fudge_up); // round *down*, but allow for 10*eps 'fudge'
firstbin = XLALRoundFrequencyDownToSFTBin ( fMin, deltaF );
if (fMax < 0)
lastbin = maxbin;
else {
lastbin = (UINT4) ceil (fMax / deltaF * fudge_down); // round *up*, but allow for 10*eps fudge
lastbin = XLALRoundFrequencyUpToSFTBin ( fMax, deltaF );
if((lastbin == 0) && (fMax != 0)) {
XLALPrintError("ERROR: last bin to read is 0 (fMax: %f, deltaF: %f)\n", fMax, deltaF);
XLALLOADSFTSERROR(XLAL_EINVAL);
......
......@@ -286,6 +286,9 @@ typedef struct tagMultiSFTCatalogView
int XLALCWGPSinRange( const LIGOTimeGPS gps, const LIGOTimeGPS* minGPS, const LIGOTimeGPS* maxGPS );
UINT4 XLALRoundFrequencyDownToSFTBin( const REAL8 freq, const REAL8 df );
UINT4 XLALRoundFrequencyUpToSFTBin( const REAL8 freq, const REAL8 df );
LALStringVector *XLALFindFiles (const CHAR *globstring);
SFTCatalog *XLALSFTdataFind ( const CHAR *file_pattern, const SFTConstraints *constraints );
......
This diff is collapsed.
......@@ -29,14 +29,17 @@ extern "C" {
/**
* \defgroup SFTutils_h Header SFTutils.h
* \ingroup lalpulsar_sft
* \author Reinhard Prix, Badri Krishnan
* \date 2005
* \brief Utility functions for handling of SFTtype and SFTVectors
* \author Reinhard Prix, Badri Krishnan, Badri Krishnan, Iraj Gholami,
* Reinhard Prix, Alicia Sintes, Karl Wette, David Keitel
* \date 2005-2020
* \brief Utility functions for handling of SFTs and associated structures
*
* The helper functions XLALCreateSFT(), XLALDestroySFT(), XLALCreateSFTVector()
* and XLALDestroySFTVector() respectively allocate and free SFT-structs and SFT-vectors.
* Similarly, XLALCreateTimestampVector() and XLALDestroyTimestampVector() allocate and free
* a bunch of GPS-timestamps.
* This module contains various helper functions to create, handle, combine,
* and destroy SFTs (Short Fourier Transforms) and related data structures,
* including SFTtype, SFTVector, SFTCatalog
* (and their multi-detector generalizations)
* as well as tools for dealing with timestamps, segments
* and ASD/PSD (Amplitude/Power Spectral Density) estimates.
*
*/
/** @{ */
......@@ -53,6 +56,7 @@ extern "C" {
#include <lal/Segments.h>
#include <lal/SFTfileIO.h>
#include <lal/StringVector.h>
#include <lal/UserInputParse.h>
/*---------- DEFINES ----------*/
......@@ -91,6 +95,22 @@ typedef struct tagMultiNoiseWeights {
BOOLEAN isNotNormalized; /**< if true: weights are saved unnormalized (divide by Sinv_Tsft to get normalized version). */
} MultiNoiseWeights;
/** common types of mathematical operations over an array */
typedef enum tagMathOpType {
MATH_OP_ARITHMETIC_SUM = 0, /**< \f$\sum_k x_k\f$ */
MATH_OP_ARITHMETIC_MEAN, /**< \f$\sum_k x_k / N\f$ */
MATH_OP_ARITHMETIC_MEDIAN, /**< \f$x_1 \leq \dots \leq x_{N/2} \leq \dots \leq x_n\f$ */
MATH_OP_HARMONIC_SUM, /**< \f$1 / \sum_k (1/x_k)\f$ */
MATH_OP_HARMONIC_MEAN, /**< \f$N / \sum_k (1/x_k)\f$ */
MATH_OP_POWERMINUS2_SUM, /**< \f$1 / \sqrt{ \sum_k (1/x_k^2) }\f$ */
MATH_OP_POWERMINUS2_MEAN, /**< \f$1 / \sqrt{ \sum_k (1/x_k^2) / N }\f$ */
MATH_OP_MINIMUM, /**< \f$\min_k(x_k)\f$ */
MATH_OP_MAXIMUM, /**< \f$\max_k(x_k)\f$ */
MATH_OP_LAST
} MathOpType;
extern const UserChoices MathOpTypeChoices;
/*---------- Global variables ----------*/
/*---------- exported prototypes [API] ----------*/
......@@ -171,6 +191,40 @@ int XLALFindTimesliceBounds ( UINT4 *iStart, UINT4 *iEnd, const LIGOTimeGPSVecto
SFTVector *XLALExtractSFTVectorWithTimestamps ( const SFTVector *sfts, const LIGOTimeGPSVector *timestamps );
MultiSFTVector *XLALExtractMultiSFTVectorWithMultiTimestamps ( const MultiSFTVector *multiSFTs, const MultiLIGOTimeGPSVector *multiTimestamps );
// compute and work with PSDs
int XLALComputePSDandNormSFTPower ( REAL8Vector **finalPSD,
MultiPSDVector **multiPSDVector,
REAL8Vector **normSFT,
MultiSFTVector *inputSFTs,
const BOOLEAN returnMultiPSDVector,
const BOOLEAN returnNormSFT,
const UINT4 blocksRngMed,
const MathOpType PSDmthopSFTs,
const MathOpType PSDmthopIFOs,
const MathOpType nSFTmthopSFTs,
const MathOpType nSFTmthopIFOs,
const BOOLEAN normalizeByTotalNumSFTs,
const REAL8 FreqMin,
const REAL8 FreqBand
);
int XLALDumpMultiPSDVector ( const CHAR *outbname, const MultiPSDVector *multiPSDVect );
int XLALCropMultiPSDandSFTVectors ( MultiPSDVector *multiPSDVect, MultiSFTVector *multiSFTVect, UINT4 firstBin, UINT4 lastBin );
REAL8FrequencySeries *XLALComputeSegmentDataQ ( const MultiPSDVector *multiPSDVect, LALSeg segment );
REAL8 XLALMathOpOverArray(const REAL8* data, const size_t length, const MathOpType optype);
REAL8 XLALGetMathOpNormalizationFactorFromTotalNumberOfSFTs ( const UINT4 totalNumSFTs, const MathOpType optypeSFTs );
int XLALWritePSDtoFilePointer ( FILE *fpOut,
REAL8Vector *PSDVect,
REAL8Vector *normSFTVect,
BOOLEAN outputNormSFT,
BOOLEAN outFreqBinEnd,
INT4 PSDmthopBins,
INT4 nSFTmthopBins,
INT4 binSize,
INT4 binStep,
REAL8 Freq0,
REAL8 dFreq
);
/** @} */
#ifdef __cplusplus
......
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