nanograv / enterprise

ENTERPRISE (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE) is a pulsar timing analysis code, aimed at noise analysis, gravitational-wave searches, and timing model analysis.
https://enterprise.readthedocs.io
MIT License
64 stars 65 forks source link

ValueError when creating PTA object with some PINT Pulsar objects #296

Closed AaronDJohnson closed 2 years ago

AaronDJohnson commented 2 years ago

For B1937+21 and J1713+0747 I am getting the following error when running

pta_pint = models.model_2a([psr], noisedict=noise_params, gamma_common=4.33).

The error starts with 07e1150673b44f8727add92890238b1564ea040c. Removing ecorr from the model fixes the error.

ValueError: ERROR: TOAs not sorted properly!
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/var/folders/2k/k752fw_d4gd91qly20lhwgg40000gn/T/ipykernel_18288/1542723239.py in <module>
----> 1 pta_pint = models.model_2a([psrs[1]], noisedict=noise_params, gamma_common=4.33, bayesephem=True)
      2 pta_t2 = models.model_2a([t2_psrs[1]], noisedict=noise_params, gamma_common=4.33, bayesephem=True)

~/mambaforge/envs/test/lib/python3.9/site-packages/enterprise_extensions/models.py in model_2a(psrs, psd, noisedict, components, n_rnfreqs, n_gwbfreqs, gamma_common, delta_common, upper_limit, bayesephem, be_type, white_vary, is_wideband, use_dmdata, select, pshift, pseed, psr_models, tm_marg, dense_like)
    569             s2 = s + white_noise_block(vary=white_vary, inc_ecorr=True,
    570                                        select=select)
--> 571             models.append(s2(p))
    572         else:
    573             s3 = s + white_noise_block(vary=white_vary, inc_ecorr=False,

~/Documents/GitHub/enterprise/enterprise/signals/signal_base.py in __init__(self, psr)
    730             self.psrname = psr.name
    731             # instantiate all the signals with a pulsar
--> 732             self._signals = [metasignal(psr) for metasignal in self._metasignals]
    733 
    734             self._residuals = psr.residuals

~/Documents/GitHub/enterprise/enterprise/signals/signal_base.py in <listcomp>(.0)
    730             self.psrname = psr.name
    731             # instantiate all the signals with a pulsar
--> 732             self._signals = [metasignal(psr) for metasignal in self._metasignals]
    733 
    734             self._residuals = psr.residuals

~/Documents/GitHub/enterprise/enterprise/signals/white_signals.py in __init__(self, psr)
    173                 nn = Umats[ct].shape[1]
    174                 U[mask, netot : nn + netot] = Umats[ct]
--> 175                 self._slices.update({key: utils.quant2ind(U[:, netot : nn + netot])})
    176                 netot += nn
    177 

~/Documents/GitHub/enterprise/enterprise/signals/utils.py in quant2ind(U)
    694         epinds = np.flatnonzero(col)
    695         if epinds[-1] - epinds[0] + 1 != len(epinds):
--> 696             raise ValueError("ERROR: TOAs not sorted properly!")
    697         inds.append(slice(epinds[0], epinds[-1] + 1))
    698     return inds

ValueError: ERROR: TOAs not sorted properly!
Hazboun6 commented 2 years ago

@AaronDJohnson the command you used here is an e_e command. Can you say a little bit about why you think this is an enterprise issue? Also, we should assign Michele since it is his commit.

AaronDJohnson commented 2 years ago

@AaronDJohnson the command you used here is an e_e command. Can you say a little bit about why you think this is an enterprise issue? Also, we should assign Michele since it is his commit.

The error also happens when models.append(s(p)) is called without e_e. I have another file which uses the same model but in "verbose" format using only enterprise (a modified 12.5 year stochastic tutorial file). By dropping out signals from this file, I was able to determine that ecorr was causing the error.

Hazboun6 commented 2 years ago

@vallis do you have any idea how that particular commit could be involved with this error?

Hazboun6 commented 2 years ago

I think the big clue here is that the errors are happening on B1937 and J1713, which are the only two pulsars with TOAs from both telescopes. Michele was working on being able to split the TOAs to look for ORFs separately at different telescopes.

vallis commented 2 years ago

Hmm... my changes to split common GPs by telescope are still in a PR and not merged yet.

The function quant2ind that fails says that it only works for sorted TOAs. Both Pulsar versions do the sorting. But perhaps having two telescopes the vector still looks unsorted?

vallis commented 2 years ago

Does this problem occur using Nihan's pickles or reading from par/tim?

AaronDJohnson commented 2 years ago

@vallis as far as I can tell it must only occur with the par/tim files from the 12.5 year data set. Nihan's pickles must have been free of this problem, because the OSU group used the newest version of enterprise (as far as I know) to analyze the pickles the other day.

vallis commented 2 years ago

Can you reproduce the problem by loading either B1937 or J1713 from par/tim, and instantiating the ECORR model?

AaronDJohnson commented 2 years ago

error_example.zip This contains the 12.5 year par and tim file for B1937 along with the following code to reproduce the error:


import glob
from enterprise.pulsar import Pulsar
from enterprise.signals import selections, parameter, white_signals

ephemeris = 'DE438'
parfiles = sorted(glob.glob('./B1937*par'))
timfiles = sorted(glob.glob('./B1937*tim'))
for p, t in zip(parfiles, timfiles):
    psr = Pulsar(p, t, ephem=ephemeris, timing_package='pint')

selection = selections.Selection(selections.by_backend)

ecorr = parameter.Constant()
ec = white_signals.EcorrKernelNoise(log10_ecorr=ecorr, selection=selection)

ec(psr)
vallis commented 2 years ago

OK, I found the problem. It's potentially dangerous. In commit 1d4975, I dropped a critical indexing step at line 266 in pulsar.py, which uses the toas sorting index to sort the backend_flags array. This resulted in scrambled selection masks when loading B1937 with PINT (but not with tempo2, which must do some sorting of its own). I just put in a PR to fix this on enterprise/master.

Hazboun6 commented 2 years ago

Nice catch @vallis !!

vallis commented 2 years ago

Thanks (although it was my fault originally). Fixed and merged in #297.