Fix handling of pulsar glitches in the known pulsar search code. (!1211) · Merge requests · lscsoft / lalsuite

Matthew Pitkin requested to merge matthew-pitkin/lalsuite:fix_glitch into master Feb 05, 2020

Description

The parameter estimation code for the known pulsar search (lalapps_pulsar_parameter_estimation_nested) can deal with parameters describing a pulsar glitch. However, the current implementation makes the unjustified assumption that the glitch(es) was/were not accounted for in the heterodyne stage of the pipeline. In reality many of the glitch parameters will be accounted for in the heterodyne stage. This fix alters the code, so that it initial calculates the phase evolution due to the glitch parameters as used during the heterodyne, and then get the difference in the case that any phase parameter are being included in the parameter estimation (in practice this is only likely to be the glitch phases).

API Changes and Justification

Backwards Compatible Changes

This change introduces no API changes
This change adds new API calls

Backwards Incompatible Changes

This change modifies an existing API
This change removes an existing API

If any of the Backwards Incompatible check boxes are ticked please provide a justification why this change is necessary and why it needs to be done in a backwards incompatible way.

Review Status

The inclusion of glitches in the code was previously looked at in the notebook here. However, in the examples it provided it only really tested lalapps_pulsar_parameter_estimation_nested was internally consistent, i.e., the injected signals and the signal recovery were done with the same bit of code, so glitches were handled identically. This did not necessarily test the reality of the full pipeline.

Here I detail more complete end-to-end tests, using a signal created with independent code and mocked up to include a glitch, and then running this data through the heterodyne and parameter estimation codes. I have done this using the un-fixed code and the code as generated from this patch.

Signal generation

I have generated a signal using lalapps_Makefakedata_v5 using this bash script, reproduced below. It generates two days worth of time series data. The signal in the second day has a frequency and phase offset added to simulate a glitch.

#!/bin/sh

exec=lalapps_Makefakedata_v5

# frame directory
frdir=frames
det=H1
channel=${det}:FAKE_GLITCH
band=128
name=fakeglitch

if [ ! -d $frdir ]; then
   mkdir $frdir
fi

sqrtSn=1e-24
dur=86400
refMJD=55000
refGPS=`python -c "import lalpulsar; print('{0:.9f}'.format(lalpulsar.TTMJDtoGPS(${refMJD})))"`
starttimefloat=`echo $refGPS-$dur | bc`
starttime=${starttimefloat%.*}

# frequency for first half of signal, with no glitch
spinfreq=29.8963954
hnaught=3.0e-24
phinaught=1.0  # initial GW phase
gwfreq=`echo 2*${spinfreq} | bc`
reftime=$(($starttime+$dur))
inj="{Alpha=0; Delta=0; Freq=$gwfreq; refTime=$refGPS; h0=$hnaught; cosi=0.1; psi=0.5; phi0=$phinaught;}"

$exec -F $frdir --outFrChannels=$channel -I $det --sqrtSX=$sqrtSn -G $starttime --duration=$dur --Band=$band --fmin=0 --injectionSources="${inj}" --outLabel=$name

# glitch frequency offset (EM values)
df=0.0001
newfreq=`echo ${gwfreq}+${df}*2 | bc`

# glitch phase offset (GW value, i.e. EM recovered value will be half of this)
dphi=2.0
newphase=`echo ${phinaught}+${dphi} | bc`

# create second half of signal after a "glitch" has occured (10 mins after ref/glitch time) adding glitch phase and frequency offset
newinj="{Alpha=0; Delta=0; Freq=$newfreq; refTime=$refGPS; h0=$hnaught; cosi=0.1; psi=0.5; phi0=$newphase;}"
starttimefloat=`echo $refGPS+600 | bc`
starttime=${starttimefloat%.*}
$exec -F $frdir --outFrChannels=$channel -I $det --sqrtSX=$sqrtSn -G $starttime --duration=$dur --Band=$band --fmin=0 --injectionSources="${newinj}" --outLabel=$name

# write out par file
parfile=J0000+0000.par
echo "PSRJ J0000+0000" > $parfile
echo "RAJ 00:00:00.0" >> $parfile
echo "DECJ 00:00:00.0" >> $parfile
echo "PEPOCH $refMJD" >> $parfile
echo "F0 $spinfreq" >> $parfile
echo "GLEP_1 $refMJD" >> $parfile
echo "GLF0_1 $df" >> $parfile
echo "GLPH_1 0" >> $parfile
echo "EPHEM DE405" >> $parfile
echo "UNITS TDB" >> $parfile

# write out a version without the glitch included
parfileng=J0000+0000_noglitch.par
echo "PSRJ J0000+0000" > $parfileng
echo "RAJ 00:00:00.0" >> $parfileng
echo "DECJ 00:00:00.0" >> $parfileng
echo "PEPOCH $refMJD" >> $parfileng
echo "F0 $spinfreq" >> $parfileng
echo "EPHEM DE405" >> $parfileng
echo "UNITS TDB" >> $parfileng

Heterodyning

I then heterodyne the data using lalapps_heterodyne_pulsar, taking account of the glitch frequency offset, but setting the glitch phase offset to zero (this will be estimated in the parameter estimation stage). The heterodying has been done with this script, reproduced below:

#!/bin/bash

# heterodyne the fake glitch data

# create frame cache file
frfiles=`ls -d $PWD/frames/*`

det=H1
psr=J0000+0000
parfile=J0000+0000.par
channel=${det}:FAKE_GLITCH

cachefile=${det}cache.lcf
segfile=${det}segments.txt

if [ -f $cachefile ]; then
  rm -f $cachefile
fi

if [ -f $segfile ]; then
  rm -f $segfile
fi

# create frame cache file and segment list file
touch $cachefile
touch $segfile
for frfile in $frfiles; do
  noext=${frfile%.*}
  gpstime=$(echo $noext | cut -d "-" -f 3)
  dur=$(echo $noext | cut -d "-" -f 4)
  echo "${det:0:1} $det $gpstime $dur file://localhost${frfile}" >> $cachefile
  echo "$gpstime $(($gpstime+$dur))" >> $segfile
done

# get the sample rate
read -r -a filearray <<< $frfiles
read -r -a samprate <<< `FrChannels ${filearray[0]}`
samplerate=${samprate[1]}

# run coarse heterodyne
coarseoutput=${PWD}/coarse_heterodyne_${psr}_${det}.txt.gz
exec=lalapps_heterodyne_pulsar
$exec -i $det -p $psr -z 0 -f $parfile -k 0.25 -s $samplerate -r 1 -d $cachefile -c $channel -o $coarseoutput -l $segfile -D 1024 -v

# run fine heterodyne
fineoutput=${PWD}/fine_heterodyne_${psr}_${det}.txt.gz
ebase=/home/matthew/repositories/lalsuite/lalpulsar/lib
ephemearth=${ebase}/earth00-40-DE405.dat.gz
ephemsun=${ebase}/sun00-40-DE405.dat.gz
$exec -i $det -p $psr -z 1 -f $parfile -k 0.25 -s 1 -r 1/60 -d $coarseoutput -o $fineoutput -l $segfile -e $ephemearth -S $ephemsun -v

# run fine heterodyne using par file that does not contain the glitch
parfileng=J0000+0000_noglitch.par
fineoutputng=${PWD}/fine_heterodyne_${psr}_${det}_noglitch.txt.gz
$exec -i $det -p $psr -z 1 -f $parfileng -k 0.25 -s 1 -r 1/60 -d $coarseoutput -o $fineoutputng -l $segfile -e $ephemearth -S $ephemsun -v

Parameter estimation

I have then run the lalapps_pulsar_parameter_estimation_nested code to estimate the standard four GW parameters (H0, COSIOTA, PHI0 and PSI) and the glitch phase offset (GLPH_1), which in this case has a value (for the EM-based parameter) of 1 radian. The current version of the code has been run with identical set up to the version based off this MR. These runs have been performed with this this script, reproduced below:

#!/bin/sh

# perform parameter estimation in various cases

det=H1
psr=J0000+0000
parfile=J0000+0000.par

# fine heterodyned data file
datafile=${PWD}/fine_heterodyne_${psr}_${det}.txt.gz

ebase=/home/matthew/repositories/lalsuite/lalpulsar/lib
ephemearth=${ebase}/earth00-40-DE405.dat.gz
ephemsun=${ebase}/sun00-40-DE405.dat.gz
ephemtime=${ebase}/tdb_2000-2040.dat.gz

# create prior file
priorfile=priors.txt
echo "H0 uniform 0 5e-23" > $priorfile
echo "PHI0 uniform 0 3.141592653589793" >> $priorfile
echo "COSIOTA uniform -1.0 1.0" >> $priorfile
echo "PSI uniform 0.0 1.5707963267948966" >> $priorfile
echo "GLPH_1 uniform 0 1" >> $priorfile

# set parameters
Nlive=1000
Ntraining=2000
outfile=nest_old.hdf5
postfile=post_old.hdf5

if [ -f $outfile ]; then
  rm -f $outfile
fi

# run parameter estimation with the old (un-fixed) code
execold=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37/bin/lalapps_pulsar_parameter_estimation_nested

$execold --detectors $det --par-file $parfile --input-files $datafile --outfile $outfile --prior-file $priorfile \
         --ephem-earth $ephemearth --ephem-sun $ephemsun --ephem-timecorr $ephemtime --Nlive $Nlive \
         --Nmcmcinitial 0 --roq --ntraining $Ntraining --verbose

ntwopexec=/cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/envs/igwn-py37/bin/lalinference_nest2pos

if [ -f $postfile ]; then
  rm -f $postfile
fi

$ntwopexec -p $postfile $outfile

# run parameter estimation with the fixed code
outfile=nest_new.hdf5
postfile=post_new.hdf5
execnew=/home/matthew/miniconda2/envs/lalsuite-fix_glitch/bin/lalapps_pulsar_parameter_estimation_nested

if [ -f $outfile ]; then
  rm -f $outfile
fi

$execnew --detectors $det --par-file $parfile --input-files $datafile --outfile $outfile --prior-file $priorfile \
         --ephem-earth $ephemearth --ephem-sun $ephemsun --ephem-timecorr $ephemtime --Nlive $Nlive \
         --Nmcmcinitial 0 --roq --ntraining $Ntraining --verbose

if [ -f $postfile ]; then
  rm -f $postfile
fi

$ntwopexec -p $postfile $outfile

Results

The results of performing parameter estimation using the old and new code can be found in the notebook here, which can be more easily viewed at this site.

I will reproduce the two main figures below. In the first, the parameter estimation on the four GW parameters and glitch phase offset is performed using the current version of the code. In it the signal amplitude recovery is wrong, with it being about half of the true value, and the glitch phase is not recovered. In the second, the parameter estimation is performed using the fix in this MR. In this case all parameters, including the GW amplitude and glitch phase offset are recovered correctly.

Notes

This code will be required for the Crab O3a search results, which we would like to have finalised soon for the paper.

cc @karl-wette @david-keitel @teviet.creighton @damir.buskulic

Edited Jul 02, 2022 by Ghost User

Fix handling of pulsar glitches in the known pulsar search code.