PyFstat / PyFstat

a python package for gravitational wave analysis with the F-statistic
https://pyfstat.readthedocs.io
MIT License
47 stars 28 forks source link

lalpulsar_Makefakedata_v5 error for high F1 at fixed Band #501

Closed viktorcikojevic closed 1 year ago

viktorcikojevic commented 1 year ago

I was trying to run the python code to generate the signals. During the execution of a python script, I get the error subprocess.CalledProcessError: Command 'lalpulsar_Makefakedata_v5 ... returned non-zero exit status 255.

When I try to run the process directly in the terminal,

lalpulsar_Makefakedata_v5 --outSingleSFT=TRUE --outSFTdir="~/viktor/T7/gravitational-waves/generated-signals/TMP-v5/signals_903553730_10373__8.671003685899436_H1_True" --outLabel="PyFstat" --IFOs="H1" --sqrtSX="9.518966822812071e-26" --SFTWindowType="tukey" --SFTWindowBeta=0.01 --timestampsFiles=/media/viktor/T7/gravitational-waves/generated-signals/TMP-v5/signals_903553730_10373__8.671003685899436_H1_True/PyFstat_timestamps_H1.csv --fmin=123.9002777777778 --Band=0.1994444444444445 --Tsft=1800 --injectionSources="/media/viktor/T7/gravitational-waves/generated-signals/TMP-v5/signals_903553730_10373__8.671003685899436_H1_True/PyFstat.cff" --ephemEarth="earth00-40-DE405.dat.gz" --ephemSun="sun00-40-DE405.dat.gz"

I got the following error message

XLAL Error - XLALCWMakeFakeData (CWMakeFakeData.c:327): Error: injection signal 0:'/media/viktor/T7/gravitational-waves/generated-signals/TMP-v5/s' needs frequency band [123.986862,213.846828]Hz, injecting into [123.900000,124.100000]Hz

The PyFstat.cff file that my python script writes reads

[TS0]
Alpha = 4.543341133226692641e+00
Delta = -2.219996285055303997e-01
h0 = 1.097792962340861523e-26
cosi = -6.328996314100563403e-01
psi = -5.365079206060429096e-01
transientWindowType = none
refTime = 1238172451.000000000
Freq = 1.240000000000000000e+02
f1dot = 8.663720666732984996e-06
f2dot = 0.000000000000000000e+00
phi0 = 5.639301306088492893e+00

When I run the same command line using the fmin = 123.986862, as suggested by the error message, everything works well.

I don't understand why, everything seems correctly set up. The python script that generates this data has the following parameters.

...
writer_kwargs = {
                  "detectors": detector, # "H1,L1",        
                  "sqrtSX": h0 * depth, # noise level
                  "Tsft": 1800,            
                  "SFTWindowType": "tukey",
                  "SFTWindowBeta": 0.01,
                  "timestamps": timestamps
                  }

  signal_parameters_generator = pyfstat.AllSkyInjectionParametersGenerator(
      priors={
          "F0": np.random.randint(50,500),
          "Band": 359/1800 , #
          "F1": lambda: 10**stats.uniform(-9, 4).rvs(),
          "F2": 0,
          "h0": h0 if draw_signal else 0, # GW strain
          **pyfstat.injection_parameters.isotropic_amplitude_priors,
      },
      seed=seed
  )
...

Can you tell me how to set up F0 correctly in the python script so that the error does not occur anymore? Thanks

dbkeitel commented 1 year ago

Hi @viktorcikojevic there's two problems here:

  1. You set "Band": 359/1800 when creating the signal_parameters_generator, but that is not a "signal parameter" in the physical sense. Band determines how broad the SFTs are, so it should more logically go into the writer_kwargs. There is no failure either from instantiating the AllSkyInjectionParametersGenerator (because it is set up to understand and return arbitrary parameters) nor the Writer (presumably because you do writer = pyfstat.Writer(**writer_kwargs, **signal_parameters) and so it doesn't matter which dict the Band comes from).
  2. Your hand-chosen Band is then actually physically too narrow for the signal you want to inject - f1dot = 8.663720666732984996e-06 is HUGE by continuous waves standards, and so over whatever duration you want to simulate (can't see your timestamps) the code internally notices that the signal would evolve over nearly 90Hz, which does not fit into the 359/1800 you're providing. If you simply pass no Band option, PyFstat will auto-estimate a Band for you, which most of the time should work correctly. Though we haven't checked too much how robust it is for really high spindowns.

So I would suggest to fix 1 just for cleaner code, and regarding 2 either you really want very strong spindown then change to auto-banding or choose a very wide Band by hand (but be aware that files will be big); or simply limit your F1 random range.

viktorcikojevic commented 1 year ago

Hi @dbkeitel, thanks for the prompt answer, now I see it more clearly.

The timestamps range between 1,238,173,136 and 1,248,557,243.

So if I understood correctly, the maximum value of f1dot should be about 0.2 / 10^6, since 0.2 is the fixed Band width, and 10^6 is approximately the difference between the end and start timestamp. Is that correct? 🙂

dbkeitel commented 1 year ago

Looks more like 10^7 to me, but otherwise that's the right logic, yes!

And if you want to include F2, that affects the frequency evolution (and hence whether the signal fits into a band) proportional to time difference squared.

viktorcikojevic commented 1 year ago

Great, thanks!