From 09e83f5f2065620293bcff63b61c4ebea6f609f2 Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Tue, 11 Jun 2024 05:52:56 -0400 Subject: [PATCH 1/3] refactor features and add hook for dataset features --- phoebe/__init__.py | 1 + phoebe/backend/universe.py | 305 +------------------------ phoebe/features/__init__.py | 2 + phoebe/features/component_features.py | 309 ++++++++++++++++++++++++++ phoebe/features/dataset_features.py | 58 +++++ phoebe/features/gaussian_processes.py | 164 ++++++++++++++ phoebe/frontend/bundle.py | 206 +++-------------- phoebe/parameters/constraint.py | 8 +- phoebe/parameters/feature.py | 21 ++ 9 files changed, 596 insertions(+), 478 deletions(-) create mode 100644 phoebe/features/__init__.py create mode 100644 phoebe/features/component_features.py create mode 100644 phoebe/features/dataset_features.py create mode 100644 phoebe/features/gaussian_processes.py diff --git a/phoebe/__init__.py b/phoebe/__init__.py index aac418216..e1e766d0f 100644 --- a/phoebe/__init__.py +++ b/phoebe/__init__.py @@ -384,6 +384,7 @@ def progressbars(self): from . import dynamics as dynamics from . import distortions as distortions from . import algorithms as algorithms +from . import features as features import libphoebe # Shortcut to building logger diff --git a/phoebe/backend/universe.py b/phoebe/backend/universe.py index 4a25ac711..c97035dd3 100644 --- a/phoebe/backend/universe.py +++ b/phoebe/backend/universe.py @@ -9,6 +9,7 @@ from phoebe.distortions import roche, rotstar from phoebe.backend import eclipse, oc_geometry, mesh, mesh_wd from phoebe.utils import _bytes +from phoebe.features import component_features import libphoebe from phoebe import u @@ -1249,7 +1250,7 @@ def from_bundle(cls, b, component, compute=None, feature_ps = b.get_feature(feature=feature, **_skip_filter_checks) if feature_ps.component != component: continue - feature_cls = globals()[feature_ps.kind.title()] + feature_cls = getattr(component_features, feature_ps.kind.title()) features.append(feature_cls.from_bundle(b, feature)) if conf.devel: @@ -3027,304 +3028,4 @@ def populate_observable(self, time, kind, dataset, **kwargs): """ for half in self._halves: - half.populate_observable(time, kind, dataset, **kwargs) - - -################################################################################ -################################################################################ -################################################################################ - - -class Feature(object): - """ - Note that for all features, each of the methods below will be called. So - changing the coordinates WILL affect the original/intrinsic loggs which - will then be used as input for that method call. - - In other words, its probably safest if each feature only overrides a - SINGLE one of the methods. Overriding multiple methods should be done - with great care. - - Each feature may or may not require recomputing a mesh, depending on the - kind of change it exacts to the mesh. For example, pulsations will require - recomputing a mesh while spots will not. By default, the mesh will be - recomputed (set in this superclass' `__init__()` method) but inherited - classes should overload `self._remeshing_required`. - """ - def __init__(self, *args, **kwargs): - pass - - @property - def _remeshing_required(self): - return True - - @property - def proto_coords(self): - """ - Override this to True if all methods (except process_coords*... those - ALWAYS expect protomesh coordinates) are expecting coordinates - in the protomesh (star) frame-of-reference rather than the - current in-orbit system frame-of-reference. - """ - return False - - def process_coords_for_computations(self, coords_for_computations, s, t): - """ - Method for a feature to process the coordinates. Coordinates are - processed AFTER scaling but BEFORE being placed in orbit. - - NOTE: coords_for_computations affect physical properties only and - not geometric properties (areas, eclipse detection, etc). If you - want to override geometric properties, use the hook for - process_coords_for_observations as well. - - Features that affect coordinates_for_computations should override - this method - """ - return coords_for_computations - - def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t): - """ - Method for a feature to process the coordinates. Coordinates are - processed AFTER scaling but BEFORE being placed in orbit. - - NOTE: coords_for_observations affect the geometry only (areas of each - element and eclipse detection) but WILL NOT affect any physical - parameters (loggs, teffs, intensities). If you want to override - physical parameters, use the hook for process_coords_for_computations - as well. - - Features that affect coordinates_for_observations should override this method. - """ - return coords_for_observations - - def process_loggs(self, loggs, coords, s=np.array([0., 0., 1.]), t=None): - """ - Method for a feature to process the loggs. - - Features that affect loggs should override this method - """ - return loggs - - def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): - """ - Method for a feature to process the teffs. - - Features that affect teffs should override this method - """ - return teffs - -class Spot(Feature): - def __init__(self, colat, longitude, dlongdt, radius, relteff, t0, **kwargs): - """ - Initialize a Spot feature - """ - super(Spot, self).__init__(**kwargs) - self._colat = colat - self._longitude = longitude - self._radius = radius - self._relteff = relteff - self._dlongdt = dlongdt - self._t0 = t0 - - @classmethod - def from_bundle(cls, b, feature): - """ - Initialize a Spot feature from the bundle. - """ - - feature_ps = b.get_feature(feature=feature, **_skip_filter_checks) - - colat = feature_ps.get_value(qualifier='colat', unit=u.rad, **_skip_filter_checks) - longitude = feature_ps.get_value(qualifier='long', unit=u.rad, **_skip_filter_checks) - - if len(b.hierarchy.get_stars())>=2: - star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks) - orbit_ps = b.get_component(component=b.hierarchy.get_parent_of(feature_ps.component), **_skip_filter_checks) - # TODO: how should this handle dpdt? - - # we won't use syncpar directly because that is defined wrt sidereal period and we want to make sure - # this translated to roche longitude correctly. In the non-apsidal motion case - # syncpar = period_anom_orb / period_star - period_anom_orb = orbit_ps.get_value(qualifier='period_anom', unit=u.d, **_skip_filter_checks) - period_star = star_ps.get_value(qualifier='period', unit=u.d, **_skip_filter_checks) - dlongdt = 2*pi * (period_anom_orb/period_star - 1) / period_anom_orb - else: - star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks) - dlongdt = star_ps.get_value(qualifier='freq', unit=u.rad/u.d, **_skip_filter_checks) - longitude += np.pi/2 - - radius = feature_ps.get_value(qualifier='radius', unit=u.rad, **_skip_filter_checks) - relteff = feature_ps.get_value(qualifier='relteff', unit=u.dimensionless_unscaled, **_skip_filter_checks) - - t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks) - - return cls(colat, longitude, dlongdt, radius, relteff, t0) - - @property - def _remeshing_required(self): - return False - - @property - def proto_coords(self): - """ - """ - return True - - def pointing_vector(self, s, time): - """ - s is the spin vector in roche coordinates - time is the current time - """ - t = time - self._t0 - longitude = self._longitude + self._dlongdt * t - - # define the basis vectors in the spin (primed) coordinates in terms of - # the Roche coordinates. - # ez' = s - # ex' = (ex - s(s.ex)) /|i - s(s.ex)| - # ey' = s x ex' - ex = np.array([1., 0., 0.]) - ezp = s - exp = (ex - s*np.dot(s,ex)) - eyp = np.cross(s, exp) - - return np.sin(self._colat)*np.cos(longitude)*exp +\ - np.sin(self._colat)*np.sin(longitude)*eyp +\ - np.cos(self._colat)*ezp - - def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): - """ - Change the local effective temperatures for any values within the - "cone" defined by the spot. Any teff within the spot will have its - current value multiplied by the "relteff" factor - - :parameter array teffs: array of teffs for computations - :parameter array coords: array of coords for computations - :t float: current time - """ - if t is None: - # then assume at t0 - t = self._t0 - - pointing_vector = self.pointing_vector(s,t) - logger.debug("spot.process_teffs at t={} with pointing_vector={} and radius={}".format(t, pointing_vector, self._radius)) - - cos_alpha_coords = np.dot(coords, pointing_vector) / np.linalg.norm(coords, axis=1) - cos_alpha_spot = np.cos(self._radius) - - filter_ = cos_alpha_coords > cos_alpha_spot - teffs[filter_] = teffs[filter_] * self._relteff - - return teffs - -class Pulsation(Feature): - def __init__(self, radamp, freq, l=0, m=0, tanamp=0.0, teffext=False, **kwargs): - self._freq = freq - self._radamp = radamp - self._l = l - self._m = m - self._tanamp = tanamp - - self._teffext = teffext - - @classmethod - def from_bundle(cls, b, feature): - """ - Initialize a Pulsation feature from the bundle. - """ - - feature_ps = b.get_feature(feature=feature, **_skip_filter_checks) - freq = feature_ps.get_value(qualifier='freq', unit=u.d**-1, **_skip_filter_checks) - radamp = feature_ps.get_value(qualifier='radamp', unit=u.dimensionless_unscaled, **_skip_filter_checks) - l = feature_ps.get_value(qualifier='l', unit=u.dimensionless_unscaled, **_skip_filter_checks) - m = feature_ps.get_value(qualifier='m', unit=u.dimensionless_unscaled, **_skip_filter_checks) - teffext = feature_ps.get_value(qualifier='teffext', **_skip_filter_checks) - - GM = c.G.to('solRad3 / (solMass d2)').value*b.get_value(qualifier='mass', component=feature_ps.component, context='component', unit=u.solMass, **_skip_filter_checks) - R = b.get_value(qualifier='rpole', component=feature_ps.component, section='component', unit=u.solRad, **_skip_filter_checks) - - tanamp = GM/R**3/freq**2 - - return cls(radamp, freq, l, m, tanamp, teffext) - - @property - def proto_coords(self): - """ - """ - return True - - def dYdtheta(self, m, l, theta, phi): - if abs(m) > l: - return 0 - - # TODO: just a quick hack - if abs(m+1) > l: - last_term = 0.0 - else: - last_term = Y(m+1, l, theta, phi) - - return m/np.tan(theta)*Y(m, l, theta, phi) + np.sqrt((l-m)*(l+m+1))*np.exp(-1j*phi)*last_term - - def dYdphi(self, m, l, theta, phi): - return 1j*m*Y(m, l, theta, phi) - - def process_coords_for_computations(self, coords_for_computations, s, t): - """ - """ - if self._teffext: - return coords_for_computations - - x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1)) - theta = np.arccos(z/r) - phi = np.arctan2(y, x) - - xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) - xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) - xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) - - new_coords = np.zeros(coords_for_computations.shape) - new_coords[:,0] = coords_for_computations[:,0] + xi_r * np.sin(theta) * np.cos(phi) - new_coords[:,1] = coords_for_computations[:,1] + xi_r * np.sin(theta) * np.sin(phi) - new_coords[:,2] = coords_for_computations[:,2] + xi_r * np.cos(theta) - - return new_coords - - def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t): - """ - Displacement equations: - - xi_r(r, theta, phi) = a(r) Y_lm (theta, phi) exp(-i*2*pi*f*t) - xi_theta(r, theta, phi) = b(r) dY_lm/dtheta (theta, phi) exp(-i*2*pi*f*t) - xi_phi(r, theta, phi) = b(r)/sin(theta) dY_lm/dphi (theta, phi) exp(-i*2*pi*f*t) - - where: - - b(r) = a(r) GM/(R^3*f^2) - """ - # TODO: we do want to displace the coords_for_observations, but the x,y,z,r below are from the ALSO displaced coords_for_computations - # if not self._teffext: - # return coords_for_observations - - x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1)) - theta = np.arccos(z/r) - phi = np.arctan2(y, x) - - xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) - xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) - xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) - - new_coords = np.zeros(coords_for_observations.shape) - new_coords[:,0] = coords_for_observations[:,0] + xi_r * np.sin(theta) * np.cos(phi) - new_coords[:,1] = coords_for_observations[:,1] + xi_r * np.sin(theta) * np.sin(phi) - new_coords[:,2] = coords_for_observations[:,2] + xi_r * np.cos(theta) - - return new_coords - - def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): - """ - """ - if not self._teffext: - return teffs - - raise NotImplementedError("teffext=True not yet supported for pulsations") + half.populate_observable(time, kind, dataset, **kwargs) \ No newline at end of file diff --git a/phoebe/features/__init__.py b/phoebe/features/__init__.py new file mode 100644 index 000000000..df7508021 --- /dev/null +++ b/phoebe/features/__init__.py @@ -0,0 +1,2 @@ +from .dataset_features import * +from .component_features import * \ No newline at end of file diff --git a/phoebe/features/component_features.py b/phoebe/features/component_features.py new file mode 100644 index 000000000..7b321ea2b --- /dev/null +++ b/phoebe/features/component_features.py @@ -0,0 +1,309 @@ +import numpy as np +import astropy.units as u + +import logging +logger = logging.getLogger("COMPONENT_FEATURES") +logger.addHandler(logging.NullHandler()) + + +__all__ = ['ComponentFeature', 'Spot', 'Pulsation'] + +_skip_filter_checks = {'check_default': False, 'check_visible': False} + +class ComponentFeature(object): + """ + Note that for all features, each of the methods below will be called. So + changing the coordinates WILL affect the original/intrinsic loggs which + will then be used as input for that method call. + + In other words, its probably safest if each feature only overrides a + SINGLE one of the methods. Overriding multiple methods should be done + with great care. + + Each feature may or may not require recomputing a mesh, depending on the + kind of change it exacts to the mesh. For example, pulsations will require + recomputing a mesh while spots will not. By default, the mesh will be + recomputed (set in this superclass' `__init__()` method) but inherited + classes should overload `self._remeshing_required`. + """ + def __init__(self, *args, **kwargs): + pass + + @classmethod + def from_bundle(cls, b, feature): + raise NotImplementedError("from_bundle must be implemented in the feature subclass") + + @property + def _remeshing_required(self): + return True + + @property + def proto_coords(self): + """ + Override this to True if all methods (except process_coords*... those + ALWAYS expect protomesh coordinates) are expecting coordinates + in the protomesh (star) frame-of-reference rather than the + current in-orbit system frame-of-reference. + """ + return False + + def process_coords_for_computations(self, coords_for_computations, s, t): + """ + Method for a feature to process the coordinates. Coordinates are + processed AFTER scaling but BEFORE being placed in orbit. + + NOTE: coords_for_computations affect physical properties only and + not geometric properties (areas, eclipse detection, etc). If you + want to override geometric properties, use the hook for + process_coords_for_observations as well. + + Features that affect coordinates_for_computations should override + this method + """ + return coords_for_computations + + def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t): + """ + Method for a feature to process the coordinates. Coordinates are + processed AFTER scaling but BEFORE being placed in orbit. + + NOTE: coords_for_observations affect the geometry only (areas of each + element and eclipse detection) but WILL NOT affect any physical + parameters (loggs, teffs, intensities). If you want to override + physical parameters, use the hook for process_coords_for_computations + as well. + + Features that affect coordinates_for_observations should override this method. + """ + return coords_for_observations + + def process_loggs(self, loggs, coords, s=np.array([0., 0., 1.]), t=None): + """ + Method for a feature to process the loggs. + + Features that affect loggs should override this method + """ + return loggs + + def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): + """ + Method for a feature to process the teffs. + + Features that affect teffs should override this method + """ + return teffs + +class Spot(ComponentFeature): + def __init__(self, colat, longitude, dlongdt, radius, relteff, t0, **kwargs): + """ + Initialize a Spot feature + """ + super(Spot, self).__init__(**kwargs) + self._colat = colat + self._longitude = longitude + self._radius = radius + self._relteff = relteff + self._dlongdt = dlongdt + self._t0 = t0 + + @classmethod + def from_bundle(cls, b, feature): + """ + Initialize a Spot feature from the bundle. + """ + + feature_ps = b.get_feature(feature=feature, **_skip_filter_checks) + + colat = feature_ps.get_value(qualifier='colat', unit=u.rad, **_skip_filter_checks) + longitude = feature_ps.get_value(qualifier='long', unit=u.rad, **_skip_filter_checks) + + if len(b.hierarchy.get_stars())>=2: + star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks) + orbit_ps = b.get_component(component=b.hierarchy.get_parent_of(feature_ps.component), **_skip_filter_checks) + # TODO: how should this handle dpdt? + + # we won't use syncpar directly because that is defined wrt sidereal period and we want to make sure + # this translated to roche longitude correctly. In the non-apsidal motion case + # syncpar = period_anom_orb / period_star + period_anom_orb = orbit_ps.get_value(qualifier='period_anom', unit=u.d, **_skip_filter_checks) + period_star = star_ps.get_value(qualifier='period', unit=u.d, **_skip_filter_checks) + dlongdt = 2 * np.pi * (period_anom_orb/period_star - 1) / period_anom_orb + else: + star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks) + dlongdt = star_ps.get_value(qualifier='freq', unit=u.rad/u.d, **_skip_filter_checks) + longitude += np.pi/2 + + radius = feature_ps.get_value(qualifier='radius', unit=u.rad, **_skip_filter_checks) + relteff = feature_ps.get_value(qualifier='relteff', unit=u.dimensionless_unscaled, **_skip_filter_checks) + + t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks) + + return cls(colat, longitude, dlongdt, radius, relteff, t0) + + @property + def _remeshing_required(self): + return False + + @property + def proto_coords(self): + """ + """ + return True + + def pointing_vector(self, s, time): + """ + s is the spin vector in roche coordinates + time is the current time + """ + t = time - self._t0 + longitude = self._longitude + self._dlongdt * t + + # define the basis vectors in the spin (primed) coordinates in terms of + # the Roche coordinates. + # ez' = s + # ex' = (ex - s(s.ex)) /|i - s(s.ex)| + # ey' = s x ex' + ex = np.array([1., 0., 0.]) + ezp = s + exp = (ex - s*np.dot(s,ex)) + eyp = np.cross(s, exp) + + return np.sin(self._colat)*np.cos(longitude)*exp +\ + np.sin(self._colat)*np.sin(longitude)*eyp +\ + np.cos(self._colat)*ezp + + def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): + """ + Change the local effective temperatures for any values within the + "cone" defined by the spot. Any teff within the spot will have its + current value multiplied by the "relteff" factor + + :parameter array teffs: array of teffs for computations + :parameter array coords: array of coords for computations + :t float: current time + """ + if t is None: + # then assume at t0 + t = self._t0 + + pointing_vector = self.pointing_vector(s,t) + logger.debug("spot.process_teffs at t={} with pointing_vector={} and radius={}".format(t, pointing_vector, self._radius)) + + cos_alpha_coords = np.dot(coords, pointing_vector) / np.linalg.norm(coords, axis=1) + cos_alpha_spot = np.cos(self._radius) + + filter_ = cos_alpha_coords > cos_alpha_spot + teffs[filter_] = teffs[filter_] * self._relteff + + return teffs + +class Pulsation(ComponentFeature): + def __init__(self, radamp, freq, l=0, m=0, tanamp=0.0, teffext=False, **kwargs): + self._freq = freq + self._radamp = radamp + self._l = l + self._m = m + self._tanamp = tanamp + + self._teffext = teffext + + @classmethod + def from_bundle(cls, b, feature): + """ + Initialize a Pulsation feature from the bundle. + """ + + feature_ps = b.get_feature(feature=feature, **_skip_filter_checks) + freq = feature_ps.get_value(qualifier='freq', unit=u.d**-1, **_skip_filter_checks) + radamp = feature_ps.get_value(qualifier='radamp', unit=u.dimensionless_unscaled, **_skip_filter_checks) + l = feature_ps.get_value(qualifier='l', unit=u.dimensionless_unscaled, **_skip_filter_checks) + m = feature_ps.get_value(qualifier='m', unit=u.dimensionless_unscaled, **_skip_filter_checks) + teffext = feature_ps.get_value(qualifier='teffext', **_skip_filter_checks) + + GM = c.G.to('solRad3 / (solMass d2)').value*b.get_value(qualifier='mass', component=feature_ps.component, context='component', unit=u.solMass, **_skip_filter_checks) + R = b.get_value(qualifier='rpole', component=feature_ps.component, section='component', unit=u.solRad, **_skip_filter_checks) + + tanamp = GM/R**3/freq**2 + + return cls(radamp, freq, l, m, tanamp, teffext) + + @property + def proto_coords(self): + """ + """ + return True + + def dYdtheta(self, m, l, theta, phi): + if abs(m) > l: + return 0 + + # TODO: just a quick hack + if abs(m+1) > l: + last_term = 0.0 + else: + last_term = Y(m+1, l, theta, phi) + + return m/np.tan(theta)*Y(m, l, theta, phi) + np.sqrt((l-m)*(l+m+1))*np.exp(-1j*phi)*last_term + + def dYdphi(self, m, l, theta, phi): + return 1j*m*Y(m, l, theta, phi) + + def process_coords_for_computations(self, coords_for_computations, s, t): + """ + """ + if self._teffext: + return coords_for_computations + + x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1)) + theta = np.arccos(z/r) + phi = np.arctan2(y, x) + + xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) + xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) + xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) + + new_coords = np.zeros(coords_for_computations.shape) + new_coords[:,0] = coords_for_computations[:,0] + xi_r * np.sin(theta) * np.cos(phi) + new_coords[:,1] = coords_for_computations[:,1] + xi_r * np.sin(theta) * np.sin(phi) + new_coords[:,2] = coords_for_computations[:,2] + xi_r * np.cos(theta) + + return new_coords + + def process_coords_for_observations(self, coords_for_computations, coords_for_observations, s, t): + """ + Displacement equations: + + xi_r(r, theta, phi) = a(r) Y_lm (theta, phi) exp(-i*2*pi*f*t) + xi_theta(r, theta, phi) = b(r) dY_lm/dtheta (theta, phi) exp(-i*2*pi*f*t) + xi_phi(r, theta, phi) = b(r)/sin(theta) dY_lm/dphi (theta, phi) exp(-i*2*pi*f*t) + + where: + + b(r) = a(r) GM/(R^3*f^2) + """ + # TODO: we do want to displace the coords_for_observations, but the x,y,z,r below are from the ALSO displaced coords_for_computations + # if not self._teffext: + # return coords_for_observations + + x, y, z, r = coords_for_computations[:,0], coords_for_computations[:,1], coords_for_computations[:,2], np.sqrt((coords_for_computations**2).sum(axis=1)) + theta = np.arccos(z/r) + phi = np.arctan2(y, x) + + xi_r = self._radamp * Y(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) + xi_t = self._tanamp * self.dYdtheta(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) + xi_p = self._tanamp/np.sin(theta) * self.dYdphi(self._m, self._l, theta, phi) * np.exp(-1j*2*np.pi*self._freq*t) + + new_coords = np.zeros(coords_for_observations.shape) + new_coords[:,0] = coords_for_observations[:,0] + xi_r * np.sin(theta) * np.cos(phi) + new_coords[:,1] = coords_for_observations[:,1] + xi_r * np.sin(theta) * np.sin(phi) + new_coords[:,2] = coords_for_observations[:,2] + xi_r * np.cos(theta) + + return new_coords + + def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None): + """ + """ + if not self._teffext: + return teffs + + raise NotImplementedError("teffext=True not yet supported for pulsations") diff --git a/phoebe/features/dataset_features.py b/phoebe/features/dataset_features.py new file mode 100644 index 000000000..2b70b1faa --- /dev/null +++ b/phoebe/features/dataset_features.py @@ -0,0 +1,58 @@ +import numpy as np +import astropy.units as u + +import logging +logger = logging.getLogger("DATASET_FEATURES") +logger.addHandler(logging.NullHandler()) + +from phoebe.parameters import FloatParameter, ParameterSet +from phoebe.parameters import constraint + +__all__ = ['DatasetFeature', 'SinusoidalThirdLight'] + +_skip_filter_checks = {'check_default': False, 'check_visible': False} + +class DatasetFeature(object): + """ + + """ + allowed_component_kinds = [None] + allowed_dataset_kinds = ['lc', 'rv', 'lp'] + def __init__(self, *args, **kwargs): + pass + + @classmethod + def from_bundle(cls, b, feature_ps): + return cls() + + @classmethod + def add_feature(self, feature, **kwargs): + raise NotImplementedError("add_feature must be implemented in the feature subclass") + + def modify_model(self, b, feature_ps, model_ps): + raise NotImplementedError("modify_model must be implemented in the feature subclass") + + + +class SinusoidalThirdLight(DatasetFeature): + allowed_dataset_kinds = ['lc'] + + @classmethod + def add_feature(self, feature, **kwargs): + params = [] + params += [FloatParameter(qualifier='amplitude', latexfmt=r'A_\mathrm{{ {feature} }}', value=kwargs.get('amplitude', 1.0), default_unit=u.W/u.m**2, description='Amplitude of the third light sinusoidal contribution')] + params += [FloatParameter(qualifier='period', latexfmt=r'P_\mathrm{{ {feature} }}', value=kwargs.get('period', 1.0), default_unit=u.d, description='Period of the third light sinusoidal contribution')] + params += [FloatParameter(qualifier='freq', latexfmt=r'f_\mathrm{{ {feature} }}', value=kwargs.get('freq', 2*np.pi/3.0), default_unit=u.rad/u.d, advanced=True, description='Orbital frequency (sidereal)')] + + constraints = [(constraint.freq, feature, 'feature')] + + return ParameterSet(params), constraints + + def modify_model(self, b, feature_ps, model_ps): + period = feature_ps.get_value(qualifier='period', context='system', unit=u.d, **_skip_filter_checks) + t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks) + for flux_param in model_ps.filter(qualifier='fluxes', **_skip_filter_checks).tolist(): + ampl = feature_ps.get_value(qualifier='amplitude', unit=flux_param.unit, **_skip_filter_checks) + times = model_ps.get_value(qualifier='times', dataset=flux_param.dataset, unit=u.d, **_skip_filter_checks) + flux_param.set_value(flux_param.get_value() + ampl * np.sin(2 * np.pi * (times - t0)) / period, ignore_readonly=True, **_skip_filter_checks) + diff --git a/phoebe/features/gaussian_processes.py b/phoebe/features/gaussian_processes.py new file mode 100644 index 000000000..22d0c9fe9 --- /dev/null +++ b/phoebe/features/gaussian_processes.py @@ -0,0 +1,164 @@ +from phoebe.parameters import * + +try: + import sklearn as _sklearn + from sklearn.gaussian_process import GaussianProcessRegressor +except ImportError: + logger.warning("scikit-learn not installed: only required for gaussian processes") + _use_sklearn = False +else: + _use_sklearn = True + +try: + import celerite2 as _celerite2 +except ImportError: + logger.warning("celerite2 not installed: only required for gaussian processes") + _use_celerite2 = False +else: + _use_celerite2 = True + +__all__ = ['handle_gaussian_processes'] + +_skip_filter_checks = {'check_default': False, 'check_visible': False} + +def handle_gaussian_processes(b, model, model_ps, enabled_features, computeparams): + + for ds in model_ps.datasets: + gp_sklearn_features = b.filter(feature=enabled_features, dataset=ds, kind='gp_sklearn', **_skip_filter_checks).features + gp_celerite2_features = b.filter(feature=enabled_features, dataset=ds, kind='gp_celerite2', **_skip_filter_checks).features + + if len(gp_sklearn_features)!=0 or len(gp_celerite2_features)!=0: + # we'll loop over components (for RVs or LPs, for example) + # get the data we need to fit the GP model + ds_ps = b.get_dataset(dataset=ds, **_skip_filter_checks) + xqualifier = {'lp': 'wavelength'}.get(ds_ps.kind, 'times') + yqualifier = {'lp': 'flux_densities', 'rv': 'rvs', 'lc': 'fluxes'}.get(ds_ps.kind) + yerrqualifier = {'lp': 'wavelength'}.get(ds_ps.kind, 'sigmas') + + _exclude_phases_enabled = computeparams.get_value(qualifier='gp_exclude_phases_enabled', dataset=ds, **_skip_filter_checks) + + if ds_ps.kind in ['lc']: + ds_comps = [None] + else: + ds_comps = ds_ps.filter(qualifier=xqualifier, check_visible=True).components + for ds_comp in ds_comps: + ds_x = ds_ps.get_value(qualifier=xqualifier, component=ds_comp, **_skip_filter_checks) + model_x = model_ps.get_value(qualifier=xqualifier, dataset=ds, component=ds_comp, **_skip_filter_checks) + ds_sigmas = ds_ps.get_value(qualifier=yerrqualifier, component=ds_comp, **_skip_filter_checks) + # ds_sigmas = ds_ps.get_value(qualifier='sigmas', component=ds_comp, **_skip_filter_checks) + # TODO: do we need to inflate sigmas by lnf? + if not len(ds_x): + # should have been caught by run_checks_compute + raise ValueError("gaussian_process requires dataset observations (cannot be synthetic only). Add observations to dataset='{}' or disable feature={}".format(ds, gp_features)) + + residuals, model_y_dstimes = b.calculate_residuals(model=model, + dataset=ds, + component=ds_comp, + return_interp_model=True, + as_quantity=False, + consider_gaussian_process=False) + model_y = model_ps.get_quantity(qualifier=yqualifier, dataset=ds, component=ds_comp, **_skip_filter_checks) + + gp_kernels = [] + alg_operations = [] + + def _load_gps(gp_kernel_classes, gp_features, ds): + + for gp in gp_features: + gp_ps = b.filter(feature=gp, context='feature', **_skip_filter_checks) + kind = gp_ps.get_value(qualifier='kernel', **_skip_filter_checks) + + kwargs = {p.qualifier: p.value for p in gp_ps.exclude(qualifier=['kernel', 'enabled']).to_list() if p.is_visible} + # TODO: replace this with getting the parameter from compute options + if _exclude_phases_enabled: + exclude_phase_ranges = computeparams.get_value(qualifier='gp_exclude_phases', dataset=ds, **_skip_filter_checks) + else: + exclude_phase_ranges = [] + + alg_operations.append(kwargs.pop('alg_operation')) + gp_kernels.append(gp_kernel_classes.get(kind)(**kwargs)) + + gp_kernel = gp_kernels[0] + for i in range(1, len(gp_kernels)): + if alg_operations[i] == 'product': + gp_kernel *= gp_kernels[i] + # print(gp_kernel) + else: + gp_kernel += gp_kernels[i] + # print(gp_kernel) + + + if len(exclude_phase_ranges) != 0: + # get t0, period and exclude_phases + ephem = b.get_ephemeris(component='binary', period='period', t0='t0_supconj') + t0 = ephem.get('t0', 0.0) + period = ephem.get('period', 1.0) + + phases = np.array(exclude_phase_ranges) + + # determine extent of data wrt t0 + i0 = int((t0 - min(ds_x))/period)-1 + i1 = int((max(ds_x-t0))/period)+1 + + x_new = ds_x + residuals_new = residuals + sigmas_new = ds_sigmas + for i in range(i0,i1+1,1): + for j in range(phases.shape[0]): + condition = (x_new < t0+(i+phases[j][0])*period) | (x_new > t0+(i+phases[j][1])*period) + x_new = x_new[condition] + residuals_new = residuals_new[condition] + sigmas_new = sigmas_new[condition] + + gp_x = x_new + gp_y = residuals_new + gp_yerr = sigmas_new + + else: + gp_x = ds_x + gp_y = residuals + gp_yerr = ds_sigmas + + return gp_kernel, gp_x, gp_y, gp_yerr + + if len(gp_sklearn_features) > 0: + gp_kernel_classes = {'constant': _sklearn.gaussian_process.kernels.ConstantKernel, + 'white': _sklearn.gaussian_process.kernels.WhiteKernel, + 'rbf': _sklearn.gaussian_process.kernels.RBF, + 'matern': _sklearn.gaussian_process.kernels.Matern, + 'rational_quadratic': _sklearn.gaussian_process.kernels.RationalQuadratic, + 'exp_sine_squared': _sklearn.gaussian_process.kernels.ExpSineSquared, + 'dot_product': _sklearn.gaussian_process.kernels.DotProduct} + + gp_kernel, gp_x, gp_y, gp_yerr = _load_gps(gp_kernel_classes, gp_sklearn_features, ds) + + gp_regressor = GaussianProcessRegressor(kernel=gp_kernel) + gp_regressor.fit(gp_x.reshape(-1,1), gp_y) + + # NOTE: .predict can also be called directly to the model times if we want to avoid interpolation altogether + gp_y = gp_regressor.predict(ds_x.reshape(-1,1), return_std=False) + + if len(gp_celerite2_features) > 0: + gp_kernel_classes = {'sho': _celerite2.terms.SHOTerm, + 'rotation': _celerite2.terms.RotationTerm, + 'matern32': _celerite2.terms.Matern32Term + } + gp_kernel, gp_x, gp_y, gp_yerr = _load_gps(gp_kernel_classes, gp_celerite2_features, ds) + + gp = _celerite2.GaussianProcess(gp_kernel, mean=0.0) + gp.compute(gp_x, yerr=gp_yerr) + gp_y = gp.predict(gp_y, t=ds_x, return_var=False) + + + # store just the GP component in the model PS as well + gp_param = FloatArrayParameter(qualifier='gps', value=gp_y, default_unit=model_y.unit, readonly=True, description='GP contribution to the model {}'.format(yqualifier)) + y_nogp_param = FloatArrayParameter(qualifier='{}_nogps'.format(yqualifier), value=model_y_dstimes, default_unit=model_y.unit, readonly=True, description='{} before adding gps'.format(yqualifier)) + if len(ds_x) != len(model_x) or not np.all(ds_x == model_x): + logger.warning("model for dataset='{}' resampled at dataset times when adding GPs".format(ds)) + model_ps.set_value(qualifier=xqualifier, dataset=ds, component=ds_comp, value=ds_x, ignore_readonly=True, **_skip_filter_checks) + + b._attach_params([gp_param, y_nogp_param], dataset=ds, check_copy_for=False, **metawargs) + + # update the model to include the GP contribution + model_ps.set_value(qualifier=yqualifier, value=model_y_dstimes+gp_y, dataset=ds, component=ds_comp, ignore_readonly=True, **_skip_filter_checks) + diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py index 149ed6168..c32557676 100644 --- a/phoebe/frontend/bundle.py +++ b/phoebe/frontend/bundle.py @@ -41,6 +41,8 @@ from phoebe.solverbackends import solverbackends as _solverbackends from phoebe.distortions import roche from phoebe.frontend import io +from phoebe.features import dataset_features +from phoebe.features.gaussian_processes import handle_gaussian_processes, _use_celerite2, _use_sklearn from phoebe.atmospheres.passbands import list_installed_passbands, list_online_passbands, get_passband, update_passband, _timestamp_to_dt from phoebe import pool as _pool from phoebe.dependencies import distl as _distl @@ -81,23 +83,6 @@ else: _can_client = True -try: - import sklearn as _sklearn - from sklearn.gaussian_process import GaussianProcessRegressor -except ImportError: - logger.warning("scikit-learn not installed: only required for gaussian processes") - _use_sklearn = False -else: - _use_sklearn = True - -try: - import celerite2 as _celerite2 -except ImportError: - logger.warning("celerite2 not installed: only required for gaussian processes") - _use_celerite2 = False -else: - _use_celerite2 = True - def _get_add_func(mod, func, return_none_if_not_found=False): if isinstance(func, str) and "." in func: @@ -12052,6 +12037,35 @@ def restore_conf(): flux_param.set_value(fluxes, ignore_readonly=True) + # handle vgamma and rv_offset + vgamma = self.get_value(qualifier='vgamma', context='system', unit=u.km/u.s, **_skip_filter_checks) + for rv_param in ml_params.filter(qualifier='rvs', kind='rv', **_skip_filter_checks).to_list(): + dataset = rv_param.dataset + component = rv_param.component + + rv_offset = self.get_value(qualifier='rv_offset', dataset=dataset, component=component, context='dataset', default=0.0, unit=u.km/u.s, **_skip_filter_checks) + + if computeparams.kind in ['phoebe', 'legacy']: + # we'll use native vgamma so ltte, etc, can be handled appropriately + rv_param.set_value(rv_param.get_value(unit=u.km/u.s)+rv_offset, ignore_readonly=True) + else: + rv_param.set_value(rv_param.get_value(unit=u.km/u.s)+vgamma+rv_offset, ignore_readonly=True) + + # handle all dataset-features (except for GPs, which need to all be handled simultaneously + # and AFTER - instead of BEFORE - dataset-scaling from pblum_mode='dataset-scaled') + enabled_features = self.filter(qualifier='enabled', compute=compute, context='compute', value=True, **_skip_filter_checks).features + for feature in enabled_features: + feature_ps = self.get_feature(feature=feature, **_skip_filter_checks) + # check first in globals()? + feature_cls = getattr(dataset_features, feature_ps.kind.title(), None) + if feature_cls is None: + # does not exist. Assume a component feature and skip. + # In the future if we add support for user-defined features + # in the namespace, we'll need to check for that as well and perhaps + # raise an error + continue + feature_cls.from_bundle(self, feature_ps).modify_model(self, feature_ps, ml_params) + # handle flux scaling for any pblum_mode == 'dataset-scaled' # or for any dataset in which pblum_mode == 'dataset-coupled' and pblum_dataset points to a 'dataset-scaled' dataset datasets_dsscaled = [] @@ -12157,20 +12171,6 @@ def _scale_fluxes_cfit(fluxes, scale_factor): flux_param.set_value(fluxes, ignore_readonly=True) - # handle vgamma and rv_offset - vgamma = self.get_value(qualifier='vgamma', context='system', unit=u.km/u.s, **_skip_filter_checks) - for rv_param in ml_params.filter(qualifier='rvs', kind='rv', **_skip_filter_checks).to_list(): - dataset = rv_param.dataset - component = rv_param.component - - rv_offset = self.get_value(qualifier='rv_offset', dataset=dataset, component=component, context='dataset', default=0.0, unit=u.km/u.s, **_skip_filter_checks) - - if computeparams.kind in ['phoebe', 'legacy']: - # we'll use native vgamma so ltte, etc, can be handled appropriately - rv_param.set_value(rv_param.get_value(unit=u.km/u.s)+rv_offset, ignore_readonly=True) - else: - rv_param.set_value(rv_param.get_value(unit=u.km/u.s)+vgamma+rv_offset, ignore_readonly=True) - ml_addl_params += [StringParameter(qualifier='comments', value=kwargs.get('comments', computeparams.get_value(qualifier='comments', default='', **_skip_filter_checks)), description='User-provided comments for this model. Feel free to place any notes here.')] self._attach_params(ml_params+ml_addl_params, check_copy_for=False, **metawargs) @@ -12178,147 +12178,9 @@ def _scale_fluxes_cfit(fluxes, scale_factor): # add any GPs (gaussian processes) to the returned model # NOTE: this has to happen after _attach_params as it uses - # several bundle methods that act on the model - enabled_features = self.filter(qualifier='enabled', compute=compute, context='compute', value=True, **_skip_filter_checks).features - - for ds in model_ps.datasets: - gp_sklearn_features = self.filter(feature=enabled_features, dataset=ds, kind='gp_sklearn', **_skip_filter_checks).features - gp_celerite2_features = self.filter(feature=enabled_features, dataset=ds, kind='gp_celerite2', **_skip_filter_checks).features - - if len(gp_sklearn_features)!=0 or len(gp_celerite2_features)!=0: - # we'll loop over components (for RVs or LPs, for example) - # get the data we need to fit the GP model - ds_ps = self.get_dataset(dataset=ds, **_skip_filter_checks) - xqualifier = {'lp': 'wavelength'}.get(ds_ps.kind, 'times') - yqualifier = {'lp': 'flux_densities', 'rv': 'rvs', 'lc': 'fluxes'}.get(ds_ps.kind) - yerrqualifier = {'lp': 'wavelength'}.get(ds_ps.kind, 'sigmas') - - _exclude_phases_enabled = computeparams.get_value(qualifier='gp_exclude_phases_enabled', dataset=ds, **_skip_filter_checks) - - if ds_ps.kind in ['lc']: - ds_comps = [None] - else: - ds_comps = ds_ps.filter(qualifier=xqualifier, check_visible=True).components - for ds_comp in ds_comps: - ds_x = ds_ps.get_value(qualifier=xqualifier, component=ds_comp, **_skip_filter_checks) - model_x = model_ps.get_value(qualifier=xqualifier, dataset=ds, component=ds_comp, **_skip_filter_checks) - ds_sigmas = ds_ps.get_value(qualifier=yerrqualifier, component=ds_comp, **_skip_filter_checks) - # ds_sigmas = ds_ps.get_value(qualifier='sigmas', component=ds_comp, **_skip_filter_checks) - # TODO: do we need to inflate sigmas by lnf? - if not len(ds_x): - # should have been caught by run_checks_compute - raise ValueError("gaussian_process requires dataset observations (cannot be synthetic only). Add observations to dataset='{}' or disable feature={}".format(ds, gp_features)) - - residuals, model_y_dstimes = self.calculate_residuals(model=model, - dataset=ds, - component=ds_comp, - return_interp_model=True, - as_quantity=False, - consider_gaussian_process=False) - model_y = model_ps.get_quantity(qualifier=yqualifier, dataset=ds, component=ds_comp, **_skip_filter_checks) - - gp_kernels = [] - alg_operations = [] - - def _load_gps(gp_kernel_classes, gp_features, ds): - - for gp in gp_features: - gp_ps = self.filter(feature=gp, context='feature', **_skip_filter_checks) - kind = gp_ps.get_value(qualifier='kernel', **_skip_filter_checks) - - kwargs = {p.qualifier: p.value for p in gp_ps.exclude(qualifier=['kernel', 'enabled']).to_list() if p.is_visible} - # TODO: replace this with getting the parameter from compute options - if _exclude_phases_enabled: - exclude_phase_ranges = computeparams.get_value(qualifier='gp_exclude_phases', dataset=ds, **_skip_filter_checks) - else: - exclude_phase_ranges = [] - - alg_operations.append(kwargs.pop('alg_operation')) - gp_kernels.append(gp_kernel_classes.get(kind)(**kwargs)) - - gp_kernel = gp_kernels[0] - for i in range(1, len(gp_kernels)): - if alg_operations[i] == 'product': - gp_kernel *= gp_kernels[i] - # print(gp_kernel) - else: - gp_kernel += gp_kernels[i] - # print(gp_kernel) - - - if len(exclude_phase_ranges) != 0: - # get t0, period and exclude_phases - ephem = self.get_ephemeris(component='binary', period='period', t0='t0_supconj') - t0 = ephem.get('t0', 0.0) - period = ephem.get('period', 1.0) - - phases = np.array(exclude_phase_ranges) - - # determine extent of data wrt t0 - i0 = int((t0 - min(ds_x))/period)-1 - i1 = int((max(ds_x-t0))/period)+1 - - x_new = ds_x - residuals_new = residuals - sigmas_new = ds_sigmas - for i in range(i0,i1+1,1): - for j in range(phases.shape[0]): - condition = (x_new < t0+(i+phases[j][0])*period) | (x_new > t0+(i+phases[j][1])*period) - x_new = x_new[condition] - residuals_new = residuals_new[condition] - sigmas_new = sigmas_new[condition] - - gp_x = x_new - gp_y = residuals_new - gp_yerr = sigmas_new - - else: - gp_x = ds_x - gp_y = residuals - gp_yerr = ds_sigmas - - return gp_kernel, gp_x, gp_y, gp_yerr - - if len(gp_sklearn_features) > 0: - gp_kernel_classes = {'constant': _sklearn.gaussian_process.kernels.ConstantKernel, - 'white': _sklearn.gaussian_process.kernels.WhiteKernel, - 'rbf': _sklearn.gaussian_process.kernels.RBF, - 'matern': _sklearn.gaussian_process.kernels.Matern, - 'rational_quadratic': _sklearn.gaussian_process.kernels.RationalQuadratic, - 'exp_sine_squared': _sklearn.gaussian_process.kernels.ExpSineSquared, - 'dot_product': _sklearn.gaussian_process.kernels.DotProduct} - - gp_kernel, gp_x, gp_y, gp_yerr = _load_gps(gp_kernel_classes, gp_sklearn_features, ds) - - gp_regressor = GaussianProcessRegressor(kernel=gp_kernel) - gp_regressor.fit(gp_x.reshape(-1,1), gp_y) - - # NOTE: .predict can also be called directly to the model times if we want to avoid interpolation altogether - gp_y = gp_regressor.predict(ds_x.reshape(-1,1), return_std=False) - - if len(gp_celerite2_features) > 0: - gp_kernel_classes = {'sho': _celerite2.terms.SHOTerm, - 'rotation': _celerite2.terms.RotationTerm, - 'matern32': _celerite2.terms.Matern32Term - } - gp_kernel, gp_x, gp_y, gp_yerr = _load_gps(gp_kernel_classes, gp_celerite2_features, ds) - - gp = _celerite2.GaussianProcess(gp_kernel, mean=0.0) - gp.compute(gp_x, yerr=gp_yerr) - gp_y = gp.predict(gp_y, t=ds_x, return_var=False) - - - # store just the GP component in the model PS as well - gp_param = FloatArrayParameter(qualifier='gps', value=gp_y, default_unit=model_y.unit, readonly=True, description='GP contribution to the model {}'.format(yqualifier)) - y_nogp_param = FloatArrayParameter(qualifier='{}_nogps'.format(yqualifier), value=model_y_dstimes, default_unit=model_y.unit, readonly=True, description='{} before adding gps'.format(yqualifier)) - if len(ds_x) != len(model_x) or not np.all(ds_x == model_x): - logger.warning("model for dataset='{}' resampled at dataset times when adding GPs".format(ds)) - model_ps.set_value(qualifier=xqualifier, dataset=ds, component=ds_comp, value=ds_x, ignore_readonly=True, **_skip_filter_checks) - - self._attach_params([gp_param, y_nogp_param], dataset=ds, check_copy_for=False, **metawargs) - - # update the model to include the GP contribution - model_ps.set_value(qualifier=yqualifier, value=model_y_dstimes+gp_y, dataset=ds, component=ds_comp, ignore_readonly=True, **_skip_filter_checks) + # several bundle methods that act on the model and also + # has to happen after dataset-scaling (for pblum_mode=='dataset-scaled') + handle_gaussian_processes(self, model, model_ps, enabled_features, computeparams) except Exception as err: restore_conf() diff --git a/phoebe/parameters/constraint.py b/phoebe/parameters/constraint.py index 3f2a1e141..3fdd8de0e 100644 --- a/phoebe/parameters/constraint.py +++ b/phoebe/parameters/constraint.py @@ -928,7 +928,7 @@ def ph_perpass(b, orbit, solve_for=None, **kwargs): _validsolvefor['freq'] = ['freq', 'period'] -def freq(b, component, solve_for=None, **kwargs): +def freq(b, component, context='component', solve_for=None, **kwargs): """ Create a constraint for frequency (either orbital or rotational) given a period. @@ -966,8 +966,8 @@ def freq(b, component, solve_for=None, **kwargs): -------- * NotImplementedError: if the value of `solve_for` is not implemented. """ - - component_ps = _get_system_ps(b, component) + # component could also be feature if context='feature' + component_ps = _get_system_ps(b, component, context=context) #metawargs = component_ps.meta #metawargs.pop('qualifier') @@ -986,7 +986,7 @@ def freq(b, component, solve_for=None, **kwargs): else: raise NotImplementedError - return lhs, rhs, [], {'component': component} + return lhs, rhs, [], {context: component} #} #{ Inter-orbit constraints diff --git a/phoebe/parameters/feature.py b/phoebe/parameters/feature.py index 909141ce3..9a9cdc610 100644 --- a/phoebe/parameters/feature.py +++ b/phoebe/parameters/feature.py @@ -2,6 +2,7 @@ from phoebe.parameters import * from phoebe.parameters import constraint +from phoebe.features.dataset_features import * from phoebe import u from phoebe import conf @@ -18,6 +19,7 @@ def _component_allowed_for_feature(feature_kind, component_kind): _allowed['gp_sklearn'] = [None] _allowed['gp_celerite2'] = [None] _allowed['gaussian_process'] = [None] # deprecated: remove in 2.5 + _allowed['sinusoidal_third_light'] = [None] return component_kind in _allowed[feature_kind] @@ -28,6 +30,7 @@ def _dataset_allowed_for_feature(feature_kind, dataset_kind): _allowed['gp_sklearn'] = ['lc', 'rv', 'lp'] _allowed['gp_celerite2'] = ['lc', 'rv', 'lp'] _allowed['gaussian_process'] = ['lc', 'rv', 'lp'] # deprecated: remove in 2.5 + _allowed['sinusoidal_third_light'] = ['lc'] return dataset_kind in _allowed[feature_kind] @@ -253,3 +256,21 @@ def gaussian_process(feature, **kwargs): """ logger.warning("gaussian_process is deprecated. Use gp_celerite2 instead.") return gp_celerite2(feature, **kwargs) + +def sinusoidal_third_light(feature, **kwargs): + """ + Add a sinusoidal third light contribution (additive) to the light curve. + + + Arguments + ---------- + * `amplitude` (float, optional, default=1.0): Amplitude of the sinusoidal term + * `period` (float, optional, default=1.0): Period of the sinusoidal term + + Returns + -------- + * (, list): ParameterSet of all newly created + objects and a list of all necessary + constraints. + """ + return SinusoidalThirdLight.add_feature(feature, **kwargs) From 5a88d218032839c0eb539158624eb4cd09bde75c Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Tue, 11 Jun 2024 06:02:59 -0400 Subject: [PATCH 2/3] build features --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index a82f71082..60f52b2dc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -98,6 +98,7 @@ packages = [ "phoebe.solverbackends", "phoebe.solverbackends.ebai", "phoebe.solverbackends.knn", + "phoebe.features", "phoebe.utils", "phoebe.helpers", "phoebe.pool", From ed40bafa49f6088c68a868d972dd501ace51c6a1 Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Tue, 11 Jun 2024 06:06:48 -0400 Subject: [PATCH 3/3] use class-defined allowed kinds --- phoebe/parameters/feature.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/phoebe/parameters/feature.py b/phoebe/parameters/feature.py index 9a9cdc610..fc05326ba 100644 --- a/phoebe/parameters/feature.py +++ b/phoebe/parameters/feature.py @@ -19,7 +19,7 @@ def _component_allowed_for_feature(feature_kind, component_kind): _allowed['gp_sklearn'] = [None] _allowed['gp_celerite2'] = [None] _allowed['gaussian_process'] = [None] # deprecated: remove in 2.5 - _allowed['sinusoidal_third_light'] = [None] + _allowed['sinusoidal_third_light'] = SinusoidalThirdLight.allowed_component_kinds return component_kind in _allowed[feature_kind] @@ -30,7 +30,7 @@ def _dataset_allowed_for_feature(feature_kind, dataset_kind): _allowed['gp_sklearn'] = ['lc', 'rv', 'lp'] _allowed['gp_celerite2'] = ['lc', 'rv', 'lp'] _allowed['gaussian_process'] = ['lc', 'rv', 'lp'] # deprecated: remove in 2.5 - _allowed['sinusoidal_third_light'] = ['lc'] + _allowed['sinusoidal_third_light'] = SinusoidalThirdLight.allowed_dataset_kinds return dataset_kind in _allowed[feature_kind]