Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

lc_geometry propose compute_phases to cover eclipses more heavily #898

Draft
wants to merge 2 commits into
base: release-2.5
Choose a base branch
from
Draft

lc_geometry propose compute_phases to cover eclipses more heavily #898

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4ee8b4d
fix bug in adopting mask_phases from lc_geometry
kecnry Jun 11, 2024
846020f
WIP: propose compute_phases in lc_geometry
kecnry Jun 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions phoebe/frontend/bundle.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13731,7 +13731,7 @@ def _get_adopt_inds_uniqueids(self, solution_ps, **kwargs):

b_uniqueids = self.uniqueids

adoptable_ps = self.get_adjustable_parameters(exclude_constrained=False) + self.filter(qualifier='mask_phases', context='dataset', **_skip_filter_checks)
adoptable_ps = self.get_adjustable_parameters(exclude_constrained=False) + self.filter(qualifier='mask_phases', context='dataset', **_skip_filter_checks) + self.filter(qualifier='compute_phases', context='dataset', **_skip_filter_checks)
if np.all([uniqueid.split('[')[0] in b_uniqueids for uniqueid in fitted_uniqueids]):
fitted_ps = adoptable_ps.filter(uniqueid=[uniqueid.split('[')[0] for uniqueid in fitted_uniqueids], **_skip_filter_checks)
else:
Expand Down Expand Up @@ -13922,8 +13922,9 @@ def adopt_solution(self, solution=None,
fitted_values = fitted_values * solution_ps.get_value(qualifier='period_factor', period_factor=kwargs.get('period_factor', None), **_skip_filter_checks)

if solver_kind in ['lc_geometry']:

print("*** adopt_uniqueids", adopt_uniqueids)
adopt_qualifiers = [self.get_parameter(uniqueid=uniqueid, **_skip_filter_checks).qualifier for uniqueid in adopt_uniqueids]
print("*** adopt_qualifiers", adopt_qualifiers)
if 'mask_phases' in adopt_qualifiers and 't0_supconj' in adopt_qualifiers:
# then we need to shift mask phases by the phase-shift introduced by the change in t0_supconj
logger.info("shifting mask_phases by phase-shift caused by change in t0_supconj")
Expand All @@ -13935,11 +13936,11 @@ def adopt_solution(self, solution=None,
t0_supconj_ind = adopt_qualifiers.index('t0_supconj')

t0_supconj_old = self.get_value(uniqueid=adopt_uniqueids[t0_supconj_ind], unit=u.d, **_skip_filter_checks)
t0_supconj_new = fitted_values[t0_supconj_ind]
t0_supconj_new = fitted_values[adopt_inds[t0_supconj_ind]]

phase_shift = self.to_phase(t0_supconj_new) - self.to_phase(t0_supconj_old)

fitted_values[mask_phases_ind] = [ph-phase_shift for ph in [ecl_ph for ecl_ph in fitted_values[mask_phases_ind]]]
fitted_values[adopt_inds[mask_phases_ind]] = [ph-phase_shift for ph in [ecl_ph for ecl_ph in fitted_values[adopt_inds[mask_phases_ind]]]]

for uniqueid, value, unit in zip(adopt_uniqueids, fitted_values[adopt_inds], fitted_units[adopt_inds]):
uniqueid, index = _extract_index_from_string(uniqueid)
Expand Down
71 changes: 65 additions & 6 deletions phoebe/solverbackends/solverbackends.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -607,6 +607,57 @@ def _get_packet_and_solution(self, b, solver, **kwargs):

return kwargs, _parameters.ParameterSet(solution_params)

@staticmethod
def propose_compute_phases(phases, fluxes, n_points, n_deriv):
from scipy.interpolate import interp1d

# derivative and the interpolation between points
dydx = np.gradient(fluxes,phases)
dydx_int = interp1d(phases, dydx, kind='linear', bounds_error=False, fill_value='extrapolate', assume_sorted=False)

#count the total length of the phased light curve and the unit dL
L = 0.0
for i in range(len(phases)-1):
L += np.sqrt( ( phases[i+1]-phases[i] )**2 + ( fluxes[i+1]-fluxes[i] )**2 )

dL = L/n_points
dL

#create a mesh of target bin boundaries with iterative correction of bin position
bin_edges = [0.0]

while True:
dx = dL/np.sqrt(1.0+dydx_int(bin_edges[-1])**2)
for j in range(10):
dx = dL/np.sqrt(1.0+dydx_int(bin_edges[-1]+dx/2.0)**2)
new_edge = bin_edges[-1]+dx
if new_edge < 1.0:
bin_edges.append(new_edge)
else:
break

bin_edges.append(1.0)

# final bin with the evaluation of total number of bins
pbin = []
fbin = []

for i in range(len(bin_edges)-1):
low = bin_edges[i]
hi = bin_edges[i+1]
try:
y = np.median(fluxes[((phases>low)&(phases<=hi))])
x = np.median(phases[((phases>low)&(phases<=hi))])
pbin.append(x)
fbin.append(y)
except:
continue

# len(pbin)

return bin_edges


def run_worker(self, b, solver, compute=None, **kwargs):
if mpi.within_mpirun:
raise NotImplementedError("mpi support for lc_geometry not yet implemented")
Expand Down Expand Up @@ -657,9 +708,15 @@ def run_worker(self, b, solver, compute=None, **kwargs):
t0_supconj_new = eb_params._t0_from_geometry(times, period=period, t0_supconj = t0_supconj_old, t0_near_times = t0_near_times)
#compute_eclipse_params(phases, fluxes, sigmas, fit_result=fit_result, diagnose=diagnose)

if analytical_model == 'two-gaussian':
analytic_fluxes = model.models
else:
analytic_fluxes = {'polyfit': model.model}

edges = eclipse_dict.get('eclipse_edges')
mask_phases = [(edges[0]-eclipse_dict.get('primary_width')*0.3, edges[1]+eclipse_dict.get('primary_width')*0.3), (edges[2]-eclipse_dict.get('secondary_width')*0.3, edges[3]+eclipse_dict.get('secondary_width')*0.3)]
fitted_params += b.filter(qualifier='mask_phases', dataset=lc_datasets, **_skip_filter_checks).to_list()
fitted_params += b.filter(qualifier='mask_phases', context='dataset', dataset=lc_datasets, **_skip_filter_checks).to_list()
fitted_params += b.filter(qualifier='compute_phases', context='dataset', dataset=lc_datasets, **_skip_filter_checks).to_list()

fitted_uniqueids = [p.uniqueid for p in fitted_params]
fitted_twigs = [p.twig for p in fitted_params]
Expand All @@ -668,6 +725,10 @@ def run_worker(self, b, solver, compute=None, **kwargs):
if fit_eclipses:
fitted_values += [eb_params.rratio, eb_params.incl]
fitted_values += [mask_phases for ds in lc_datasets]
proposed_compute_phases = self.propose_compute_phases(model.phases, analytic_fluxes[model.best_fit['func']],
n_points=100,
n_deriv=10/eclipse_dict.get('primary_width')/2.)
fitted_values += [proposed_compute_phases for ds in lc_datasets]

fitted_units = [u.d.to_string(),
u.dimensionless_unscaled.to_string(),
Expand All @@ -676,7 +737,8 @@ def run_worker(self, b, solver, compute=None, **kwargs):
u.dimensionless_unscaled.to_string()]
if fit_eclipses:
fitted_units += [u.dimensionless_unscaled.to_string(), u.deg.to_string()]
fitted_units += [u.dimensionless_unscaled.to_string() for ds in lc_datasets]
fitted_units += [u.dimensionless_unscaled.to_string() for ds in lc_datasets] # mask_phases
fitted_units += [u.dimensionless_unscaled.to_string() for ds in lc_datasets] # compute_phases

return_ = [{'qualifier': 'primary_width', 'value': eclipse_dict.get('primary_width')},
{'qualifier': 'secondary_width', 'value': eclipse_dict.get('secondary_width')},
Expand All @@ -697,10 +759,7 @@ def run_worker(self, b, solver, compute=None, **kwargs):
]

if kwargs.get('expose_model', True):
if analytical_model == 'two-gaussian':
analytic_fluxes = model.models
else:
analytic_fluxes = {'polyfit': model.model}

return_ += [{'qualifier': 'analytic_phases', 'value': model.phases},
{'qualifier': 'analytic_fluxes', 'value': analytic_fluxes},
{'qualifier': 'analytic_best_model', 'value': model.best_fit['func']}
Expand Down
Loading