Inverse problem - WEB server results not reproducible in python

[mcregogarcia]

I have managed to get access to the web server for 2 days in the last six months. During those 2 days I managed to solve wit lcgeom solver a light curve and got some parameters for it.

The curves were not perfectly matching, but they were similar enough to start workign with the solution.

since the server has been unavailalbe for me ever since I tried to follow the invers problem tutorial and reproduce it in python, with no success.

the script runs without giving me any errors, but I do not get any decent LC from it. Is this something I can get help with?

The script is as follows:

#################SET-UP

import phoebe
from phoebe import u #units
import numpy as np
import matplotlib.pyplot as plt
from reportlab.lib.pagesizes import letter
from reportlab.pdfgen import canvas
#plt.rc('font', family='serif', size=14, serif='STIXGeneral')
#plt.rc('mathtext', fontset='stix')

logger = phoebe.logger()

#################LOAD LC

lcdata = np.loadtxt('data/LC3.data', delimiter='\t')

lctimes = lcdata[:,0]
fluxes = lcdata[:,1]
fsigmas = lcdata[:,2]

#################PLOT LC

plt.plot(lctimes,fluxes)
plt.show()

### UP TO HERE IT WORKS FINE

#################GENERATE DATA

b = phoebe.default_binary()
b.add_constraint('semidetached', 'primary')
b['requiv@constraint@primary']
b['requiv_max@constraint@primary']

b.set_value('latex_repr', component='binary', value='orb')
b.set_value('latex_repr', component='primary', value='1')
b.set_value('latex_repr', component='secondary', value='2')

#################SET PARAMETERS
#THESE ARE THE VALUES I GOT AS SOLUTION USING THE WEB INTERFACE.

b.set_value(qualifier='period', context='component', component='binary', value='6.45895 d')
b.set_value(qualifier='teff', context='component', component='primary', value='15700.0 K')
b.set_value(qualifier='teff', context='component', component='secondary', value='5000.0 K')
b.set_value(qualifier='ld_mode_bol', context='component', component='primary', value='manual')
b.set_value(qualifier='gravb_bol', context='component', component='primary', value='0.9')
b.set_value(qualifier='irrad_frac_refl_bol', context='component', component='primary', value='1')
b.set_value(qualifier='sma', context='component', component='binary', value='26.0 solRad')
b.set_value(qualifier='requiv', context='component', component='secondary', value='1.0 solRad')
b.set_value(qualifier='syncpar', context='component', component='primary', value='14.00')

#################ADD LC DATASETS

b.add_dataset('lc',
        compute_phases=phoebe.linspace(0,1,101),
                times=lctimes, 
                fluxes=fluxes, 
                sigmas=fsigmas, 
                dataset='lc01',
                overwrite=True)
b.set_value_all('ld_mode', 'manual')
b.set_value_all('atm', 'blackbody')
b.set_value(qualifier='ld_mode_bol', context='component', component='primary', value='manual')

#################COMPUTE MODEL

b.run_compute(irrad_method='none')

#################INITIALIZE BUNDLE

b = phoebe.default_binary()
b.add_dataset('lc', times=lctimes, fluxes=fluxes)
b.set_value_all('ld_mode', 'manual')
b.set_value_all('ld_mode_bol', 'manual')
b.set_value_all('atm', 'blackbody')
b.set_value('pblum_mode', 'dataset-scaled')

b.run_compute(model='default')
_ = b.plot(x='times', show=True)

#################ADOPT A SOLUTION METHOD LC ESTIMATOR

b.add_solver('estimator.lc_geometry', lc_datasets='lc01', solver='lcgeom')
print(b.filter(solver='lcgeom'))
b.run_solver('lcgeom', solution='lcgeom_solution', max_computations=15000, detach=True, overwrite=True)
print(b.filter(solution='lcgeom_solution'))
plt.plot(b.filter(solution='lcgeom_solution'), color='r', label='solution')
plt.plot(lctimes, fluxes, color='g', label= 'data')
plt.show()

Thanks in advance and apologies for the inconveniences.

[kecnry]

Apologies about the web server being down so often - we are working on a more permanent solution, but that may take a little while yet.

Without the data, its hard to reproduce and tell what is going on here. I do see that you run an lc geometry estimator, though, but do not actually adopt its results and don't re-run the forward model. You're also resetting the bundle (b) by calling phoebe.default_binary() again, meaning you lose all the parameters that you set from your previous work in the web UI. Try putting a b.plot(show=True) after your b.run_compute(irrad_method='none') line and see if that looks how you'd expect. If you then want to run an estimator or optimizer, just continue with that bundle instead of creating a new one so that you can build on top of existing results.

[mcregogarcia]

miguel.zip

Hi,

thanks for the reply. I have enclosed the bundle ( OPCION1.py), the script (miguelY.py) and the source date ( LC3.data)

when I fix it now I get it this>>>

#################SET-UP

import phoebe
from phoebe import u #units
import numpy as np
import matplotlib.pyplot as plt
from reportlab.lib.pagesizes import letter
from reportlab.pdfgen import canvas
#plt.rc('font', family='serif', size=14, serif='STIXGeneral')
#plt.rc('mathtext', fontset='stix')

logger = phoebe.logger()

#################LOAD LC

lcdata = np.loadtxt('data/LC3.data', delimiter='\t')

lctimes = lcdata[:,0]
fluxes = lcdata[:,1]
fsigmas = lcdata[:,2]

#################PLOT LC

plt.plot(lctimes,fluxes)
plt.show()

#################GENERATE DATA

b = phoebe.default_binary()
b.add_constraint('semidetached', 'primary')
b['requiv@constraint@primary']
b['requiv_max@constraint@primary']

b.set_value('latex_repr', component='binary', value='orb')
b.set_value('latex_repr', component='primary', value='1')
b.set_value('latex_repr', component='secondary', value='2')

#################SET PARAMETERS
#THESE ARE THE VALUES I GOT AS SOLUTION USING THE WEB INTERFACE.

b.set_value(qualifier='period', context='component', component='binary', value='6.45895 d')
b.set_value(qualifier='teff', context='component', component='primary', value='15700.0 K')
b.set_value(qualifier='teff', context='component', component='secondary', value='5000.0 K')
b.set_value(qualifier='ld_mode_bol', context='component', component='primary', value='manual')
b.set_value(qualifier='gravb_bol', context='component', component='primary', value='0.9')
b.set_value(qualifier='irrad_frac_refl_bol', context='component', component='primary', value='1')
b.set_value(qualifier='sma', context='component', component='binary', value='26.0 solRad')
b.set_value(qualifier='requiv', context='component', component='secondary', value='1.0 solRad')
b.set_value(qualifier='syncpar', context='component', component='primary', value='14.00')

#################ADD LC DATASETS

b.add_dataset('lc',
        compute_phases=phoebe.linspace(0,1,101),
                times=lctimes, 
                fluxes=fluxes, 
                sigmas=fsigmas, 
                dataset='lc01',
                overwrite=True)
b.set_value_all('ld_mode', 'manual')
b.set_value_all('atm', 'blackbody')
b.set_value(qualifier='ld_mode_bol', context='component', component='primary', value='manual')

#################COMPUTE MODEL

b.run_compute(irrad_method='none')
b.plot(x='phase',show=True)

### UP TO HERE IT WORKS FINE and it represent the same curve I was getting in the web interface.

#################ADOPT A SOLUTION METHOD LC ESTIMATOR

b.add_solver('estimator.lc_geometry', lc_datasets='lc01', solver='lcgeom')
print(b.filter(solver='lcgeom'))
b.run_solver('lcgeom', solution='lcgeom_solution', max_computations=15000, detach=True, overwrite=True)
print(b.filter(solution='lcgeom_solution'))
plt.plot(b.filter(solution='lcgeom_solution'), color='r', label='solution')
plt.plot(lctimes, fluxes, color='g', label= 'data')
plt.show()
b.plot(x='phases', ls='-', legend=True, show=True)

the problem is that I see that the system is not adapting any model, but giving me straight lines all the time... I do get the following message though: so I am not sure if maybe I do need to fix something with my phoebe installation or is still something off with my script.

thanks again and apologies for all the troubles

[mcregogarcia]

for some clarity, the bundle is the bundle used in the web server, and the first image is the result from the lc estimator in the web version, vs the rest of the pictures what the script gives me

[kecnry]

It doesn't look like you're loading a bundle at any point... I'm assuming OPCION1.py was exported from the web UI, but that is just meant as a log of commands and I wouldn't suggest running that directly (I'm guessing that file is where the error message is coming from, not from the script you attached in your message?).

the problem is that I see that the system is not adapting any model

You need to call run_compute manually whenever you want a new synthetic forward-model to be generated. Plotting only plots the models that have already been computed.

but giving me straight lines all the time

I suspect they're not actually straight lines (try b.plot(context='model', show=True) to only show the model)... but rather because the range on your observation fluxes is much larger than the model (which by default will be roughly 0-2). Phoebe will never produce negative fluxes, so you'll also probably need to convert the units on your input data or renormalize before importing them.

If you're just trying to modify that or a single example, I would suggest going back and going through the main tutorials to get familiar with how the python package of phoebe works as it does have a bit of a learning curve before people generally feel comfortable using it to fit a new system.

[mcregogarcia]

Thanks. I have fixed a couple of problems

#################SET-UP

import phoebe
from phoebe import u #units
import numpy as np
import matplotlib.pyplot as plt
from reportlab.lib.pagesizes import letter
from reportlab.pdfgen import canvas

logger = phoebe.logger()

#################LOAD LC

lcdata = np.loadtxt('data/LC9.data', delimiter='\t')

lctimes = lcdata[:,0]
fluxes = lcdata[:,1]
fsigmas = lcdata[:,2]

#################PLOT LC

plt.plot(lctimes,fluxes)
plt.show()

#################GENERATE DATA

b = phoebe.default_binary()
b.add_constraint('semidetached', 'primary')
b['requiv@constraint@primary']
b['requiv_max@constraint@primary']

b.set_value('latex_repr', component='binary', value='orb')
b.set_value('latex_repr', component='primary', value='1')
b.set_value('latex_repr', component='secondary', value='2')

#################SET PARAMETERS
#THESE ARE THE VALUES I GOT AS SOLUTION USING THE WEB INTERFACE.

b.set_value(qualifier='period', context='component', component='binary', value='6.45895 d')
b.set_value(qualifier='teff', context='component', component='primary', value='15700.0 K')
b.set_value(qualifier='teff', context='component', component='secondary', value='5000.0 K')
b.set_value(qualifier='ld_mode_bol', context='component', component='primary', value='manual')
b.set_value(qualifier='gravb_bol', context='component', component='primary', value='0.9')
b.set_value(qualifier='irrad_frac_refl_bol', context='component', component='primary', value='1')
b.set_value(qualifier='sma', context='component', component='binary', value='26.0 solRad')
b.set_value(qualifier='requiv', context='component', component='secondary', value='1.0 solRad')
b.set_value(qualifier='syncpar', context='component', component='primary', value='14.00')

#################ADD LC DATASETS

b.add_dataset('lc',
        compute_phases=phoebe.linspace(0,1,101),
                times=lctimes, 
                fluxes=fluxes, 
                sigmas=fsigmas, 
                dataset='lc02',
                overwrite=True)
b.set_value_all('ld_mode', 'manual')
b.set_value_all('atm', 'blackbody')
b.set_value(qualifier='ld_mode_bol', context='component', component='primary', value='manual')

#################COMPUTE MODEL

b.run_compute(irrad_method='none')
b.plot(x='times', context='model', show=True)

#################ADOPT A SOLUTION METHOD LC ESTIMATOR

b.add_solver('estimator.lc_geometry', lc_datasets='lc01', solver='lcgeom')
print(b.filter(solver='lcgeom'))
b.run_solver('lcgeom', solution='lcgeom_solution', max_computations=15000, detach=True, overwrite=True)
print(b.filter(solution='lcgeom_solution'))
plt.plot(b.filter(solution='lcgeom_solution'), color='r', label='solution')
plt.plot(lctimes, fluxes, color='g', label= 'data')
plt.show()
b.plot(x='times', ls='-', legend=True, show=True)

So now> I have loaded my light curve the system creates a model based on my defined parameters. the system applies the lc solver to the model that phoebe creates with the parameters, instead of the loaded light curve.....

are there any way to fix this? I have reviewed the tutorials, but this I cannot see how to fix.

Thanks in advance. Miguel

[kecnry]

Calling run_solver only proposes new values to change. You should then manually inspect them and decide whether to adopt some or all of the values using adopt_solver, and then will need to call run_compute again to create a new forward model based on the new values of your system parameters (see the solvers tutorial for details)

[kecnry]

Feel free to re-open or create a new issue if you're still running into any issues.