nespinoza / juliet

A versatile modelling tool for transiting and non-transiting (single and multiple) exoplanetary systems
MIT License
53 stars 31 forks source link

KeyError: 'a_p1' #30

Closed TakuNishiumi closed 4 years ago

TakuNishiumi commented 4 years ago

Hello, again I executed below program. I used a_p1 instead of rho, then I got this error. Do you have any ideas to solve this?

Program:

import juliet
times, fluxes, fluxes_error = {},{},{}

# First define the priors:
priors = {}

# Same priors as for the transit-only fit, but we now add the GP priors:
params = ['P_p1','t0_p1','r1_p1','r2_p1','q1_M1g','q2_M1g','ecc_p1','omega_p1', 'a_p1',\
          'mdilution_M1g', 'mflux_M1g', 'sigma_w_M1g', \
          'GP_sigma_M1g', 'GP_rho_M1g']

dists = ['fixed','normal','uniform','uniform','uniform','uniform','fixed','fixed','fixed', \
        'fixed', 'normal', 'loguniform', \
         'loguniform', 'loguniform']

hyperps = [p, [t0,0.1], [0.,1], [0.,1.], [0., 1.], [0., 1.], e, 0.0, a, \
          1.0, [0.,0.1], [0.1, 1000.], \
           [1e-6, 1e6], [1e-3, 1e3]]

# Populate the priors dictionary:
for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

times['M1g'], fluxes['M1g'], fluxes_error['M1g'] = ts[0], rands+flux, df["g_err"].values

dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, \
                      yerr_lc = fluxes_error, GP_regressors_lc = times, out_folder = '55Cnc_GP_juliet3', verbose = True)

results = dataset.fit(dynesty_nthreads = 6)

Error:

     Transit fit detected for instrument  M1g
     >> ecc,omega parametrization detected for lc planet p1
     >> (b,p) parametrization detected for lc planet p1

PyMultinest installation not detected. Forcing dynesty as the sampler.

74582it [22:36, 34.21it/s, batch: 3 | bound: 374 | nc: 25 | ncall: 1316971 | eff(%):  5.663 | loglstar: 5873.061 < 5880.250 < 5879.175 | logz: 5757.302 +/-    nan | stop:  0.985]             

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-189-b2729c2d4e7e> in <module>
     27                       yerr_lc = fluxes_error, GP_regressors_lc = times, out_folder = '55Cnc_GP_juliet3', verbose = True)
     28 
---> 29 results = dataset.fit(dynesty_nthreads = 6)

~/anaconda3/envs/py35/lib/python3.6/site-packages/juliet/fit.py in fit(self, use_dynesty, dynamic, dynesty_bound, dynesty_sample, dynesty_nthreads, n_live_points, ecclim, delta_z_lim, pl, pu)
    716         return fit(self, use_dynesty = use_dynesty, dynamic = dynamic, dynesty_bound = dynesty_bound, dynesty_sample = dynesty_sample, \
    717                    dynesty_nthreads = dynesty_nthreads, n_live_points = n_live_points, ecclim = ecclim, delta_z_lim = delta_z_lim, \
--> 718                    pl = pl, pu = pu)
    719 
    720     def __init__(self,priors = None, input_folder = None, t_lc = None, y_lc = None, yerr_lc = None, \

~/anaconda3/envs/py35/lib/python3.6/site-packages/juliet/fit.py in __init__(self, data, use_dynesty, dynamic, dynesty_bound, dynesty_sample, dynesty_nthreads, n_live_points, ecclim, delta_z_lim, pl, pu, ta)
   1321             if not os.path.exists(self.out_folder+'posteriors.dat'):
   1322                 outpp = open(self.out_folder+'posteriors.dat','w')
-> 1323                 writepp(outpp,out, data.priors)
   1324 
   1325         # Save all results (posteriors) to the self.results object:

~/anaconda3/envs/py35/lib/python3.6/site-packages/juliet/utils.py in writepp(fout, posteriors, priors)
    528                     a = ((posteriors['posterior_samples']['rho']*G*((priors['P_'+planet]['hyperparameters']*24.*3600.)**2))/(3.*np.pi))**(1./3.)
    529             else:
--> 530                 a = posteriors['posterior_samples']['a_'+planet]
    531             inc_inv_factor = (b/a)*ecc_factor
    532             inc = np.arccos(inc_inv_factor)*180./np.pi

KeyError: 'a_p1'
nespinoza commented 4 years ago

Hi,

That's very weird, because it seems the fit finished to run, but the problem is happening when trying to write the posteriors.dat file. Do you have more than one pickle file in the output folder?

Best, N.

TakuNishiumi commented 4 years ago

My output folder have these files.

_dynesty_NS_posteriors.pkl  lc.dat          priors.dat
GP_lc_regressors.dat        posteriors.dat
TakuNishiumi commented 4 years ago

I realized that at the first execution, the below error was returned.

     Transit fit detected for instrument  M1g
     >> ecc,omega parametrization detected for lc planet p1
     >> (b,p) parametrization detected for lc planet p1
PyMultinest installation not detected. Forcing dynesty as the sampler.

501it [00:00, 1281.76it/s, batch: 0 | bound: 0 | nc: 1 | ncall: 501 | eff(%): 100.000 | loglstar:   -inf <   -inf <    inf | logz:   -inf +/-    nan | dlogz:  0.000 >  0.010]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-193-04d4a9319b6a> in <module>
     28                       yerr_lc = fluxes_error, GP_regressors_lc = times, out_folder = '55Cnc_GP_juliet4', verbose = True)
     29 
---> 30 results = dataset.fit(dynesty_nthreads = 6)

~/anaconda3/envs/py35/lib/python3.6/site-packages/juliet/fit.py in fit(self, use_dynesty, dynamic, dynesty_bound, dynesty_sample, dynesty_nthreads, n_live_points, ecclim, delta_z_lim, pl, pu)
    716         return fit(self, use_dynesty = use_dynesty, dynamic = dynamic, dynesty_bound = dynesty_bound, dynesty_sample = dynesty_sample, \
    717                    dynesty_nthreads = dynesty_nthreads, n_live_points = n_live_points, ecclim = ecclim, delta_z_lim = delta_z_lim, \
--> 718                    pl = pl, pu = pu)
    719 
    720     def __init__(self,priors = None, input_folder = None, t_lc = None, y_lc = None, yerr_lc = None, \

~/anaconda3/envs/py35/lib/python3.6/site-packages/juliet/fit.py in __init__(self, data, use_dynesty, dynamic, dynesty_bound, dynesty_sample, dynesty_nthreads, n_live_points, ecclim, delta_z_lim, pl, pu, ta)
   1181                         sampler = dynesty.DynamicNestedSampler(self.loglike, self.prior, self.data.nparams, nlive = self.n_live_points, \
   1182                                                               bound = self.dynesty_bound, sample = self.dynesty_sample, pool=executor, queue_size=nthreads)
-> 1183                         sampler.run_nested()
   1184                         results = sampler.results
   1185                 out['dynesty_output'] = results

~/anaconda3/envs/py35/lib/python3.6/site-packages/dynesty/dynamicsampler.py in run_nested(self, nlive_init, maxiter_init, maxcall_init, dlogz_init, logl_max_init, n_effective_init, nlive_batch, wt_function, wt_kwargs, maxiter_batch, maxcall_batch, maxiter, maxcall, maxbatch, n_effective, stop_function, stop_kwargs, use_stop, save_bounds, print_progress, print_func, live_points)
   1668                                               print_progress=print_progress,
   1669                                               print_func=print_func,
-> 1670                                               stop_val=stop_val)
   1671                     ncall, niter, logl_bounds, results = passback
   1672                 elif logl_bounds[1] != np.inf:

~/anaconda3/envs/py35/lib/python3.6/site-packages/dynesty/dynamicsampler.py in add_batch(self, nlive, wt_function, wt_kwargs, maxiter, maxcall, logl_bounds, save_bounds, print_progress, print_func, stop_val)
   1768                 lnz, lnzerr = res.logz[-1], res.logzerr[-1]
   1769                 if logl_bounds is None:
-> 1770                     logl_bounds = wt_function(res, wt_kwargs)
   1771                 for results in self.sample_batch(nlive_new=nlive,
   1772                                                  logl_bounds=logl_bounds,

~/anaconda3/envs/py35/lib/python3.6/site-packages/dynesty/dynamicsampler.py in weight_function(results, args, return_weights)
    149     nsamps = len(logz)
    150     bounds = np.arange(nsamps)[weight > maxfrac * max(weight)]
--> 151     bounds = (min(bounds) - lpad, min(max(bounds) + lpad, nsamps - 1))
    152     if bounds[0] < 0:
    153         logl_min = -np.inf

ValueError: min() arg is an empty sequence

and I re-executed the same cell, the program started to sampling and I got the "KEY Error" I described.

Sincerely, Taku

nespinoza commented 4 years ago

Interesting. This seems to be a dynesty issue really, but perhaps we can do something about it. If instead of doing results = dataset.fit(dynesty_nthreads = 6) you try results = dataset.fit(use_dynesty = True, dynesty_nthreads = 6) do you get the same result in the first run?

TakuNishiumi commented 4 years ago

Thank you. I tried it and got the same ValueError: min() arg is an empty sequence in the first run.

nespinoza commented 4 years ago

Hard to debug here, really. Could you send me an e-mail with the dataset + script so I can run it on my machine to check (nespinoza at stsci.edu)?

Thanks!

TakuNishiumi commented 4 years ago

Yes, I will send you. Thank you.

nespinoza commented 4 years ago

Hi @TakuNishiumi,

I'm closing this issue, as it turned to not be a juliet problem. For readers having questions about this: this was related to some of the parameters being ingested as fixed to juliet having unphysical values, so the code didn't produce any samples. Perhaps I'll set a warning about such fixed parameters in the future --- but this is very low priority right now. Will create an "enhancement" issue for future reference about this.

Best, Néstor

TakuNishiumi commented 4 years ago

Dear @nespinoza

Thank you very much for your help. @TakuNishiumi