Iterations for emcee

[AstroJLin]

Hi Kyle,

1. Since I need a lot of iterations, I try to use b.export_solver and b.import_solution on my laptop. However, it seems does not work and the phoebe show no such file or directory.

My code as follow:
b.save('inverse_paper_examples_before_emcee.bundle')
b.export_solver('inverse_paper_examples_run_emcee.py',
solver='emcee_solver',
niters=2000, progress_every_niters=100,
nwalkers=16,
solution='emcee_sol',
log_level='warning',
pause=True)
import matplotlib.pyplot as plt
plt.rc('font', family='serif', size=14, serif='STIXGeneral')
plt.rc('mathtext', fontset='stix')
import phoebe
import numpy as np
logger = phoebe.logger('error')
b = phoebe.load('inverse_paper_examples_before_emcee.bundle')
b.import_solution('inverse_paper_examples_run_emcee.py.out', solution='emcee_sol')

2. Alternatively, I try to use b.run_solver('emcee_solver', niters=3000, nwalkers=16, solution='emcee_sol') on my laptop, it seems does not work after about 10 iterations.

How should I solve this problem?

Many Thanks

[kecnry]

b.export_solver(...) will create a python script named "inverse_paper_examples_run_ecmee.py" which you then need to run (on another machine, if you'd like) in order to create the "inverse_paper_examples_run_emcee.py.out" file. My suggestion would be to run on a cluster within MPI or multiprocessing and set the progress_every_niters parameter to some reasonable number (say 10) so that you can monitor the progress while its running and don't risk losing everything if the job doesn't complete or gets stuck.

This tutorial might be of some help as well: http://phoebe-project.org/docs/2.3/tutorials/export_solver

[kecnry]

Adding the distribution only defines the distribution objects themselves. You then need to reference them in a parameter (such as init_from and/or priors in the case of mcmc). init_from is required to be set, as it determines how to draw the walkers for the first iteration. priors is optional, and affects the cost function itself. It looks here like you have not actually used the defined distribution for either of these "uses".

[AstroJLin]

Hi Kyle,
Thank you for your reply again.

Actually, I have set a priors after I define the distribution. However, it still throw a errors that is no valid distribution in init_from as shown Figure.

You means that the init_from must be required to be set. Otherwise it won't work if I only set a priors? For example, Here I need to set the init_from on the effective temperature (T) and the radius (R) at first If I want to adopt Gaussian priors on the effective temperature (T) and the radius (R)?

Is the code like this?

b.add_solver('sampler.emcee', solver='emcee_solver')
b.set_value('compute', solver='emcee_solver', value='fastcompute')

b.add_distribution({'teff@primary': phoebe.gaussian(4016,35), 'requiv@primary':phoebe.gaussian(0.68,0.01)}, distribution='ball_around_guess')
b.set_value('init_from', 'ball_around_guess')

b.add_distribution({'teff@primary':phoebe.gaussian(4016,35),'requiv@primary':phoebe.gaussian(0.68,0.01)},
distribution='priors_from_external_source')
b.set_value('priors', 'priors_from_external_source')

b.set_value('nwalkers', solver='emcee_solver', value=12)
b.set_value('niters', solver='emcee_solver', value=250)
b.run_solver('emcee_solver', solution='emcee_sol')

[kecnry]

Yes, init_from is required, priors is optional. You were getting the error because you left init_from blank. You can set init_from and priors to the same distribution label if you would like, but they can also be different (and priors can include distributions on constrained parameters that are not included in init_from).

In your case here, there is no need to create two separate distributions if you want to initialize from the same distribution that you're using for your priors, but your code should work as-is.

[AstroJLin]

Hi Kyle,

Thank you very much your patience.

If I set init_from and priors to the same distribution label, then those parameters from init_from and priors will be shown in the posterior distributions. However, If they has different distribution label, then those parameters from priors will not be shown in the posterior distributions though priors is affects the cost function itself. Is that right?

[kecnry]

That's correct - you'll get posteriors from all parameters in init_from, whether they're in priors or not. And if you put a prior on a parameter not in init_from, it will not be included as MCMC is not sampling over that parameter.

You can always propagate the posterior distributions to see how those posteriors affect the constrained parameters with priors.

[AstroJLin]

Hi Kyle,

Thank you very much for your reply again.

Here, I have two questions.

1、When I try to run a script on a cluster, it does not work. However, when I run this script on my laptop, it can work. The error is as follow

I don't know what caused this error.

2、I have the CSS light curve, which is in V magnitudes. I want to know whether passband='Johnson:V' [b.add_dataset('lc', times=times, fluxes=fluxes, sigmas=sigmas_lc,passband='Johnson:V')] should be given. If I do not set passband='Johnson:V', how much will this affect the results of the model? In addition, to make the model match the data, should I try to set the value of the pblum [b.set_value('pblum', component='primary', value=0.12*np.pi)] or use normalized flux? In other words, which of these two methods will be more conducive to get the correct parameters through the model.

[kecnry]

1. That is a known bug that I'm working to fix and will release once I have a fix. In the meantime, downgrading astropy on the cluster should avoid the bug. Be on the lookout for the next bugfix release (and I'll try to remember to update back here once its released).

2. You'll need to convert to fluxes (phoebe does not currently support magnitudes). If you don't set a passband, it will default to Johnson:V. The contribution of the choice of passband to the light curve is quite situational dependent, so if you're unsure of the passband set, you could try switching between different close passbands to see how much affect it has on the synthetic model for your particular system. As far as pblum, if you have calibrated absolute fluxes, then using pblum_mode='absolute' and providing/fitting for the luminosity and distance gives you the best information. Otherwise, your fluxes and passband luminosities will be arbitrarily scaled.

Hope this helps!

[kecnry]

I believe the bug you ran into should be fixed by 2.3.51 and 2.3.52. Please update and if you do still run into something similar, file a bug report and let me know the versions of phoebe, astropy, and numpy that you're using so I can better reproduce and fix it. Thanks!

[AstroJLin]

Hi Kyle,

Thank you very much for your reply and help.

I have updated the versions of phoebe to 2.3.51 and the bug has disappeared. However, it runs as expensive as my laptop when I run the script of inverse_paper_examples_run_ecmee.py on a cluster. Therefore, I want to ask you how to set number of thread or process to run faster on a cluster.

As far as pblum, if I use pblum_mode=dataset-scaled, will it have a significant effect on some parameters (such as the radii and temperatures of the star and so on) derived from the model.

[kecnry]

PHOEBE itself should pick up on the number of processors sent via mpirun/mpiexec or available to multiprocessing and distribute the job intelligently (either by-time within the forward model or per-walker within mcmc). If you're using a scheduler like slurm, just provide the number of processors in your submission script and run the script within mpirun/mpiexec to use MPI, or without to use multiprocessing.

Using dataset-scaled can cause walkers to "freeze" within MCMC, so I would probably suggest switching to another mode and "marginalizing" over the luminosity once you get to MCMC. But datset-scaled can be quite convenient while optimizing or hand-tweaking the model. If you either use dataset-scaled OR non-physical luminosities (to match normalized fluxes, for example), then its just important to remember that your luminosities, distances, and fluxes will all have arbitrary units.

[AstroJLin]

Hi Kyle,

Thank you very much for your patience.

You means that the luminosity, distances and fluxes will all have arbitrary units when I use dataset-scaled Or non-physical luminosity.

Actually, I want to know whether these parameters (the binary mass ratio, the orbital inclination, radius and temperature of the star) produces inaccuracy when I use the model to synthetic observable with pblum_mode=dataset-scaled.

[AstroJLin]

Hi Kyle,

Actually, my last question is that normalized light curve must be required for the phoebe. In other words, If I use dataset-scaled but not a normalized light curve, is it have an effect on the optimization parameters?

Thanks

[kecnry]

No, normalized fluxes are not necessary. dataset-scaled uses a quick least-squares run to determine a scale factor for each individual forward-model, but does not care what that scale factor is or means. Once you start fitting for pblum, the input fluxes will then affect the fitted value for pblum... but if they're not calibrated absolute fluxes, then the value of pblum is arbitrary anyways, so you won't care about the actual physical value. Other parameters should remain unaffected either way.

[AstroJLin]

Hi Kyle,
Thank you very much for your reply.

I have folded the light curve into a phase format outside of phoebe. Now, I fixed an orbital period value and tried to fit the folded light curve with a model. However, the result of the fitting is very strange. (

).
Could you give me some suggestions.

Thanks

[kecnry]

PHOEBE takes times as input, but does not support phases (intentionally). If you only have access to phased data, then I would suggest trying to unphase it to a single orbital cycle with the best estimate for the orbital period available (and if the orbital period changes during the analysis, you may want to update that mapping and/or be transparent about the assumptions in any publications). If you do have time-series data, I would suggest leaving the times as-is and set compute_phases in the dataset if you're trying to reduce the number of time-points necessary for each forward model. This will then interpolate linearly in phase space automatically when computing the cost function.

[AstroJLin]

Hi Kyle,
Thank you very much for your patience.

Actually, you mean that the phase conversion is performed within the phoebe. Here, as long as I input time-series data and set an orbital period and give a value of t0_supconj in the model, the phoebe will perform the phase conversion. The code can be write as follow:
b = phoebe.default_binary()
b.set_value('t0_supconj',2457853.185245582)
b.set_value('period', component='binary', value=0.5556686*u.d)
b.set_value('ecc', 0.02)
b.set_value('incl@binary', 75)
b.set_value('q', 0.5)
b.set_value('teff@primary', 4816)
b.add_dataset(kind='lc', times=times, fluxes=fluxes, sigmas=sigmas_lc,passband='Johnson:V')
b.run_compute()

Is that right?

[AstroJLin]

In addition, when I try to use my model to synthetic observable with b.set_value('pblum', component='primary', value=0.12*np.pi) rather than pblum_mode=dataset-scaled, the orbital inclination is very small compared with pblum_mode=dataset-scaled.

I don't know what caused this difference?

[kecnry]

If you're not using dataset-scaled, you'll want to include pblum as a free parameter. Here the inclination is likely artificially low to try to increase the out-of-eclipse fluxes to match the observations (its using inclination to account for the scale since you're not giving it the freedom to with luminosities)

[AstroJLin]

Hi Kyle,

Thank you very much for your suggestion, I will do with your responses!