Merge pull request #702 from kevin218/catwoman

Adding catwoman_tr to model transit limb-asymmetries
kevin218 · Nov 7, 2024 · 7483bc0 · 7483bc0
2 parents fde2b69 + 05aa571
commit 7483bc0
Show file tree

Hide file tree

Showing 15 changed files with 144 additions and 13 deletions.
diff --git a/demos/HST/S5_wfc3_template.ecf b/demos/HST/S5_wfc3_template.ecf
@@ -10,7 +10,7 @@ allapers        False  # Run S5 on all of the apertures considered in S4? Otherw
 
 fit_par         ./S5_fit_par_wfc3_template.epf  # What fitting epf do you want to use?
 fit_method      [dynesty]              #options are: lsq, emcee, dynesty (can list multiple types separated by commas)
-run_myfuncs     [batman_tr,polynomial,hstramp] # options are: batman_tr, batman_ecl, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
+run_myfuncs     [batman_tr,polynomial,hstramp] # options are: batman_tr, batman_ecl, catwoman_tr, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
 compute_ltt     False                   # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect)
 
 # Manual clipping in time

diff --git a/demos/JWST/S5_template.ecf b/demos/JWST/S5_template.ecf
@@ -13,7 +13,7 @@ manual_clip     None    # A list of lists specifying the start and end integrati
 
 fit_par         ./S5_fit_par_template.epf   # What fitting epf do you want to use?
 fit_method      [dynesty]                   #options are: lsq, emcee, dynesty (can list multiple types separated by commas)
-run_myfuncs     [batman_tr,polynomial]      # options are: batman_tr, batman_ecl, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
+run_myfuncs     [batman_tr,polynomial]      # options are: batman_tr, batman_ecl, catwoman_tr, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
 compute_ltt     False                       # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect)
 force_positivity False                      # Optional boolean for sinusoid_pc and poet_pc models.  Set True to force positive phase variations.
 

diff --git a/docs/media/S5_template.ecf b/docs/media/S5_template.ecf
@@ -13,7 +13,7 @@ manual_clip     None    # A list of lists specifying the start and end integrati
 
 fit_par         ./S5_fit_par_template.epf   # What fitting epf do you want to use?
 fit_method      [dynesty]                   #options are: lsq, emcee, dynesty (can list multiple types separated by commas)
-run_myfuncs     [batman_tr,polynomial]      # options are: batman_tr, batman_ecl, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
+run_myfuncs     [batman_tr,polynomial]      # options are: batman_tr, batman_ecl, catwoman_tr, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
 compute_ltt     False                       # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect)
 force_positivity False                      # Optional boolean for sinusoid_pc and poet_pc models.  Set True to force positive phase variations.
 

diff --git a/docs/source/ecf.rst b/docs/source/ecf.rst
@@ -887,7 +887,7 @@ run_myfuncs
 '''''''''''
 Determines the astrophysical and systematics models used in the Stage 5 fitting.
 For standard numpy functions, this can be one or more (separated by commas) of the following:
-[batman_tr, batman_ecl, fleck_tr, poet_tr, poet_ecl, sinusoid_pc, quasilambert_pc, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, GP].
+[batman_tr, batman_ecl, catwoman_tr, fleck_tr, poet_tr, poet_ecl, sinusoid_pc, quasilambert_pc, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, GP].
 For theano-based differentiable functions, this can be one or more of the following:
 [starry, sinusoid_pc, quasilambert_pc, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth],
 where starry replaces both the batman_tr and batman_ecl models and offers a more complicated phase variation model than sinusoid_pc that accounts for eclipse mapping signals.
@@ -915,6 +915,17 @@ manual_clip
 '''''''''''
 Optional. A list of lists specifying the start and end integration numbers for manual removal. E.g., to remove the first 20 data points specify [[0,20]], and to also remove the last 20 data points specify [[0,20],[-20,None]]. If you want to clip the 10th integration, this would be index 9 since python uses zero-indexing. And the manual_clip start and end values are used to slice a numpy array, so they follow the same convention of *inclusive* start index and *exclusive* end index. In other words, to trim the 10th integrations, you would set manual_clip to [[9,10]].
 
+Catwoman Convergence Parameters
+'''''''''''''''''''''''''''''''
+The following two parameters can help in the case of convergence issues when using catwoman_tr.
+
+catwoman_max_err
+^^^^^^^^^^^^^^^^
+The ``max_err`` parameter used by catwoman; defaults to ``1.0``. For more information, see the relevant location of catwoman's `readthedocs page <https://catwoman.readthedocs.io/en/latest/tutorial.html#error-tolerance>`_ or the relevant location of `catwoman's API <https://catwoman.readthedocs.io/en/latest/API.html#catwoman.TransitModel>`_.
+
+catwoman_fac
+^^^^^^^^^^^^
+The ``fac`` parameter used by catwoman; defaults to ``None``. For more information, see the relevant location of catwoman's `readthedocs page <https://catwoman.readthedocs.io/en/latest/tutorial.html#error-tolerance>`_ or the relevant location of `catwoman's API <https://catwoman.readthedocs.io/en/latest/API.html#catwoman.TransitModel>`_.
 
 Limb Darkening Parameters
 '''''''''''''''''''''''''
@@ -1091,6 +1102,8 @@ Available fitting parameters are:
    - Transit and Eclipse Depth Parameters
       - ``rp`` or ``rprs`` - planet-to-star radius ratio, for the transit models.
       - ``fp`` or ``fpfs`` - planet-to-star flux ratio, for the eclipse models.
+      - ``rp2`` or ``rprs2`` - an additional planet-to-star radius ratio for use with the catwoman transit model to model transit limb-asymmetries.
+      - ``phi`` - the angle (in degrees) of the line separating the semi-circles defined by ``rp`` and ``rp2`` in the catwoman transit model. If ``phi`` is set to 90 degrees (the parameter's default value), the ``rp`` is the trailing hemisphere and ``rp2`` is the leading hemisphere. If ``phi`` is set to 0, then ``rp`` is the northern hemisphere and ``rp2`` is the southern hemisphere.
       When fitting for multiple planets, add ``_pl#`` after the parameter (e.g., ``rprs``, ``rprs_pl1``, ``rprs_pl2``, etc.). This also applies to the planetaty orbital parameters below. Also be sure to set the ``num_planets`` parameter in your ECF (not EPF) to specify the number of planets being modelled simultaneously.
    - Orbital Parameters
       - ``per`` - orbital period (in days)

diff --git a/environment.yml b/environment.yml
@@ -36,6 +36,7 @@ dependencies:
     - tqdm
     - pip:
         - astraeus@git+https://github.com/kevin218/Astraeus@main
+        - catwoman@git+https://github.com/taylorbell57/catwoman@master
         - celerite2 # Needed for GP
         - crds<12 # Upper limit needed to avoid bugs with crds.get_context_name
         - exotic-ld==3.2.0 # Lower limit needed for updated JWST sensitivity files, upper limit needed for breaking changes

diff --git a/environment_pymc3.yml b/environment_pymc3.yml
@@ -39,6 +39,7 @@ dependencies:
     - pip:
         - arviz==0.12.1
         - astraeus@git+https://github.com/kevin218/Astraeus@main
+        - catwoman@git+https://github.com/taylorbell57/catwoman@master
         - celerite2 # Needed for GP
         - crds<12 # Upper limit needed to avoid bugs with crds.get_context_name
         - exoplanet

diff --git a/setup.cfg b/setup.cfg
@@ -29,6 +29,7 @@ install_requires =
     astroquery
     batman-package
     bokeh<3.0
+    catwoman@git+https://github.com/taylorbell57/catwoman@master
     ccdproc
     celerite2 # Needed for GP
     corner

diff --git a/src/eureka/S5_lightcurve_fitting/models/AstroModel.py b/src/eureka/S5_lightcurve_fitting/models/AstroModel.py
@@ -68,6 +68,9 @@ def __init__(self, model, pid=0, channel=0, eval=True):
         self.t0 = None
         self.rprs = None
         self.rp = None
+        self.rprs2 = None
+        self.rp2 = None
+        self.phi = 90.
         self.inc = None
         self.ars = None
         self.a = None
@@ -158,6 +161,23 @@ def __init__(self, model, pid=0, channel=0, eval=True):
             if eval:
                 value = value.value
             setattr(self, 'rp', value)
+        # Allow for rp2 or rprs2
+        if (self.rprs2 is None) and ('rp2' in model.parameters.dict.keys()):
+            item0 = 'rp2' + self.pid_id
+            if model.parameters.dict[item0][1] == 'free':
+                item0 += self.channel_id
+            value = getattr(parameterObject, item0)
+            if eval:
+                value = value.value
+            setattr(self, 'rprs2', value)
+        if (self.rp2 is None) and ('rprs2' in model.parameters.dict.keys()):
+            item0 = 'rprs2' + self.pid_id
+            if model.parameters.dict[item0][1] == 'free':
+                item0 += self.channel_id
+            value = getattr(parameterObject, item0)
+            if eval:
+                value = value.value
+            setattr(self, 'rp2', value)
         # Allow for a or ars
         if (self.ars is None) and ('a' in model.parameters.dict.keys()):
             item0 = 'a' + self.pid_id

diff --git a/src/eureka/S5_lightcurve_fitting/models/CatwomanModel.py b/src/eureka/S5_lightcurve_fitting/models/CatwomanModel.py
@@ -0,0 +1,34 @@
+from functools import partial
+try:
+    import catwoman
+except ImportError:
+    print("Could not import catwoman. Functionality may be limited.")
+
+from .BatmanModels import BatmanTransitModel
+
+
+class CatwomanTransitModel(BatmanTransitModel):
+    """Transit Model"""
+    def __init__(self, **kwargs):
+        """Initialize the transit model
+
+        Parameters
+        ----------
+        **kwargs : dict
+            Additional parameters to pass to
+            eureka.S5_lightcurve_fitting.models.Model.__init__().
+            Can pass in the parameters, longparamlist, nchan, and
+            paramtitles arguments here.
+        """
+        # Inherit from BatmanTransitModel class
+        super().__init__(**kwargs)
+        self.name = 'catwoman transit'
+        # Define transit model to be used
+        self.transit_model = partial(catwoman.TransitModel,
+                                     max_err=kwargs['max_err'],
+                                     fac=kwargs['fac'])
+
+        if ('rp2' not in self.longparamlist[0]
+                and 'rprs2' not in self.longparamlist[0]):
+            raise AssertionError('You must include an rp2 parameter in your '
+                                 'EPF when using catwoman.')
diff --git a/src/eureka/S5_lightcurve_fitting/models/__init__.py b/src/eureka/S5_lightcurve_fitting/models/__init__.py
@@ -1,6 +1,7 @@
 from .Model import Model, CompositeModel
 from .AstroModel import AstroModel, PlanetParams
 from .BatmanModels import BatmanTransitModel, BatmanEclipseModel
+from .CatwomanModel import CatwomanTransitModel
 from .FleckModel import FleckTransitModel
 from .CentroidModel import CentroidModel
 from .DampedOscillator import DampedOscillatorModel

diff --git a/src/eureka/S5_lightcurve_fitting/s5_fit.py b/src/eureka/S5_lightcurve_fitting/s5_fit.py
@@ -569,6 +569,7 @@ def fit_channel(meta, time, flux, chan, flux_err, eventlabel, params,
     else:
         BatmanTransitModel = m.BatmanTransitModel
         BatmanEclipseModel = m.BatmanEclipseModel
+        CatwomanTransitModel = m.CatwomanTransitModel
         PoetTransitModel = m.PoetTransitModel
         PoetEclipseModel = m.PoetEclipseModel
         PoetPCModel = m.PoetPCModel
@@ -670,6 +671,27 @@ def fit_channel(meta, time, flux, chan, flux_err, eventlabel, params,
                                        nints=lc_model.nints,
                                        num_planets=meta.num_planets)
         modellist.append(t_eclipse)
+    if 'catwoman_tr' in meta.run_myfuncs:
+        t_transit = CatwomanTransitModel(parameters=params,
+                                         fmt='r--', log=log, time=time,
+                                         time_units=time_units,
+                                         freenames=freenames,
+                                         longparamlist=lc_model.longparamlist,
+                                         nchannel=chanrng,
+                                         nchannel_fitted=nchannel_fitted,
+                                         fitted_channels=fitted_channels,
+                                         paramtitles=paramtitles,
+                                         ld_from_S4=meta.use_generate_ld,
+                                         ld_from_file=meta.ld_file,
+                                         ld_coeffs=ldcoeffs,
+                                         recenter_ld_prior=meta.recenter_ld_prior,  # noqa: E501
+                                         compute_ltt=meta.compute_ltt,
+                                         multwhite=lc_model.multwhite,
+                                         nints=lc_model.nints,
+                                         num_planets=meta.num_planets,
+                                         fac=meta.catwoman_fac,
+                                         max_err=meta.catwoman_max_err)
+        modellist.append(t_transit)
     if 'fleck_tr' in meta.run_myfuncs:
         t_transit = m.FleckTransitModel(parameters=params,
                                         fmt='r--', log=log, time=time,

diff --git a/src/eureka/S5_lightcurve_fitting/s5_meta.py b/src/eureka/S5_lightcurve_fitting/s5_meta.py
@@ -110,6 +110,10 @@ def set_defaults(self):
             # Set this to False if not relevant
             self.recenter_ld_prior = getattr(self, 'recenter_ld_prior', False)
 
+        # Catwoman convergence-aiding parameters
+        self.catwoman_fac = getattr(self, 'catwoman_fac', None)
+        self.catwoman_max_err = getattr(self, 'catwoman_max_err', 1.0)
+
         # General fitter, fitparams CSV file to resume from
         self.old_fitparams = getattr(self, 'old_fitparams', None)
 

diff --git a/src/eureka/S6_planet_spectra/s6_spectra.py b/src/eureka/S6_planet_spectra/s6_spectra.py
@@ -224,7 +224,8 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None):
                     continue
 
                 # Manipulate fitted values if needed
-                if meta.y_param_basic in ['rp^2', 'rprs^2']:
+                if (meta.y_param_basic[:2] == 'rp' and
+                        meta.y_param_basic[-2:] == '^2'):
                     meta = compute_transit_depth(meta)
                 elif meta.y_param_basic in ['1/r1', '1/r3']:
                     meta = compute_timescale(meta)
@@ -234,11 +235,12 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None):
 
                 if meta.y_label is None:
                     # Provide some default formatting
-                    if meta.y_param_basic in ['rp^2', 'rprs^2']:
+                    if (meta.y_param_basic[:2] == 'rp' and
+                            meta.y_param_basic[-2:] == '^2'):
                         # Transit depth
                         suffix = planetSuffix+channelSuffix
                         meta.y_label = '$(R_{\\rm p'+suffix+'}/R_{\\rm *})^2$'
-                    elif meta.y_param_basic in ['rp', 'rprs']:
+                    elif meta.y_param_basic[:2] == 'rp':
                         # Radius ratio
                         suffix = planetSuffix+channelSuffix
                         meta.y_label = '$R_{\\rm p'+suffix+'}/R_{\\rm *}$'
@@ -361,8 +363,7 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None):
 
                 # Should we also make the scale_height version of the figure?
                 make_fig6301 = (meta.isplots_S6 >= 3 and meta.has_fig6301reqs
-                                and meta.y_param in ['rp', 'rp^2',
-                                                     'rprs', 'rprs^2'])
+                                and meta.y_param[:2] == 'rp')
                 if make_fig6301:
                     # Make spectrum plot with scale height on the 2nd y-axis
                     scale_height = compute_scale_height(meta, log)
@@ -414,7 +415,8 @@ def parse_s5_saves(meta, log, fit_methods, channel_key='shared'):
     errs
         The uncertainties from a sampling algorithm like dynesty or emcee.
     """
-    if meta.y_param_basic in ['rp^2', 'rprs^2']:
+    if (meta.y_param_basic[:2] == 'rp' and
+            meta.y_param_basic[-2:] == '^2'):
         y_param = meta.y_param[:-2]
     elif meta.y_param_basic in ['1/r1', '1/r3']:
         y_param = meta.y_param[2:]
@@ -1567,7 +1569,8 @@ def compute_scale_height(meta, log):
     """
     if meta.planet_Rad is None:
         meta.planet_Rad = meta.spectrum_median
-        if meta.y_param_basic in ['rp^2', 'rprs^2']:
+        if (meta.y_param_basic[:2] == 'rp' and
+                meta.y_param_basic[-2:] == '^2'):
             meta.planet_Rad = np.sqrt(meta.planet_Rad)
         meta.planet_Rad = np.nanmean(meta.planet_Rad)
         meta.planet_Rad *= (meta.star_Rad*constants.R_sun /
@@ -1829,7 +1832,8 @@ def transit_latex_table(meta, log):
     rows = int(np.ceil(nvals/meta.ncols))
 
     # Figure out the labels for the columns
-    if meta.y_param_basic in ['rp^2', 'rprs^2']:
+    if (meta.y_param_basic[:2] == 'rp' and
+            meta.y_param_basic[-2:] == '^2'):
         suffix = getPlanetSuffix(meta)+getChannelSuffix(meta)
         colhead = '\\colhead{Transit Depth'+suffix+'}'
     elif meta.y_param_basic in ['rp', 'rprs']:
@@ -1851,7 +1855,7 @@ def transit_latex_table(meta, log):
         out += "CC|"
     out = out[:-1]+"}\n"
     # Give the table a caption based on the tabulated data
-    if meta.y_param_basic in ['rp', 'rp^2', 'rprs', 'rprs^2']:
+    if meta.y_param_basic[:2] == 'rp':
         out += "\\tablecaption{\\texttt{Eureka!}'s Transit Spectroscopy "
         out += "Results \\label{tab:eureka_transit_spectra}}\n"
     elif meta.y_param_basic in ['fp', 'fpfs']:

diff --git a/src/eureka/lib/citations.py b/src/eureka/lib/citations.py
@@ -735,6 +735,34 @@
                     adsnote = {Provided by the SAO/NASA Astrophysics Data System}
                 }"""
                ],
+    "catwoman": [r"""@ARTICLE{catwoman:joss,
+                    author = {{Jones}, Kathryn and {Espinoza}, N{\'e}stor},
+                    title = "{catwoman: A transit modelling Python package for asymmetric light curves}",
+                    journal = {The Journal of Open Source Software},
+                    keywords = {C, exoplanets, Python, transit, C++, astronomy, Astrophysics - Instrumentation and Methods for Astrophysics, Astrophysics - Earth and Planetary Astrophysics},
+                    year = 2020,
+                    month = nov,
+                    volume = {5},
+                    number = {55},
+                    eid = {2382},
+                    pages = {2382},
+                    doi = {10.21105/joss.02382},
+                    archivePrefix = {arXiv},
+                    eprint = {2106.15643},
+                    primaryClass = {astro-ph.IM},
+                    adsurl = {https://ui.adsabs.harvard.edu/abs/2020JOSS....5.2382J},
+                    adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+                  }""",
+                 r"""@software{catwoman:ascl,
+                    author = {{Jones}, Kathryn and {Espinoza}, N{\'e}stor},
+                    title = "{catwoman: Transit modeling Python package for asymmetric light curves}",
+                    howpublished = {Astrophysics Source Code Library, record ascl:2108.007},
+                    year = 2021,
+                    month = aug,
+                    eid = {ascl:2108.007},
+                    adsurl = {https://ui.adsabs.harvard.edu/abs/2021ascl.soft08007J},
+                    adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+                  }"""],
     "fleck": [r"""@software{2022ascl.soft03009M,
                 author = {{Morris}, Brett M.},
                 title = "{fleck: Fast starspot rotational modulation light curves}",

diff --git a/src/eureka/lib/util.py b/src/eureka/lib/util.py
@@ -895,6 +895,8 @@ def make_citations(meta, stage=None):
         # check if batman or GP is being used for transit/eclipse modeling
         if "batman_tr" in meta.run_myfuncs or "batman_ecl" in meta.run_myfuncs:
             other_cites.append("batman")
+        if "catwoman_tr" in meta.run_myfuncs:
+            other_cites.append("catwoman")
         if "fleck_tr" in meta.run_myfuncs:
             other_cites.append("fleck")
         if "starry" in meta.run_myfuncs: