Is there an phoebe example of fitting the light curve affected by star spots with mcmc method?

[dingxu6207]

I want to fit a light curve affected by star spots, but I do not know how to do it. Is there an example? Many thanks!

[kecnry]

Unfortunately there is no example for that exact scenario, but you should be able to piece together what you need from a few different examples:

• binary with spots shows how to define spots parameters and attach them to individual stars
• emcee sampler tutorial shows how to initialize an emcee run. You'll just reference the spot parameters that you want to include in your initializing (and optionally priors) distributions. The distributions tutorial and solvers tutorial will also be helpful if you haven't read them already.

Keep in mind that spots are terribly degenerate, so even if you get decent looking results, I would caution against claiming it as a physical or constrained solution.

[dingxu6207]

Thank you very much, I have solved it.

[aliyyanurr]

Hello, I am also exploring an observational binary with spots right now and am a bit clueless. What is the 'spot feature'? Could the spots be optimized with NMS by calling 'spot' to the fit_parameters=[] when adding NMS solvers? I randomly tried this and run the optimizer outside (export_compute, remote running using nohup), but when the status is complete, there is no '.py.out' file but the nohup.out file size is large (~200GB)

Or, when optimizing the spots, should I call each of the parameter twigs (like "colat@spot01@feature") in the fit_parameter=[] (and by using MCMC, when initializing the MCMC and creating the distribution of each spot parameter)?

[kecnry]

Yes, you'll need to reference the individual parameter names that you want to adjust in either fit_parameters (for optimizers) or through the initializing distributions (for samplers). If that gives you any problems, feel free to send a minimal script or output errors and we'll try to help!

[aliyyanurr]

Thank you. I tried adding spot individual parameters (relteff@spot01, relteff@spot01) with another parameters to be optimized using NMS in this script:

b.add_compute('ellc', compute='fastcompute', overwrite=True)
b.add_solver('optimizer.nelder_mead', solver = 'nm_solver_spot', fit_parameters=['teffratio', 'incl@binary', 'q', 'per0', 'relteff@spot01@feature', 'relteff@spot02@feature'], compute='fastcompute', overwrite=True)

and export the solver to be computed in my univ computational facility
b.export_solver(script_fname = 'nelder_mead_spot.py', solver = 'nm_solver_spot', maxiter=1000, solution='nm_sol_spot')

In terminal, I used nohup command:
nohup python3 nelder_mead_spot.py

I ever used export_solver and nohup in my previous project and everything went nicely, the program ran well, there was '.py.out' file output, and I can use this by calling b.import_solution.

But when optimizing the spots just now (it takes ~1 hour), the nohup status is completed, but there is no '.py.out' file. Besides, there is 'nohup.out' with the size of 25 GB. There is no output of the optimization process that I could use. I wonder if the computation process takes so too much space so it stopped before it actually finished.

I'll also send the bundle here
FFVul_optimize_spot_1.bundle.zip .

[kecnry]

I'm guessing that something is failing and so it never gets to write the output. My suggestion would be to try running it in the original thread (use run_solver directly without export_solver or nohup) so that its easier to see logger warnings/errors. If that doesn't work or you don't know what the error means, let me know and I'll try to reproduce on my machine.

[aliyyanurr]

I tried to run the b.run_solver before, but its behavior is different with another run_solver I did. There is no progress/warning/error printed (i set the phoebe_logger = error), the progress bar also doesn't change. In terminal where i opened the jupyter notebook, there are many values and numbers printed and change rapidly (my another run_solver didn't do that). It took a long time, I don't know how long it takes, so I interrupted the run_solver...

[aliyyanurr]

Hello, I'd like to inform you and ask about my previous progress and questions.

I tried to b.run_solver the spot NMS optimizer on my jupyter notebook with this script:

#defining spot
b.add_feature('spot',
              component='primary',
              feature='spot01',
              relteff=0.8, radius=10, colat=75, long=25)

b.add_feature('spot',
              component='secondary',
              feature='spot02',
              relteff=1.3, radius=20, colat=105, long=10)

#adding solver
b.add_solver('optimizer.nelder_mead', solver = 'nm_solver_spot01',
             fit_parameters=['colat@spot01@feature', 'long@spot01@feature'],
             compute='fastcompute', overwrite=True)

#run solver
b.run_solver(solver='nm_solver_spot01', maxiter=200, solution='nm_sol_spot01')

but it takes a long time (more than 12 hours; I'm running PHOEBE at my university computational service server), so I always halt the process before it finishes. I tried the b.export_solver with the same code above but the b.run_solver is changed to:

b.export_solver('25jan22_nmsolver_spot01.py', 
                solver='nm_solver_spot01',
                maxiter=2000,
                solution='nm_solver_spot01_sol')

I vary the max iter number from 1000 to 10000 (which differs at the computing time, around 2 to 4 days). The b.export_solver seems to work since it returns *.py.out file, but when I tried to import and see the solution with this code:

b.import_solution(fname = 'nelder_mead_spot_p01.py.out', solution = 'nm_solver_spot_p01_sol')
print(b.get_solution('nm_solver_spot_p01_sol'))

the output of all varying num of max iter is:

...
R  message@nm_solver_spot01_p0...: Maximum number of iterations has been exceeded.
R  nfev@nm_solver_spot01_p01@n...: 39999
R  niter@nm_solver_spot01_p01@...: 10000
R  success@nm_solver_spot01_p0...: False
...

My first question is: did I assign the parameters and optimization settings wrong? Since it seems like all the optimization attempts failed.

Then I tried to use MCMC to sample the spots (like this forum's questions) since dingxu6207 said he/she has solved it. Here's how I prepare the MCMC run:

#adding prior dist
b.add_distribution({'colat@spot01@feature': phoebe.gaussian_around(1),
                   'long@spot01@feature': phoebe.gaussian_around(1),
                   'radius@spot01@feature': phoebe.gaussian_around(1),
                   'relteff@spot01@feature': phoebe.gaussian_around(0.3),
                   'colat@spot02@feature': phoebe.gaussian_around(1),
                   'long@spot02@feature': phoebe.gaussian_around(1),
                   'radius@spot02@feature': phoebe.gaussian_around(1),
                   'relteff@spot02@feature': phoebe.gaussian_around(0.2)},
                  distribution = 'spot_prior', overwrite_all=True)

#adding mcmc solver
b.add_solver('sampler.emcee',
             init_from='spot_prior',
             compute='fastcompute',
             solver='mcmc_spot')

#saving bundle and export solver
b.save('7feb22_mcmc_spot.bundle')
b.export_solver('7feb22_mcmc_spot.py', 
                solver='mcmc_spot',
                niters=5000, progress_every_niters=100, 
                nwalkers=16,
                solution='mcmc_spot_sol',
                log_level='warning',
                pause=True)

It takes around 2-3 days and seems to work since it returns *.py.out file. But when I import the solution
b.import_solution('7feb22_mcmc_spot.py.out', solution='spot_mcmc_1') and want to make the lnprobability plot _ = b.plot(solution='spot_mcmc_1', style='lnprobability', burnin=0, thin=1, show=True), it returns an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-65c7e1699869> in <module>
----> 1 _ = b.plot(solution='spot_mcmc_1', style='lnprobability', burnin=0, thin=1, show=True)

/opt/anaconda3/lib/python3.7/site-packages/phoebe/parameters/parameters.py in plot(self, twig, **kwargs)
   5662                                           subplot_grid=subplot_grid,
   5663                                           animate_callback=animate_callback,
-> 5664                                           **kwargs)
   5665             except Exception as err:
   5666                 self.clf()

/opt/anaconda3/lib/python3.7/site-packages/phoebe/parameters/parameters.py in _show_or_save(self, save, show, animate, fig, draw_sidebars, draw_title, tight_layout, subplot_grid, **kwargs)
   5769             afig = self.gcf()
   5770             if not len(afig.axes):
-> 5771                 raise ValueError("Nothing could be found to plot.  Check all arguments.")
   5772 
   5773 

ValueError: Nothing could be found to plot.  Check all arguments.

My second question is the same as the first question: did I assign the parameters and MCMC sampler settings wrong? I don't understand the error; why is it nothing could be found to plot?

I'll attach the solution py.out file here in ZIP. Thank you!
FF Vul.zip

[kecnry]

For your optimization run, the returned solution says that it failed with a message from the optimizer that 'Maximum number of iterations has been exceeded.'. This means that even in 10000 iterations, it failed to achieve the set tolerances. There's a chance the original model fails, but you're not exposing the original lnprobability so I can't easily tell (more on this later for MCMC). But assuming that's not the case, you can increase the number of iterations further, but chances are it still might never fully converge. Alternatively, you can alter the tolerances or give more freedom to your forward model. If you're willing to add a bit of overhead (as it re-runs the forward model every N iterations and exports the file), you can also use progress_every_niters to ensure that you'll at least get the "best" result at the end, even if the tolerances aren't reached.

As for MCMC, it looks like you were exposing failed samples, which is great news! Try plotting those (with style='failed'). It looks like many of them are because the returned lnprobability was nan (this might actually be the same problem with the optimizer). So let's try to figure out why that is:

• Does your forward model run with the face values? What does b.calculate_lnlikelihood return if you pass that a model computed at the face values?
• Does it run if you sample across your init_from distributions (by passing your init_from distribution to sample_from in run_compute)?
• Do you have any priors that may be rejecting your initial samples? If often like making corner plots of both init_from and priors to make sure they aren't in disagreement.

[kittipongwangnok]

May I have the example of FF Vul light curve fitting? The previous zip file that you uploaded seems missing the data of LC_FFVul_B.txt. Thank you very much.