rodluger / starry

Tools for mapping stars and planets.
https://starry.readthedocs.io
MIT License
138 stars 32 forks source link

"Hot jupiter phase curve" notebook with non-fixed planet radius #280

Closed LucaNap closed 3 years ago

LucaNap commented 3 years ago

OS: Windows 10, with latest version of Starry on Anaconda

I have encountered this issue while trying to fit my own data so I tried to replicate it into the Hot Jupiter notebook in order to make it simple to debug: basically I am trying to find the planet radius injected into the system (r=0.1) by using in the pymc3 model:

p_radius = pm.Normal("p_radius", mu=0.1, sd=0.02)

and

b = starry.Secondary( starry.Map(ydeg=1, udeg=0, amp=10 ** log_amp, inc=90.0, obl=0.0), m=0.0, r=p_radius, prot=1.0, porb=1.0, )

and the notebook runs fine up until I try to evaluate the posteriors over the model parameters. At this point, for some reason (which might be more related to pymc3 rather than starry itself), the chains get stuck at 0% and the following error is returned:

BrokenPipeError Traceback (most recent call last) D:\Anaconda3\lib\multiprocessing\connection.py in _recv_bytes(self, maxsize) 300 try: --> 301 ov, err = _winapi.ReadFile(self._handle, bsize, 302 overlapped=True)

BrokenPipeError: [WinError 109] The pipe has been ended

During handling of the above exception, another exception occurred:

EOFError Traceback (most recent call last)

in 1 with model: ----> 2 trace = pmx.sample( 3 tune=250, 4 draws=500, 5 start=map_soln, ~\AppData\Roaming\Python\Python38\site-packages\pymc3_ext\sampling\sampling.py in sample(draws, tune, model, step_kwargs, warmup_window, adapt_window, cooldown_window, initial_accept, target_accept, gamma, k, t0, parameter_groups, adapt_type, recompute_interval, regularization_steps, regularization_variance, **kwargs) 95 ) 96 ---> 97 return pm.sample(draws=draws, tune=tune, model=model, step=step, **kwargs) D:\Anaconda3\lib\site-packages\pymc3\sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs) 557 _print_step_hierarchy(step) 558 try: --> 559 trace = _mp_sample(**sample_args, **parallel_args) 560 except pickle.PickleError: 561 _log.warning("Could not pickle model, sampling singlethreaded.") D:\Anaconda3\lib\site-packages\pymc3\sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs) 1475 try: 1476 with sampler: -> 1477 for draw in sampler: 1478 trace = traces[draw.chain - chain] 1479 if trace.supports_sampler_stats and draw.stats is not None: D:\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py in __iter__(self) 477 478 while self._active: --> 479 draw = ProcessAdapter.recv_draw(self._active) 480 proc, is_last, draw, tuning, stats, warns = draw 481 self._total_draws += 1 D:\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py in recv_draw(processes, timeout) 349 idxs = {id(proc._msg_pipe): proc for proc in processes} 350 proc = idxs[id(ready[0])] --> 351 msg = ready[0].recv() 352 353 if msg[0] == "error": D:\Anaconda3\lib\multiprocessing\connection.py in recv(self) 248 self._check_closed() 249 self._check_readable() --> 250 buf = self._recv_bytes() 251 return _ForkingPickler.loads(buf.getbuffer()) 252 D:\Anaconda3\lib\multiprocessing\connection.py in _recv_bytes(self, maxsize) 319 except OSError as e: 320 if e.winerror == _winapi.ERROR_BROKEN_PIPE: --> 321 raise EOFError 322 else: 323 raise EOFError:

The same goes for the star radius, masses and so on. Am I doing anything wrong?

rodluger commented 3 years ago

This is probably not a starry issue. See here for a similar issue using pymc3: https://github.com/pymc-devs/pymc3/issues/3388. It actually seems to be specifically an issue for multiprocessing in Jupyter notebooks on Windows: https://github.com/pymc-devs/pymc3/issues/3388#issuecomment-468636046

One of their suggestions is to use cores=1 to disable multiprocessing.

LucaNap commented 3 years ago

This is probably not a starry issue. See here for a similar issue using pymc3: pymc-devs/pymc3#3388. It actually seems to be specifically an issue for multiprocessing in Jupyter notebooks on Windows: pymc-devs/pymc3#3388 (comment)

One of their suggestions is to use cores=1 to disable multiprocessing.

Spot on! Also, apparently sometimes it's useful to use the option pickle_backend='dill'. Just for the record, I am having a couple of problems with the planetary radius sometimes going negative in the sampling (thus interrupting the process), but the issue is probably going to be fixed by:

p_radius = pm.Bound(pm.Normal, lower=0.01, upper=0.2)("p_radius", mu=0.1, sd=0.05)

Thank you

rodluger commented 3 years ago

Yes, that will probably do it. Just keep in mind that you might have to explicitly pass a testval when using pm.Bound:

p_radius = pm.Bound(pm.Normal, lower=0.01, upper=0.2)("p_radius", mu=0.1, sd=0.05, testval=0.1)

I remember an issue (which may have been fixed) that it defaults to zero, which could cause trouble. But I could be wrong!