TTV fit in several instruments - issue when planet is not transiting in one of the instruments

[manustalport]

Hi,

I encountered an issue with juliet when trying to fit multi-planet TTVs using several photometric instruments. I don't know if it is an aspect missing from the code, or if it just reflects my bad use of it. Below is a brief description of the issue. Note that prior to having this problem, I had to slightly modify fit.py so to allow for multi-planet TTV fits, in the same way as reported in a former post.

Let's suppose I have 3 transiting planets and want to fit TTVs for all of them, using the T parametrization. Suppose also the photometric data are spread in 4 different instruments INS1, 2, 3, and 4. I would define the transit time priors (T_pi_instrument_n) for all transits of the 3 planets, and would not indicate t0 and P priors, given that they should be computed from the transit time posteriors. In one of the instruments (say INS2), one of the planets does not transit but the others do (say planet 3 is not transiting, for example). Running the code, fit.py returns an error asking me to provide t0 and P priors for planet 3, given that no T_p3_INS2_n prior was provided. In other words, the code expects to be given transit time priors for each planet one wants to fit TTVs on AND for each instrument. More technically, in the function 'generate_lc_model' of the fit.py file, the following condition will be satisfied for planet 3 on instrument INS2:

if not self.dictionary[instrument]['TTVs'][i]['status']:
    t0, P = parameter_values['t0_p' + str(i)], parameter_values['P_p' + str(i)]

From my understanding, this error should also take place in the case of a single transiting planet if it doesn't transit in one of the instruments. So I wanted to report this small problem (which I manually fixed in fit.py), and I would be happy to hear if there are ways to ask juliet to ignore certain planets in certain instruments. I guess I could also set t0 and P priors for the problematic planets, but wouldn't it bias the later estimation of t0 and P from the transit time posteriors?

Thanks!

Manu

[nespinoza]

Hi @manustalport,

I just released a new version of juliet. My hunch is that your issue should still be there, but can you double check?

If I'm understanding correctly, you want for a particular instrument on a particular night to not fit the transit. Can't you bypass this by making that particular transit data its own instrument, linking all instrumental parameters, and not defining that particular instrument's limb-darkening coefficients? If you don't set the limb darkening coefficients for an instrument, juliet assumes you just want the photometry.

Let me know if I misunderstood your question. A working example (even if simulated) would be great.

N.

[manustalport]

Hi Nestor,

Thank you for your reply. It already helps me to clarify certain things, as indeed in the case of 1 planet I could just omit the limb darkening coefficients on specific instruments.

However, I think this will not work for my toy example here above unfortunately. Because in my case, for the specific instrument where the issue arises, I want to fit for the transits of the other planets. So I cannot omit the limb darkening coefficients. Like, it is a situation of a system of planets observed in multiple photometric instruments, and it happens that one of the planets does not transit in one of the instruments (I don't want to ignore its transit, it is just not there). However in the fit.py code, I understood that for TTV fitting with Tparametrization, we HAVE TO indicate at least one transit time prior (T_pi_instrument_n) for each TTV planet and for EVERY instrument where limb darkening coefficients were indicated. Hence the issue I got, since juliet raises the error that no transit time was provided for a certain (planet, instrument) duo and therefore asks to provide (t0,P) priors for this planet instead. How would you specify to juliet that a certain planet is not transiting in a certain instrument, and that it's normal if no transit time prior is provided for this specific instrument and specific planet (but also knowing that the other planets transit in this instrument)? I wasn't so clear about this in my former post, sorry about that. Please let me know if this is still unclear.

I also tested the new version, and got indeed the same error message. I reproduced this error with a simplified simulated example, consisting of a system of 2 planets and 3 instruments. Planet b transits in all the instruments, while planet c only transits in instruments 1 and 3. Attached are the data files, and here below is the code which leads to the error message. SimulatedLC_2pl_ins1.txt SimulatedLC_2pl_ins2.txt SimulatedLC_2pl_ins3.txt

import os 
import numpy as np 
import juliet

t_INS1 = np.genfromtxt('SimulatedLC_2pl_ins1.txt', usecols=(0)) 
f_INS1 = np.genfromtxt('SimulatedLC_2pl_ins1.txt', usecols=(1)) 
ferr_INS1 = np.genfromtxt('SimulatedLC_2pl_ins1.txt', usecols=(2)) 

t_INS2 = np.genfromtxt('SimulatedLC_2pl_ins2.txt', usecols=(0)) 
f_INS2 = np.genfromtxt('SimulatedLC_2pl_ins2.txt', usecols=(1)) 
ferr_INS2 = np.genfromtxt('SimulatedLC_2pl_ins2.txt', usecols=(2)) 

t_INS3 = np.genfromtxt('SimulatedLC_2pl_ins3.txt', usecols=(0)) 
f_INS3 = np.genfromtxt('SimulatedLC_2pl_ins3.txt', usecols=(1)) 
ferr_INS3 = np.genfromtxt('SimulatedLC_2pl_ins3.txt', usecols=(2)) 

# Put data arrays into dictionaries so we can fit it with juliet:
times, fluxes, fluxes_error = {},{},{}

times['INS1'], fluxes['INS1'], fluxes_error['INS1'] = t_INS1,f_INS1,ferr_INS1
times['INS2'], fluxes['INS2'], fluxes_error['INS2'] = t_INS2,f_INS2,ferr_INS2
times['INS3'], fluxes['INS3'], fluxes_error['INS3'] = t_INS3,f_INS3,ferr_INS3

# Name of the parameters to be fit:
params = ['r1_p1','r2_p1','ecc_p1','omega_p1',\
          'r1_p2','r2_p2','ecc_p2','omega_p2',\
          'q1_INS1','q2_INS1','mdilution_INS1','mflux_INS1','sigma_w_INS1',\
          'q1_INS2','q2_INS2','mdilution_INS2','mflux_INS2','sigma_w_INS2',\
          'q1_INS3','q2_INS3','mdilution_INS3','mflux_INS3','sigma_w_INS3',\
            'rho']

# Distributions:
dists = ['uniform','uniform','fixed','fixed',\
         'uniform','uniform','fixed','fixed',\
         'uniform','uniform','fixed', 'normal', 'loguniform',\
         'uniform','uniform','fixed', 'normal', 'loguniform',\
         'uniform','uniform','fixed', 'normal', 'loguniform',\
            'loguniform'] 

# Hyperparameters
hyperps = [[0.,1.], [0.,1.], 0.0, 90.,\
           [0.,1.], [0.,1.], 0.0, 90.,\
           [0.,1.], [0.,1.], 1.0, [0.,0.1], [0.1, 1000.],\
           [0.,1.], [0.,1.], 1.0, [0.,0.1], [0.1, 1000.],\
           [0.,1.], [0.,1.], 1.0, [0.,0.1], [0.1, 1000.],\
                   [100., 10000.]]

##### Let's add the individual transit times  
# --- Planet b  -> P = 4.4 days 
# INS1
params = params + ['T_p1_INS1_0'] 
dists = dists + ['normal'] 
hyperps = hyperps + [[3.9,0.02]] 
# INS2
params = params + ['T_p1_INS2_1'] + ['T_p1_INS2_2']
dists = dists + ['normal'] + ['normal'] 
hyperps = hyperps + [[8.3,0.02]] + [[12.7,0.02]] 
# INS3
params = params + ['T_p1_INS3_5'] 
dists = dists + ['normal']  
hyperps = hyperps + [[25.9,0.02]] 

# --- Planet c  -> P = 7.0 days
# INS1 
params = params + ['T_p2_INS1_0'] 
dists = dists + ['normal'] 
hyperps = hyperps + [[0.45,0.02]] 
# INS2
# == No Transit of planet c == 
# INS3
params = params + ['T_p2_INS3_4'] 
dists = dists + ['normal'] 
hyperps = hyperps + [[28.45,0.02]] 

# Build the prior dictionary with the above information:
priors = juliet.utils.generate_priors(params,dists,hyperps)

dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, yerr_lc = fluxes_error, \
                      out_folder = 'TTV_2pl_3ins_test')

results = dataset.fit(sampler='dynesty', nthreads = 4, n_live_points=500)

Thank you! Manu

[nespinoza]

Hi @manustalport,

Thanks for providing this example! I'll have a look at this in a few weeks.

N.

[melissa-hobson]

Hi @manustalport , I am having the same issue, and you say you had a manual fix, if I understood your first post correctly? Could you share it?

[manustalport]

Hi @melissa-hobson

I'm happy to share with you the temporary modified version of fit.py that attempts to bypass the issue reported in this post. However, I recommend you to do a sanity check on a well-known system prior to using it. The file is attached below: fit_modif_AllTTV.py.zip

It contains two modifications compared to the original file: 1) It fixes the indexing problem for multi-planet TTV fits (as reported in a former post). Lines 646-650 become:

if pri[0:2] == 'dt' or pri[0:2] == 'T_':
    planet_number, instrument, ntransit = pri.split('_')[1:]
    if pri[0:2] == 'T_':
        dictionary[inames[i]]['TTVs'][int(planet_number[1:])]['parametrization'] = 'T'

2) To fix my issue reported above, I simply set the TTV ON by default. Line 622 becomes:

dictionary[inames[i]]['TTVs'][pi]['status'] = True

Please note that as a consequence, the program fits for TTV on all the planets. It won't work if you set linear ephemerides priors on one of the planets. That's why it is only temporary fix. The above dictionary[inames[i]]['TTVs'][pi]['status'] (loop on all instruments and planets) should be set to True only if in a multi-planet system, a planet is missing from an instrument AND if one wants to fit TTV on all planets. Alternatively, you could manually set to True the TTV status of the problematic (planet, instrument) duo. It is a "system-specific" fix in both cases. Anyway, I hope this helps!

[melissa-hobson]

Hi Manu, thanks, that looks like it should work for my case!

[nespinoza]

Solved in #110