Contact systems in LC Estimators return errors

[jjuvan]

Hello,

I am currently attempting extract various parameters of observed eclipsing systems by following the LC Estimators documentation. I have been able to run all three models with an EA system, but am having issues with both EB and EW systems. For each star, their phases, magnitudes (which are converted into fluxes within the program), and errors are given to the program (data can be found here: EA System, EB System, EW System). The code can be found below, while images to each system can be found here: EA System, EB System, EW System

def run_phoebe(data, if_contact):
    
    import phoebe
    from phoebe import u
    import numpy as np
    
    if if_contact == 'detached':
        b = phoebe.default_binary()
    elif if_contact == 'contact':
        b = phoebe.default_contact_binary()
        
    time = []
    flux = []
    sigma = []
    for obs in data:
        time.append(obs[0])
        obs_flux = 3060 * (10**(-26)) * (10 ** (-obs[1]/2.5))
        flux.append(obs_flux)
        upper = 3060 * (10**(-26)) * (10 ** (-(obs[1]-obs[2])/2.5))
        lower = 3060 * (10**(-26)) * (10 ** (-(obs[1]+obs[2])/2.5))
        sigma.append(upper-lower)
            
    b.add_dataset('lc', times=time, fluxes=flux, sigmas=sigma)
    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(irrad_method='none')
    _ = b.plot(x='phase', ls='-', legend=True, show=True)
        
    def flip_constraints(model_type):
        if model_type == 'knn' or model_type == 'mlp':
            try:
                b.flip_constraint('esinw', solve_for='ecc')
                b.flip_constraint('ecosw', solve_for='per0')
            except:
                placeholder = False
        if model_type == 'lc_geom':
            try: 
                b.flip_constraint('ecc', solve_for='esinw')
                b.flip_constraint('per0', solve_for='ecosw')
            except:
                placeholder = False
    
        #default
        #b.run_compute(model='default')
        #print('Default ran successfully \n\n')

    if if_contact == 'detached':
        b.flip_constraint('requivsumfrac', solve_for='requiv@primary')
    elif if_contact == 'contact':
        b.flip_constraint('requivsumfrac', solve_for='sma')
    b.flip_constraint('teffratio', solve_for='teff@secondary')

    #mlp
    try:
        model = 'mlp'
        b.add_solver('estimator.ebai', ebai_method='mlp', solver='ebai_mlp')
        b['phase_bin@ebai_mlp'] = False
        b.run_solver('ebai_mlp', solution='ebai_mlp_solution')
        flip_constraints(model)
        b.adopt_solution('ebai_mlp_solution')
        b.run_compute(model='ebai_mlp_model')
        print()
        print('MLP ran successfully \n\n')
        print(b.adopt_solution('ebai_mlp_solution'))
        print('\n\n')
    
    except Exception as e:
        print(e)
        print('*** ERROR: MLP model could not run on current data, trying KNN *** \n\n')
    
    #knn
    try:
        model = 'knn'
        b.add_solver('estimator.ebai', ebai_method='knn', solver='ebai_knn')
        b['phase_bin@ebai_knn'] = False
        b.run_solver(solver='ebai_knn', solution='ebai_knn_solution')
        flip_constraints(model)
        b.adopt_solution('ebai_knn_solution')
        b.run_compute(model='ebai_knn_model')
        print()
        print('KNN ran successfully \n\n')
        print(b.adopt_solution('ebai_knn_solution'))
        print('\n\n')
    except Exception as e:
        print(e)
        print('*** ERROR: KNN model could not run on current data, trying LC_GEOM *** \n\n')
    
    
    #lc_geom
    try:
        model = 'lc_geom'
        b.add_solver('estimator.lc_geometry', solver='lcgeom')
        b['phase_bin@lcgeom'] = False
        b.run_solver('lcgeom', solution='lcgeom_solution')
        flip_constraints(model)
        b.adopt_solution('lcgeom_solution')
        b.run_compute(model='lcgeom_model')
        print()
        print('LC_GEOM ran successfully')
        print(b.adopt_solution('lcgeom_solution'))
        print('\n\n')
    except Exception as e:
        print(e)
        print('*** ERROR: LC_GEOM model could not run on current data, plotting data *** \n\n')
    
    _ = b.plot(x='phase', ls='-', legend=True, show=True)

import numpy as np
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("indata")
parser.add_argument("contact")

args = parser.parse_args()
indata = args.indata
contact_status = args.contact

use = (np.loadtxt(indata)).tolist()
plot = run_phoebe(use,contact_status)

The EA system can run each model (MLP, KNN, and LC_GEOM) successfully using the default_binary bundle (see image here), however when using the default_contact_binary bundle (and only the KNN model) for the EW system (since using the default_binary bundle causes the periastron to overflow), I receive the following error:

If instead I attempt to use the EB system, I get this error instead:

I have attempted to use the default_contact_binary bundle on simulated data, however the results for the KNN model did not agree with the data. I have the code for the simulated contact system and its plot below.

import phoebe
from phoebe import u
import numpy as np
b = phoebe.default_contact_binary()
#for contact: ecc@binary=0, 
b.set_value('ecc', 0)
b.set_value('per0', 0.5)
b.set_value('teff@primary', 7000)
b.set_value('teff@secondary', 6000)
b.set_value('sma@binary', 7)
b.set_value('incl@binary', 80)
b.set_value('q', 0.3)
b.set_value('t0_supconj', 0.1)
b.set_value('requiv@primary', 3.5)

# add an lc dataset
b.add_dataset('lc', compute_phases=phoebe.linspace(0,1,201))
b.set_value_all('ld_mode', 'manual')
b.set_value_all('ld_mode_bol', 'manual')
b.set_value_all('atm', 'blackbody')

#compute the model
b.run_compute(irrad_method='none')


# extract the arrays from the model that we'll use as observables in the next step
times = b.get_value('times', context='model', dataset='lc01')
# here we're adding noise to the fluxes as well to make the fake data more "realistic"
np.random.seed(0) # to ensure reproducibility with added noise
fluxes = b.get_value('fluxes', context='model', dataset='lc01') + np.random.normal(size=times.shape) * 0.005
sigmas_lc = np.ones_like(times) * 0.04

b = phoebe.default_binary()
b.add_dataset('lc', times=times, fluxes=fluxes, sigmas=sigmas_lc)
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='phase', show=True)

b.flip_constraint('requivsumfrac',solve_for='sma')
b.flip_constraint('teffratio',solve_for='teff@primary')
b.flip_constraint('esinw', solve_for='ecc')
b.flip_constraint('ecosw', solve_for='per0')

model = 'knn'
b.add_solver('estimator.ebai', ebai_method='knn', solver='ebai_knn')
b['phase_bin@ebai_knn'] = False
b.run_solver(solver='ebai_knn', solution='ebai_knn_solution')
#print(b)
b.adopt_solution('ebai_knn_solution')
b.run_compute(model='ebai_knn_model')
print()
print('KNN ran successfully \n\n')
print(b.adopt_solution('ebai_knn_solution'))
print('\n\n')
_ = b.plot(x='phase', ls='-', legend=True, show=True)

Is there any way to correct for either error with the observed EW and EB data? Thank you in advance and sorry for the inconvenience.

[aprsa]

LC estimators are principally designed for detached binaries. They may work on semi-detached but are unlikely to do a good job for contact binaries. Existing LC estimators can certainly be generalized, or a dedicated contact-binary estimator could be implemented, but that's not on our immediate radar. PRs are always welcome, though! :)