diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 1f35376..733b1b4 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -52,9 +52,10 @@ jobs: flag-name: Unit Test upload_pypi: - if: github.event_name == 'push' + if: startsWith(github.ref, 'refs/tags/v') runs-on: ubuntu-latest permissions: + contents: read id-token: write # or, alternatively, upload to PyPI on every tag starting with 'v' (remove on: release above to use this) # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v') @@ -72,13 +73,6 @@ jobs: pip install -r requirements.txt pip install -e . python setup.py sdist bdist_wheel - - - uses: actions/download-artifact@v4 - with: - # unpacks all CIBW artifacts into dist/ - path: dist - merge-multiple: true - - uses: pypa/gh-action-pypi-publish@release/v1 with: repository-url: https://test.pypi.org/legacy/ diff --git a/.markdownlint.json b/.markdownlint.json new file mode 100644 index 0000000..286140e --- /dev/null +++ b/.markdownlint.json @@ -0,0 +1,14 @@ +{ + "default": true, + "MD013": false, + "MD022": false, + "MD024": false, + "MD032": false, + "MD041": false, + "MD043": { + "headings": [ + "# CHANGELOG", + "+" + ] + } +} diff --git a/.readthedocs.yml b/.readthedocs.yml index 1620807..5f924bd 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -16,6 +16,7 @@ build: sphinx: builder: html configuration: docs/source/conf.py + fail_on_warning: true # If using Sphinx, optionally build your docs in additional formats such as PDF # formats: @@ -24,4 +25,6 @@ sphinx: # Optionally declare the Python requirements required to build your docs python: install: - - requirements: docs/source/requirements.txt + - requirements: docs/source/requirements.txt + - method: pip + path: . diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..f1fe39f --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,137 @@ +# CHANGELOG + +All notable changes to WDPhotTools are documented in this file. + +The 0.0.x line is treated as the beta series. Entries below are listed newest +to oldest and aligned to the 10 beta tags (`0.0.4` to `v0.0.13`). + +## [Unreleased] - 2026-07-07 + +- New feature: Added canonical photometry API support in `WDfitter.fit` with + `photometry`, `photometry_errors`, and `photometry_space`. +- New feature: Added fitting in flux space alongside magnitude space. +- New feature: Added `self.best_fit_photometry` in the selected + `photometry_space`. +- New feature: Added opt-in interpolator extrapolation controls for atmosphere + and cooling model readers (default remains disabled). +- New feature: Added sanitisation for extrapolated values to avoid aphysical + outputs (`NaN`, `inf`, and impossible negatives) in affected paths. +- Behaviour: Kept in-range interpolation and boundary behaviour unchanged when + extrapolation is disabled. +- API change (breaking): Removed legacy fit inputs `mags`, `mag_errors`, + `fluxes`, and `flux_errors`. +- API change (breaking): Standardised `self.fitting_params` to only + `photometry`, `photometry_errors`, and `photometry_space`. +- API change (breaking): Removed legacy output aliases `best_fit_mag` and + `best_fit_flux`. +- Logging: Replaced touched `print` calls with structured `logging`. +- Logging: Added INFO/WARNING logs at fitter orchestration boundaries. +- Documentation: Updated README, examples, and docstrings to the canonical + photometry API. +- Documentation: Updated RTD configuration/build and migration guidance. +- Documentation: Added extrapolation guidance and a success-rate table across + 10% to 50% extrapolation levels. +- Tests: Updated fitter tests to canonical API calls. +- Tests: Added deterministic flux-vs-magnitude parity tests. +- Tests: Added deterministic regression checks against `origin/main`. +- Tests: Reduced runtime for the new parity/regression matrix. +- Tests: Added extrapolation safety and boundary-behaviour coverage for + atmosphere/cooling interpolators. +- Examples: Updated example scripts for the new API/workflow expectations and + checked near-identical behaviour against `v0.0.13` where deterministic. +- Documentation: Added interpolation scheme explanation + +## [v0.0.13] - 2026-01-14 + +- Release: Versioned release `v0.0.13` (tag commit `e81da92`). + +## [v0.0.12] - 2025-12-27 + +- Fix: Corrected conversion constant usage in fitter-related numerical paths. +- Fix: Ensured `_integrand` returns `float` for improved compatibility. +- Compatibility: Improved posture for the NumPy 1/2 transition. +- Compatibility: Marked Windows with Python 3.13/3.14 CI as allowed failures + during transition. +- CI/tooling: Updated GitHub Actions workflow and checkout handling. +- CI/tooling: Updated setup metadata and test-pypi workflow on `dev`. +- CI/tooling: Performed general housekeeping before release. + +## [v0.0.11] - 2025-09-22 + +- Performance: Significantly improved WDLF computation runtime. +- Performance: Refactored basis interpolation and integration helper paths. +- Fix: Corrected wrong-array update bug and additional minor defects. +- Fix: Cleaned redundant `.keys()` usage. +- Dependency/packaging: Replaced `pkg_resources` with `importlib`. +- Documentation/tooling: Updated RTD configuration/docs plus formatting and CI. +- Tests: Updated tests to use `Agg` backend where needed. +- Tests: Fixed test filename issues. + +## [v0.0.10] - 2025-05-13 + +- New feature: Added support for user-provided priors in fitting workflows. +- Validation/fix: Added and updated tests for prior handling. +- Validation/fix: Fixed missing `z_min` and `z_max` handling in reddening + cases. +- Validation/fix: Decoupled objective functions from fitter orchestration. +- Validation/fix: Avoided interpolation-grid breakage in log-scale age grids. +- CI/tooling: Updated test environments and weekly build checks. +- CI/tooling: Updated pre-commit line-length configuration to 120 chars. +- CI/tooling: Applied multiple housekeeping, style, and markdown fixes. + +## [v0.0.9] - 2023-08-11 + +- Fix: Corrected example script syntax issues. +- Fix: Corrected interpolation variable-allowlist omissions. +- Fix: Corrected plotting type and element-wise comparison edge cases. +- Fix: Added `> 0` safety checks in density computation paths. + +## [v0.0.8] - 2023-05-07 + +- New feature: Added and validated fitter mass-estimation test coverage. +- Fix: Corrected fitting-mass bug. +- Fix: Corrected independent-variable name comparison robustness. +- Fix: Migrated from deprecated `interp2d` to RBF interpolation. +- Fix: Standardised fitter method outputs to scalar floats. +- Documentation: Updated README and RTD content plus citation references. +- Documentation: Added clarifications on uncertainty/error estimation behaviour. + +## [v0.0.7] - 2022-12-04 + +- New feature: Stored `number_density` as object properties before CSV export. +- Fix: Corrected sign bug in `dL/dt` number-density computation. +- Fix: Added protections for log-scale y-axis with non-positive data. +- Fix: Corrected negative-zero and additional plotting/cleanup issues. +- Fix: Corrected default `Rv` from `0.0` to `3.1`. +- CI/tooling: Added Python 3.11 in tests. +- CI/tooling: Added code-analysis workflow and broad tidying. + +## [v0.0.6] - 2022-10-24 + +- New feature: Added linear fractional extinction model support. +- New feature: Added RA/Dec-aware fitter updates for extinction workflows. +- New feature: Extended reddening API for linearly interpolated extinction. +- Compatibility: Lowered minimum `astropy` version to support Python 3.7. +- Compatibility: Added Python 3.10 and Windows CI/test coverage. +- Tests: Extended tests for new reddening/extinction paths. +- Tests: Reduced MCMC test walkers/steps/burn-in and relaxed tolerance for + stability. +- Fix: Corrected dependent-variable distance handling edge cases. + +## [v0.0.5] - 2022-09-12 + +- Refactor: Moved package code to `src/` layout. +- Refactor: Removed `autograd` usage. +- Fix: Corrected minimization-function bug(s). +- Fix: Updated outdated examples. +- Documentation: Added docs and README notes for retrieving fitted solutions. +- Documentation: Added uncertainty notes and plotting updates. + +## [v0.0.4] - 2022-08-09 + +- Compatibility: Adapted to SciPy 1.9 stricter data-type requirements. +- Compatibility: Added SciPy-version handling for RegularGridInterpolator + behaviour. +- Fix: Corrected NaN-extinction behaviour when reddening is not fitted. +- Fix: Corrected K01 IMF normalization bug. +- Fix: Corrected README/docs typos and fitter docstrings. diff --git a/README.md b/README.md index 4d9c9bf..81177e0 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,9 @@ For best performance with RBFInterpolator, use SciPy 1.9+. Documentation and more examples can be found at [Read the Docs](https://wdphottools.readthedocs.io/en/latest/). +For extrapolation behaviour and success-rate estimates from 10%-50% beyond the +grid span, see the +[interpolator extrapolation table](https://wdphottools.readthedocs.io/en/latest/background/extrapolation.html). ## Attribution @@ -241,8 +244,8 @@ ftr = WDfitter() ftr.fit( atmosphere="H", filters=["g_ps1", "r_ps1", "i_ps1", "z_ps1", "y_ps1", "G3", "G3_BP", "G3_RP", "J_mko", "H_mko", "K_mko"], - mags=[21.1437, 19.9678, 19.4993, 19.2981, 19.1478, 20.0533, 20.7883, 19.1868, 19.45-0.91, 19.96-1.39, 20.40-1.85], - mag_errors=[0.0321, 0.0229, 0.0083, 0.0234, 0.0187, 0.006322, 0.118615, 0.070880, 0.05, 0.03, 0.05], + photometry=[21.1437, 19.9678, 19.4993, 19.2981, 19.1478, 20.0533, 20.7883, 19.1868, 19.45-0.91, 19.96-1.39, 20.40-1.85], + photometry_errors=[0.0321, 0.0229, 0.0083, 0.0234, 0.0187, 0.006322, 0.118615, 0.070880, 0.05, 0.03, 0.05], independent=["Teff", "logg"], initial_guess=[4000.0, 7.5], distance=71.231, @@ -264,6 +267,46 @@ ftr.show_corner_plot( ) ``` +### Fitting in relative flux space + +`WDfitter.fit` also supports relative flux fitting through `photometry_space="flux"`. + +```python +import numpy as np +from WDPhotTools.fitter import WDfitter + +filters = ["G3", "G3_BP", "G3_RP"] +magnitudes = np.array([10.882, 10.853, 10.946], dtype=float) +magnitude_errors = np.array([0.02, 0.02, 0.02], dtype=float) +flux_photometry = 10.0 ** (-0.4 * magnitudes) +flux_photometry_errors = flux_photometry * np.log(10.0) / 2.5 * magnitude_errors + +ftr = WDfitter() +ftr.fit( + atmosphere="H", + filters=filters, + photometry=flux_photometry, + photometry_errors=flux_photometry_errors, + photometry_space="flux", + independent=["Teff"], + initial_guess=[13000.0], + logg=7.5, + distance=10.0, + distance_err=0.1, +) +ftr.show_best_fit(display=False) +``` + +### API migration (v0.0.13 -> this branch) + +| Old API name | New canonical name | +| --- | --- | +| `mags` | `photometry` | +| `mag_errors` | `photometry_errors` | +| `fluxes` | `photometry` + `photometry_space="flux"` | +| `flux_errors` | `photometry_errors` + `photometry_space="flux"` | +| `best_fit_mag` / `best_fit_flux` | `best_fit_photometry` | + ### Reddening model The default setup assumes the provided reddening is the total amount at the diff --git a/docs/source/_static/interpolation_rbf_fixed_start.png b/docs/source/_static/interpolation_rbf_fixed_start.png new file mode 100644 index 0000000..a488209 Binary files /dev/null and b/docs/source/_static/interpolation_rbf_fixed_start.png differ diff --git a/docs/source/_static/interpolation_rbf_multistart.png b/docs/source/_static/interpolation_rbf_multistart.png new file mode 100644 index 0000000..f10aa41 Binary files /dev/null and b/docs/source/_static/interpolation_rbf_multistart.png differ diff --git a/docs/source/background/cooling.rst b/docs/source/background/cooling.rst index f28463e..df85434 100644 --- a/docs/source/background/cooling.rst +++ b/docs/source/background/cooling.rst @@ -2,6 +2,6 @@ White Dwarf Cooling =================== -Most of the internal energy of a WD is the residual heat from the progenitor once it passed the planetary nebula phase. However, there are various physical processes that can provide an appreciable amount of energy and govern the cooling rate of a WD at different stages. Following the time sequence in which the physical processes that has direct effects to the photo-luminosity: (1) in the first :math:`10^{8}-10^{9}` years, shell burning of hydrogen via pp-chain can contribute up to 30% of the total luminosity (`Renedo et al. 2010 `_). (2) Neutrino losses -- contribute to a significant fraction of energy loss in the early time of WDs when they were still hot, in the case of massive WDs, neutrino bremsstrahlung effect must also be taken into account (`Martin, Georg \& Achim 1994 `_, `Itoh et al. 1996 `_). (3) Gravitational settling of :math:`^{22}`Ne in intermediate to massive WDs releases sufficient gravitational potential energy to prolong the cooling times (`Deloye \& Bildsten 2002 `_, `Althaus et al. 2010 `_). The heavier :math:`^{22}`Ne relative to the environment that is dominated by carbon, oxygen and nitrogen leads to a slow settling towards the core. This effect is the most obvious in the old and metal-rich systems, such as NGC 6791 (`Bedin et al. 2008 `_, `Garvia-Berro et al. 2010 `_). (4) In the late time of the WD evolution, convection plays a significant role in slowing down the cooling. As temperature decreases, the convective zone grows deeper into the interior and eventually reaches the degenerate core (see Figure 11 from `Althaus, Corsico, Isern, & García-Berro, 2010 `_). This efficiently replenish the energy radiated away from the photosphere, thus this process known as the convective coupling, modifies the relations between the WD luminosity and core temperature (`D'Antona \& Mazzitelli 1989 `_, `Fontaine, Brassard \& Bergeron 2001 `_). (5) Crystallisation occurs as the non-degenerate ions evolve from gas to fluid and eventually solid. The liquid-solid transition releases latent heat that slows down the cooling process. This also couples with the release of gravitational energy associated with changes in the carbon-oxygen profile (`Salaris et al. 1997 `_) when the heavier oxygen-rich crystals displace carbon as a result of gravitational settling. Depending on the changes in the carbon-oxygen abundance profile, and the choice of phase diagram of a carbon-oxygen mixture, it modifies the rate of cooling and this specific effect is colloquially known as the (6) Phase Separation effect. (7) Coulomb Interactions modify the thermodynamical properties of the ionic gas, in particular the specific heat. Its strength is determined by the Coulomb coupling parameters. At first, the parameter is small, it slowly increases as an WD cools and the ions begin to change from gas to liquid and eventually form lattice. This releases latent heat that contribute to ~5% of the total luminosity (`Shaviv \& Kovetz 1976 `_). At late time, few modes of the lattice are excited, the heat capacity drops according to the Debye law, this results in enhanced cooling. This process kicks in after :math:`10^9` yr for a :math:`1.0\,\odot` WD and over a Hubble time for a :math:`0.5\,\odot` WD. See below the plot of the luminosity as a function of cooling age. +Most of the internal energy of a WD is the residual heat from the progenitor once it passed the planetary nebula phase. However, there are various physical processes that can provide an appreciable amount of energy and govern the cooling rate of a WD at different stages. Following the time sequence in which the physical processes that has direct effects to the photo-luminosity: (1) in the first :math:`10^{8}-10^{9}` years, shell burning of hydrogen via pp-chain can contribute up to 30% of the total luminosity (`Renedo et al. 2010 `_). (2) Neutrino losses -- contribute to a significant fraction of energy loss in the early time of WDs when they were still hot, in the case of massive WDs, neutrino bremsstrahlung effect must also be taken into account (`Martin, Georg and Achim 1994 `_, `Itoh et al. 1996 `_). (3) Gravitational settling of :math:`^{22}\mathrm{Ne}` in intermediate to massive WDs releases sufficient gravitational potential energy to prolong the cooling times (`Deloye and Bildsten 2002 `_, `Althaus et al. 2010 `_). The heavier :math:`^{22}\mathrm{Ne}` relative to the environment that is dominated by carbon, oxygen and nitrogen leads to a slow settling towards the core. This effect is the most obvious in the old and metal-rich systems, such as NGC 6791 (`Bedin et al. 2008 `_, `Garvia-Berro et al. 2010 `_). (4) In the late time of the WD evolution, convection plays a significant role in slowing down the cooling. As temperature decreases, the convective zone grows deeper into the interior and eventually reaches the degenerate core (see Figure 11 from `Althaus, Corsico, Isern, and García-Berro, 2010 `_). This efficiently replenish the energy radiated away from the photosphere, thus this process known as the convective coupling, modifies the relations between the WD luminosity and core temperature (`D'Antona and Mazzitelli 1989 `_, `Fontaine, Brassard and Bergeron 2001 `_). (5) Crystallisation occurs as the non-degenerate ions evolve from gas to fluid and eventually solid. The liquid-solid transition releases latent heat that slows down the cooling process. This also couples with the release of gravitational energy associated with changes in the carbon-oxygen profile (`Salaris et al. 1997 `_) when the heavier oxygen-rich crystals displace carbon as a result of gravitational settling. Depending on the changes in the carbon-oxygen abundance profile, and the choice of phase diagram of a carbon-oxygen mixture, it modifies the rate of cooling and this specific effect is colloquially known as the (6) Phase Separation effect. (7) Coulomb Interactions modify the thermodynamical properties of the ionic gas, in particular the specific heat. Its strength is determined by the Coulomb coupling parameters. At first, the parameter is small, it slowly increases as an WD cools and the ions begin to change from gas to liquid and eventually form lattice. This releases latent heat that contribute to ~5% of the total luminosity (`Shaviv and Kovetz 1976 `_). At late time, few modes of the lattice are excited, the heat capacity drops according to the Debye law, this results in enhanced cooling. This process kicks in after :math:`10^9` yr for a :math:`1.0\,\odot` WD and over a Hubble time for a :math:`0.5\,\odot` WD. See below the plot of the luminosity as a function of cooling age. -.. figure:: ../_static/DA_cooling_model_from_plotter.png \ No newline at end of file +.. figure:: ../_static/DA_cooling_model_from_plotter.png diff --git a/docs/source/background/extrapolation.rst b/docs/source/background/extrapolation.rst new file mode 100644 index 0000000..6f58d7b --- /dev/null +++ b/docs/source/background/extrapolation.rst @@ -0,0 +1,55 @@ +================================= +Interpolator Extrapolation Limits +================================= + +`AtmosphereModelReader.interp_am` and +`CoolingModelReader.compute_cooling_age_interpolator` keep +``allow_extrapolation=False`` by default. + +When extrapolation is enabled: + +- queries are limited to at most 50% beyond each axis span, +- non-finite outputs are replaced by smooth fallback interpolation, +- physically impossible outputs are sanitised: + + - atmosphere ``Teff/logg/mass/age`` are forced to finite positive values, + - cooling age is forced to finite positive values, + - cooling rate (``dL/dt``) is forced to finite non-positive values. + +Estimated success rates below were measured with 1000 random samples at each +extrapolation fraction and report the fraction of *extrapolated* outputs that +remained physical and finite. + +Atmosphere Interpolator (DA, ``dependent="Teff"``, ``independent=["logg","Mbol"]``) +-------------------------------------------------------------------------------------- + ++----------------------+---------+----------+ +| Extrapolation Level | CT | RBF | ++======================+=========+==========+ +| 10% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 20% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 30% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 40% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 50% | 1.000 | 1.000 | ++----------------------+---------+----------+ + +Cooling Interpolator (age + cooling-rate physical checks) +---------------------------------------------------------- + ++----------------------+---------+----------+ +| Extrapolation Level | CT | RBF | ++======================+=========+==========+ +| 10% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 20% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 30% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 40% | 1.000 | 1.000 | ++----------------------+---------+----------+ +| 50% | 1.000 | 1.000 | ++----------------------+---------+----------+ diff --git a/docs/source/background/interpolation_scheme.rst b/docs/source/background/interpolation_scheme.rst new file mode 100644 index 0000000..c4b97eb --- /dev/null +++ b/docs/source/background/interpolation_scheme.rst @@ -0,0 +1,70 @@ +================================= +Interpolation Scheme Diagnostics +================================= + +The ``kernel`` argument used with ``extinction_convolved=False`` selects the +radial-basis interpolation scheme for the tabulated reddening profile. The +available ``linear``, ``cubic``, and ``quintic`` kernels are useful numerical +approximations. They do not make pivot-wavelength extinction equivalent to +filter-by-filter extinction convolution. + +Higher-order kernels need particular care. Cubic and quintic kernels introduce +more curvature into the fitted surface than the linear kernel. Combined with a +non-linear photometric fit, this can make a poor starting point converge to a +different local minimum. + +Fixed Starting Point Diagnostic +=============================== + +The diagnostic below fitted 10,000 extinction-stratified GF21 DA candidates +with a fixed initial point of ``(Teff, logg) = (10000 K, 8.0)``. Each panel +compares a convolved-extinction fit with an interpolated-extinction fit using +the named RBF kernel. The linear result follows the reference relation closely. +The cubic result has a broader high-temperature branch and the quintic result +has a pronounced low-temperature branch. Their reported Pearson correlations +with the GF21 temperatures are 0.9416 and 0.5338, respectively, compared with +0.9863 for the linear fit. + +.. figure:: ../_static/interpolation_rbf_fixed_start.png + :width: 100% + :alt: Fixed-start comparison of convolved and RBF-interpolated reddening fits. + + Fixed-start RBF diagnostic. Cubic and quintic interpolated fits can select + erroneous solution branches when fitting starts far from the best solution. + +This diagnostic is a failure-mode check, not a measurement of a kernel-only +effect. High-order RBF profile shape, sparse reddening-table sampling, and +optimizer initialization can all contribute. It does show that a single, +generic starting point is insufficient validation for cubic or quintic fits. + +Multiple Starting Points +======================== + +For each source, use several physically plausible starting points and retain +the finite solution with the lowest :math:`\chi^2`. For example, the diagnostic +used ``(10000 K, 8.0)``, ``(25000 K, 7.0)``, and ``(70000 K, 8.0)``. This is a +deterministic multistart search, not posterior sampling. + +The 1,000-source diagnostic below applies that procedure. The selected +solutions give closely matched temperature relations for the three kernels. +This does not prove that interpolated extinction is equivalent to convolved +extinction; it shows that the obvious incorrect local-minimum solutions were +removed before comparing interpolation schemes. + +.. figure:: ../_static/interpolation_rbf_multistart.png + :width: 100% + :alt: Multistart comparison of convolved and RBF-interpolated reddening fits. + + RBF diagnostic after choosing the lowest finite :math:`\chi^2` solution from + three starting points per source and kernel. + +Recommended Practice +==================== + +* Prefer ``extinction_convolved=True`` for science results. It evaluates + extinction through the model spectrum and passband. +* Treat ``extinction_convolved=False`` as a pivot-wavelength approximation. +* For cubic or quintic RBF kernels, evaluate multiple initial guesses and + inspect both fitted parameters and :math:`\chi^2` before accepting a result. + + diff --git a/docs/source/conf.py b/docs/source/conf.py index 0253470..a405dbb 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -10,8 +10,10 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # +import configparser import os import sys +from pathlib import Path sys.path.insert(0, os.path.abspath("../..")) sys.path.insert(0, os.path.abspath("../../src")) @@ -19,9 +21,18 @@ # -- Project information ----------------------------------------------------- project = "WDPhotTools" -copyright = "2020-2025, Marco C Lam" +copyright = "2020-2026, Marco C Lam" author = "Marco C Lam" -__version__ = "0.0.10" + + +def _load_version(): + parser = configparser.ConfigParser() + setup_cfg = Path(__file__).resolve().parents[2] / "setup.cfg" + parser.read(setup_cfg) + return parser.get("metadata", "version", fallback="0.0.0") + + +__version__ = _load_version() # The full version, including alpha/beta/rc tags version = __version__ @@ -50,7 +61,7 @@ exclude_patterns = [] # The suffix of source filenames. -source_suffix = ".rst" +source_suffix = {".rst": "restructuredtext"} # The encoding of source files. source_encoding = "utf-8-sig" @@ -74,11 +85,8 @@ master_doc = "index" on_rtd = os.environ.get("READTHEDOCS", None) == "True" -if on_rtd: # only import and set the theme if we're building docs locally - import sphinx_rtd_theme - +if on_rtd: # Use RTD theme on Read the Docs builds. html_theme = "sphinx_rtd_theme" - html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] html_static_path = ["_static"] else: html_theme = "alabaster" diff --git a/docs/source/index.rst b/docs/source/index.rst index 0c0e76e..df14503 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -27,6 +27,8 @@ User Guide background/atmosphere background/cooling + background/extrapolation + background/interpolation_scheme background/photometric_fit background/wdlf @@ -108,4 +110,4 @@ Indices and tables * :ref:`genindex` * :ref:`modindex` -* :ref:`search` \ No newline at end of file +* :ref:`search` diff --git a/docs/source/tutorials/fitting.rst b/docs/source/tutorials/fitting.rst index c49a6ba..bbba470 100644 --- a/docs/source/tutorials/fitting.rst +++ b/docs/source/tutorials/fitting.rst @@ -32,19 +32,22 @@ Below is the default setup of the `fit` method, we will go through the 8 cases o self, atmosphere=["H", "He"], filters=["G3", "G3_BP", "G3_RP"], - mags=[None, None, None], - mag_errors=[1.0, 1.0, 1.0], + photometry=None, + photometry_errors=None, + photometry_space="magnitude", allow_none=False, distance=None, distance_err=None, extinction_convolved=True, kernel="cubic", - rv=0.0, + rv=3.1, ebv=0.0, + ra=None, + dec=None, independent=["Mbol", "logg"], initial_guess=[10.0, 8.0], logg=8.0, - atmosphere_interpolator="RBF", + atmosphere_interpolator="CT", reuse_interpolator=False, method="minimize", nwalkers=100, @@ -53,11 +56,13 @@ Below is the default setup of the `fit` method, we will go through the 8 cases o progress=True, refine=False, refine_bounds=[5.0, 95.0], + prior=log_dummy_prior, kwargs_for_RBF={}, kwargs_for_CT={}, kwargs_for_minimize={}, kwargs_for_least_squares={}, kwargs_for_emcee={}, + allow_extrapolation=False, ): ... @@ -66,6 +71,41 @@ Below is the default setup of the `fit` method, we will go through the 8 cases o Below shows the example of how the fitter can be configured to fit differently, the argument that is worth particular mentioning is the `independent`. It refers to the mathemtical term *independent variable*, which are the parameters to be fitted. It accepts a list of 1 or 2 strings, one from `["Mbol", "Teff"]`, and one from `["logg", "mass"]`. When distance is to be fitted, it is not configured here, but instead a `None` should be provided to the `distance` argument. See further down for the examples. When a distance is provided, **its uncertainty has to be provided**. +To fit in relative flux space, set `photometry_space="flux"` and provide `photometry` and `photometry_errors`. + +.. code:: python + + ftr.fit( + atmosphere=["H"], + filters=["G3", "G3_BP", "G3_RP"], + photometry=[f_g, f_bp, f_rp], + photometry_errors=[f_g_err, f_bp_err, f_rp_err], + photometry_space="flux", + distance=21.0, + distance_err=0.1, + independent=["Mbol"], + initial_guess=[10.0], + logg=7.81, + ) + +API migration notes (v0.0.13 -> this branch): + +.. list-table:: + :header-rows: 1 + + * - Old API name + - New canonical name + * - ``mags`` + - ``photometry`` + * - ``mag_errors`` + - ``photometry_errors`` + * - ``fluxes`` + - ``photometry`` with ``photometry_space="flux"`` + * - ``flux_errors`` + - ``photometry_errors`` with ``photometry_space="flux"`` + * - ``best_fit_mag`` / ``best_fit_flux`` + - ``best_fit_photometry`` + **Case 1** When a log(g) is not included in the `independent` list, it will assume a fixed surface gravitiy as provided by `logg`, which is defaulted to 8.0, in this case we want to fit for the bolometric magnitude with a surface gravity of 7.81 for a DA at 21.0 pc with a reddening of `E(B-V) = 0.1` and `Rv` of 3.1 where the extinction is computed by convolving the filter profiles with the DA spectra. The magnitudes and uncertainties in the Gaia eDR3 are some variables `a`, `b`, and `c`. @@ -77,8 +117,8 @@ When a log(g) is not included in the `independent` list, it will assume a fixed ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=21.0, distance_err=0.1, extinction_convolved=True, @@ -97,8 +137,8 @@ Compared to case 1, this fits for the logg, so we needs to add `"logg"` to the ` ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=21.0, distance_err=0.1, extinction_convolved=True, @@ -117,8 +157,8 @@ Compared to case 1, this fits for the distance, but we need to change two things ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=None, extinction_convolved=True, kernel="cubic", @@ -136,8 +176,8 @@ This requires a very simple change, compared to case 1, we change `ebv` to 0.0, ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=21.0, distance_err=0.1, ebv=0.0, @@ -153,8 +193,8 @@ This is a combination of case 3 and 4, and on top, if we opt to use the other in ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=None, ebv=0.0, independent=["Mbol"], @@ -172,8 +212,8 @@ This is a combination of case 2 and 4. We are also demonstrating how to modify t ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=20.1, distance_err=0.1, ebv=0.0, @@ -198,8 +238,8 @@ This is the setup that is the most likely to fail because it is fitting 3 unknow ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=None, rv=3.1, ebv=0.1, @@ -220,8 +260,8 @@ This is the same as case 7 except the reddening is not considered (ebv is set to ftr.fit(atmosphere=["H"], filters=["G3", "G3_BP", "G3_RP"], - mags=[a, b, c], - mag_errors=[a_err, b_err, c_err], + photometry=[a, b, c], + photometry_errors=[a_err, b_err, c_err], distance=None, ebv=0.0, independent=["Mbol", "logg"], diff --git a/example/compare_ps_cooling_rates.py b/example/compare_ps_cooling_rates.py index afea6c9..56ed783 100644 --- a/example/compare_ps_cooling_rates.py +++ b/example/compare_ps_cooling_rates.py @@ -10,6 +10,7 @@ from WDPhotTools import atmosphere_model_reader as amr from WDPhotTools import theoretical_lf +from example_utils import get_example_output_dir try: @@ -17,6 +18,7 @@ except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() # Salaris' model with phase separation wdlf = theoretical_lf.WDLF() @@ -120,4 +122,4 @@ plt.suptitle("log(dL/dt) contour plot") plt.subplots_adjust(wspace=0.0, left=0.075, right=0.975) -plt.savefig(os.path.join(HERE, "example_output", "compare_ps_cooling_rates.png")) +plt.savefig(os.path.join(OUTPUT_DIR, "compare_ps_cooling_rates.png")) diff --git a/example/converting_from_gaia_to_sdss.py b/example/converting_from_gaia_to_sdss.py index e46ad3e..1c7e513a 100644 --- a/example/converting_from_gaia_to_sdss.py +++ b/example/converting_from_gaia_to_sdss.py @@ -10,12 +10,15 @@ from matplotlib import pyplot as plt from WDPhotTools import atmosphere_model_reader as amr +from example_utils import get_example_output_dir try: HERE = os.path.dirname(os.path.realpath(__file__)) except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() + _hdu = fits.open(os.path.join(HERE, "GaiaEDR3_WD_SDSSspec.fits"))[1] data = _hdu.data size = len(data) @@ -168,7 +171,7 @@ plt.subplots_adjust(top=0.95, bottom=0.075, left=0.08, right=0.975, wspace=0.0, hspace=0.1) -plt.savefig(os.path.join(HERE, "example_output", "gaia_to_sdss_cc_diagram.png")) +plt.savefig(os.path.join(OUTPUT_DIR, "gaia_to_sdss_cc_diagram.png")) # Plot the residual of the catalogue value to converted values # Top row is the residual in u from the interpolation of G_BP, G & G_RP @@ -240,7 +243,7 @@ plt.subplots_adjust(top=0.95, bottom=0.07, left=0.08, right=0.975, wspace=0.0, hspace=0.15) -plt.savefig(os.path.join(HERE, "example_output", "gaia_to_sdss_ugr_residual.png")) +plt.savefig(os.path.join(OUTPUT_DIR, "gaia_to_sdss_ugr_residual.png")) # Plot the residual of the catalogue value to derived values from (logg, Teff) @@ -271,4 +274,4 @@ plt.subplots_adjust(top=0.95, bottom=0.075, left=0.075, right=0.975, wspace=0.0, hspace=0.15) -plt.savefig(os.path.join(HERE, "example_output", "gaia_to_sdss_ugr_residual_logg_teff.png")) +plt.savefig(os.path.join(OUTPUT_DIR, "gaia_to_sdss_ugr_residual_logg_teff.png")) diff --git a/example/example_utils.py b/example/example_utils.py new file mode 100644 index 0000000..3ad12bb --- /dev/null +++ b/example/example_utils.py @@ -0,0 +1,23 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""Shared helpers for example scripts.""" + +import os + + +def get_example_output_dir(): + """ + Return the output directory used by example scripts. + + The default location is `/example_output`. Set the environment + variable `WDPHOTTOOLS_EXAMPLE_OUTPUT_DIR` to override this location. + """ + + here = os.path.dirname(os.path.realpath(__file__)) + output_dir = os.environ.get( + "WDPHOTTOOLS_EXAMPLE_OUTPUT_DIR", + os.path.join(here, "example_output"), + ) + os.makedirs(output_dir, exist_ok=True) + return output_dir diff --git a/example/fitting_atmospheric_properties.py b/example/fitting_atmospheric_properties.py index 12381d1..57493df 100644 --- a/example/fitting_atmospheric_properties.py +++ b/example/fitting_atmospheric_properties.py @@ -6,6 +6,7 @@ import os from WDPhotTools.fitter import WDfitter +from example_utils import get_example_output_dir try: @@ -13,6 +14,7 @@ except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() ftr = WDfitter() @@ -31,7 +33,7 @@ "H_mko", "K_mko", ], - mags=[ + photometry=[ 21.1437, 19.9678, 19.4993, @@ -44,7 +46,7 @@ 19.96 - 1.39, 20.40 - 1.85, ], - mag_errors=[ + photometry_errors=[ 0.0321, 0.0229, 0.0083, @@ -68,7 +70,7 @@ ftr.show_best_fit( display=False, savefig=True, - folder=os.path.join(HERE, "example_output"), + folder=OUTPUT_DIR, filename="PSOJ1801p6254", ) @@ -92,7 +94,7 @@ "H_mko", "K_mko", ], - mags=[ + photometry=[ 21.1437, 19.9678, 19.4993, @@ -105,7 +107,7 @@ 19.96 - 1.39, 20.40 - 1.85, ], - mag_errors=[ + photometry_errors=[ 0.0321, 0.0229, 0.0083, @@ -129,7 +131,7 @@ ftr.show_best_fit( display=False, savefig=True, - folder=os.path.join(HERE, "example_output"), + folder=OUTPUT_DIR, filename="PSOJ1801p6254_least_squares", ) @@ -151,7 +153,7 @@ "H_mko", "K_mko", ], - mags=[ + photometry=[ 21.1437, 19.9678, 19.4993, @@ -164,7 +166,7 @@ 19.96 - 1.39, 20.40 - 1.85, ], - mag_errors=[ + photometry_errors=[ 0.0321, 0.0229, 0.0083, @@ -192,14 +194,14 @@ atmosphere="H", display=False, savefig=True, - folder=os.path.join(HERE, "example_output"), + folder=OUTPUT_DIR, filename="PSOJ1801p6254_emcee", ) ftr.show_corner_plot( figsize=(10, 10), display=True, savefig=True, - folder=os.path.join(HERE, "example_output"), + folder=OUTPUT_DIR, filename="PSOJ1801p6254_emcee_corner", kwarg={ "quantiles": [0.158655, 0.5, 0.841345], diff --git a/example/generate_all_plots.py b/example/generate_all_plots.py index f19bc83..d172bcb 100644 --- a/example/generate_all_plots.py +++ b/example/generate_all_plots.py @@ -6,7 +6,9 @@ import numpy as np from WDPhotTools import theoretical_lf +from example_utils import get_example_output_dir +OUTPUT_DIR = get_example_output_dir() wdlf = theoretical_lf.WDLF() @@ -20,7 +22,7 @@ cooling_model_use_mag=False, imf_log=True, display=True, - folder="example_output", + folder=OUTPUT_DIR, ext=["png", "pdf"], savefig=True, ) diff --git a/example/generate_burst_and_exponential_sfr_wdlf.py b/example/generate_burst_and_exponential_sfr_wdlf.py index b7a59dc..896150b 100644 --- a/example/generate_burst_and_exponential_sfr_wdlf.py +++ b/example/generate_burst_and_exponential_sfr_wdlf.py @@ -9,6 +9,7 @@ from matplotlib import pyplot as plt from WDPhotTools import theoretical_lf +from example_utils import get_example_output_dir try: @@ -16,6 +17,7 @@ except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() wdlf = theoretical_lf.WDLF() @@ -40,7 +42,7 @@ mag=mag, passband="G3", save_csv=True, - folder=os.path.join(HERE, "example_output"), + folder=OUTPUT_DIR, ) ax1.plot(mag, np.log10(burst_density), label=f"{age / 1e9:.2f} Gyr") @@ -52,7 +54,7 @@ mag=mag, passband="G3", save_csv=True, - folder=os.path.join(HERE, "example_output"), + folder=OUTPUT_DIR, ) ax2.plot(mag, np.log10(decay_density), label=f"{age / 1e9:.2f} Gyr") @@ -65,10 +67,8 @@ ax1.set_title("Star Formation History: 1 Gyr Burst") fig1.savefig( os.path.join( - HERE, - "example_output", - "burst_C16_C08_montreal_co_da_20_" - "montreal_co_da_20_montreal_co_da_20.png", + OUTPUT_DIR, + "burst_C16_C08_montreal_co_da_20_" "montreal_co_da_20_montreal_co_da_20.png", ) ) @@ -81,9 +81,7 @@ ax2.set_title("Star Formation History: Exponential Decay") fig2.savefig( os.path.join( - HERE, - "example_output", - "decay_C16_C08_montreal_co_da_20_" - "montreal_co_da_20_montreal_co_da_20.png", + OUTPUT_DIR, + "decay_C16_C08_montreal_co_da_20_" "montreal_co_da_20_montreal_co_da_20.png", ) ) diff --git a/example/generate_constant_sfr_wdlf.py b/example/generate_constant_sfr_wdlf.py index da29070..dfb8100 100644 --- a/example/generate_constant_sfr_wdlf.py +++ b/example/generate_constant_sfr_wdlf.py @@ -9,6 +9,7 @@ from matplotlib import pyplot as plt from WDPhotTools import theoretical_lf +from example_utils import get_example_output_dir try: @@ -16,6 +17,7 @@ except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() wdlf = theoretical_lf.WDLF() wdlf.compute_cooling_age_interpolator() @@ -29,9 +31,7 @@ for i, age in enumerate(age_list): # Constant SFR wdlf.set_sfr_model(age=age) - _, constant_density = wdlf.compute_density( - mag=mag, save_csv=True, folder=os.path.join(HERE, "example_output") - ) + _, constant_density = wdlf.compute_density(mag=mag, save_csv=True, folder=OUTPUT_DIR) ax1.plot(mag, np.log10(constant_density), label=f"{age / 1e9:.2f} Gyr") ax1.legend() @@ -43,9 +43,7 @@ ax1.set_title("Star Formation History: Constant") fig1.savefig( os.path.join( - HERE, - "example_output", - "constant_C16_C08_montreal_co_da_20_" - "montreal_co_da_20_montreal_co_da_20.png", + OUTPUT_DIR, + "constant_C16_C08_montreal_co_da_20_" "montreal_co_da_20_montreal_co_da_20.png", ) ) diff --git a/example/generate_cooling_model_with_plotter.py b/example/generate_cooling_model_with_plotter.py index b07b364..26f627a 100644 --- a/example/generate_cooling_model_with_plotter.py +++ b/example/generate_cooling_model_with_plotter.py @@ -4,10 +4,13 @@ """Plot the default cooling model""" from WDPhotTools import plotter +from example_utils import get_example_output_dir + +OUTPUT_DIR = get_example_output_dir() plotter.plot_cooling_model( mass=[0.2, 0.4, 0.6, 0.8, 1.0], savefig=True, - folder="example_output", + folder=OUTPUT_DIR, filename="DA_cooling_model_from_plotter", ) diff --git a/example/generate_cooling_tracks.py b/example/generate_cooling_tracks.py index 20e69ca..66c8b56 100644 --- a/example/generate_cooling_tracks.py +++ b/example/generate_cooling_tracks.py @@ -9,6 +9,7 @@ import numpy as np from WDPhotTools.atmosphere_model_reader import AtmosphereModelReader +from example_utils import get_example_output_dir try: @@ -16,6 +17,7 @@ except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() atm = AtmosphereModelReader() @@ -43,4 +45,4 @@ plt.ylabel("G / mag") plt.title("DA Cooling tracks") plt.tight_layout() -plt.savefig(os.path.join(HERE, "example_output", "DA_cooling_tracks.png")) +plt.savefig(os.path.join(OUTPUT_DIR, "DA_cooling_tracks.png")) diff --git a/example/generate_cooling_tracks_with_plotter.py b/example/generate_cooling_tracks_with_plotter.py index 252d0d2..0d97b10 100644 --- a/example/generate_cooling_tracks_with_plotter.py +++ b/example/generate_cooling_tracks_with_plotter.py @@ -1,22 +1,16 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -"""Create the default cooling tracks""" - -import os +"""Create the default cooling tracks.""" from WDPhotTools import plotter +from example_utils import get_example_output_dir - -try: - HERE = os.path.dirname(os.path.realpath(__file__)) - -except NameError: - HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() plotter.plot_atmosphere_model( invert_yaxis=True, savefig=True, - folder="example_output", + folder=OUTPUT_DIR, filename="DA_cooling_tracks_from_plotter", ) diff --git a/example/generate_wdlf_compare_sfr.py b/example/generate_wdlf_compare_sfr.py index 750d61b..98ee249 100644 --- a/example/generate_wdlf_compare_sfr.py +++ b/example/generate_wdlf_compare_sfr.py @@ -9,12 +9,14 @@ import numpy as np from WDPhotTools import theoretical_lf +from example_utils import get_example_output_dir try: HERE = os.path.dirname(os.path.realpath(__file__)) except NameError: HERE = os.path.dirname(os.path.realpath(__name__)) +OUTPUT_DIR = get_example_output_dir() wdlf = theoretical_lf.WDLF() wdlf.set_ifmr_model("C08") wdlf.compute_cooling_age_interpolator() @@ -22,9 +24,7 @@ mag = np.arange(0, 20.0, 0.1) age_list = 1e9 * np.arange(2, 15, 2) -fig1, (ax1, ax2, ax3) = plt.subplots( - 3, 1, sharex=True, sharey=True, figsize=(10, 15) -) +fig1, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, sharey=True, figsize=(10, 15)) for i, age in enumerate(age_list): # Constant SFR @@ -58,8 +58,7 @@ plt.savefig( os.path.join( - HERE, - "example_output", + OUTPUT_DIR, "wdlf_compare_sfr.png", ) ) diff --git a/src/WDPhotTools/atmosphere_model_reader.py b/src/WDPhotTools/atmosphere_model_reader.py index 3b40094..00fd6fd 100644 --- a/src/WDPhotTools/atmosphere_model_reader.py +++ b/src/WDPhotTools/atmosphere_model_reader.py @@ -315,6 +315,9 @@ def interp_am( Keyword argument for the interpolator. See `scipy.interpolate.RBFInterpolator`. kwargs_for_CT: dict (Default: {'fill_value': -np.inf, 'tol': 1e-10, 'maxiter': 100000}) Keyword argument for the interpolator. See `scipy.interpolate.CloughTocher2DInterpolator`. + allow_extrapolation: bool (Default: False) + If True, allow smooth extrapolation up to 50% outside each interpolation axis span. Extrapolated values + are sanitised to avoid non-finite and physically impossible outputs for positive-definite parameters. Returns ------- @@ -338,6 +341,8 @@ def interp_am( } _kwargs_for_CT.update(**kwargs_for_CT) + max_extrapolation_fraction = 0.5 + # DA atmosphere if atmosphere.lower() in ["h", "hydrogen", "da"]: model = self.model_da @@ -352,6 +357,38 @@ def interp_am( "provided {}.format(atmosphere.lower())" ) + _dependent_values = np.asarray(model[dependent], dtype=float).reshape(-1) + _finite_dependent = _dependent_values[np.isfinite(_dependent_values)] + _finite_positive_dependent = _finite_dependent[_finite_dependent > 0.0] + _dependent_floor = ( + float(np.nanmin(_finite_positive_dependent)) + if _finite_positive_dependent.size + else float(np.finfo(float).tiny) + ) + _dependent_fallback = float(np.nanmedian(_finite_dependent)) if _finite_dependent.size else 0.0 + _positive_dependent = {"Teff", "logg", "mass", "age"} + + def _clip_with_fraction(values, value_min, value_max, fraction): + span = max(float(value_max - value_min), float(np.finfo(float).eps)) + return np.clip(values, value_min - span * fraction, value_max + span * fraction) + + def _sanitize_output(values): + out = np.asarray(values, dtype=float).reshape(-1) + if dependent in _positive_dependent: + out[~np.isfinite(out)] = _dependent_floor + out = np.maximum(out, _dependent_floor) + else: + out[~np.isfinite(out)] = _dependent_fallback + return out + + def _sanitize_extrapolated(values, mask_oob): + out = np.asarray(values, dtype=float).reshape(-1) + mask_oob = np.asarray(mask_oob, dtype=bool).reshape(-1) + sanitize_mask = mask_oob | (~np.isfinite(out)) + if np.any(sanitize_mask): + out[sanitize_mask] = _sanitize_output(out[sanitize_mask]) + return out + independent = np.asarray(independent, dtype=object).reshape(-1) independent_list = self.column_key @@ -373,14 +410,19 @@ def interp_am( independent = np.array([independent_list[_independent_arg_0], independent_list[_independent_arg_1]]) - arg_0 = model[independent[0]] - arg_1 = model[independent[1]] + arg_0 = np.asarray(model[independent[0]], dtype=float).reshape(-1) + arg_1_raw = np.asarray(model[independent[1]], dtype=float).reshape(-1) - arg_1_min = np.nanmin(arg_1) - arg_1_max = np.nanmax(arg_1) + arg_1_min = float(np.nanmin(arg_1_raw)) + arg_1_max = float(np.nanmax(arg_1_raw)) + arg_1_floor = ( + float(np.nanmin(arg_1_raw[arg_1_raw > 0.0])) if np.any(arg_1_raw > 0.0) else float(np.finfo(float).tiny) + ) if independent[1] in ["Teff", "age"]: - arg_1 = np.log10(arg_1) + arg_1 = np.log10(np.maximum(arg_1_raw, arg_1_floor)) + else: + arg_1 = arg_1_raw if interpolator.lower() == "ct": # Interpolate with the scipy CloughTocher2DInterpolator @@ -389,12 +431,37 @@ def interp_am( model[dependent], **_kwargs_for_CT, ) + _rbf_extrapolator = None + if allow_extrapolation: + _rbf_extrapolator = RBFInterpolator( + np.stack((arg_0, arg_1), -1), + model[dependent], + **_kwargs_for_RBF, + ) def atmosphere_interpolator(_x): + _x_raw = np.asarray(_x).reshape(-1).astype(float) + mask_oob = (_x_raw < arg_1_min) | (_x_raw > arg_1_max) + if allow_extrapolation: + _x_eval = _clip_with_fraction(_x_raw, arg_1_min, arg_1_max, max_extrapolation_fraction) + else: + _x_eval = np.clip(_x_raw, arg_1_min, arg_1_max) + if independent[1] in ["Teff", "age"]: - _x = np.log10(_x) + _x_eval = np.log10(np.maximum(_x_eval, arg_1_floor)) - return _atmosphere_interpolator(logg, _x) + _logg = np.full(_x_eval.size, logg, dtype=float) + out = np.asarray(_atmosphere_interpolator(_logg, _x_eval)).reshape(-1) + + if allow_extrapolation: + bad = ~np.isfinite(out) + if np.any(bad): + out[bad] = np.asarray( + _rbf_extrapolator(np.column_stack((_logg[bad], _x_eval[bad]))) + ).reshape(-1) + out = _sanitize_extrapolated(out, mask_oob) + + return out elif interpolator.lower() == "rbf": # Interpolate with the scipy RBFInterpolator @@ -405,18 +472,23 @@ def atmosphere_interpolator(_x): ) def atmosphere_interpolator(_x): - _x_arr = np.asarray(_x).reshape(-1).astype(float) - length = _x_arr.size + _x_raw = np.asarray(_x).reshape(-1).astype(float) + length = _x_raw.size _logg = np.full(length, logg, dtype=float) + mask_oob = (_x_raw < arg_1_min) | (_x_raw > arg_1_max) if not allow_extrapolation: - _x_arr[_x_arr < arg_1_min] = arg_1_min - _x_arr[_x_arr > arg_1_max] = arg_1_max + _x_eval = np.clip(_x_raw, arg_1_min, arg_1_max) + else: + _x_eval = _clip_with_fraction(_x_raw, arg_1_min, arg_1_max, max_extrapolation_fraction) if independent[1] in ["Teff", "age"]: - _x_arr = np.log10(_x_arr) + _x_eval = np.log10(np.maximum(_x_eval, arg_1_floor)) - return _atmosphere_interpolator(np.column_stack((_logg, _x_arr))) + out = np.asarray(_atmosphere_interpolator(np.column_stack((_logg, _x_eval)))).reshape(-1) + if allow_extrapolation: + out = _sanitize_extrapolated(out, mask_oob) + return out else: raise ValueError("Interpolator should be CT or RBF, {interpolator} is given.") @@ -429,19 +501,86 @@ def atmosphere_interpolator(_x): independent = np.array([independent_list[_independent_arg_0], independent_list[_independent_arg_1]]) - arg_0 = model[independent[0]] - arg_1 = model[independent[1]] + arg_0_raw = np.asarray(model[independent[0]], dtype=float).reshape(-1) + arg_1_raw = np.asarray(model[independent[1]], dtype=float).reshape(-1) - arg_0_min = np.nanmin(arg_0) - arg_0_max = np.nanmax(arg_0) - arg_1_min = np.nanmin(arg_1) - arg_1_max = np.nanmax(arg_1) + arg_0_min = float(np.nanmin(arg_0_raw)) + arg_0_max = float(np.nanmax(arg_0_raw)) + arg_1_min = float(np.nanmin(arg_1_raw)) + arg_1_max = float(np.nanmax(arg_1_raw)) + arg_0_floor = ( + float(np.nanmin(arg_0_raw[arg_0_raw > 0.0])) if np.any(arg_0_raw > 0.0) else float(np.finfo(float).tiny) + ) + arg_1_floor = ( + float(np.nanmin(arg_1_raw[arg_1_raw > 0.0])) if np.any(arg_1_raw > 0.0) else float(np.finfo(float).tiny) + ) if independent[0] in ["Teff", "age"]: - arg_0 = np.log10(arg_0) + arg_0 = np.log10(np.maximum(arg_0_raw, arg_0_floor)) + else: + arg_0 = arg_0_raw if independent[1] in ["Teff", "age"]: - arg_1 = np.log10(arg_1) + arg_1 = np.log10(np.maximum(arg_1_raw, arg_1_floor)) + else: + arg_1 = arg_1_raw + + def _prepare_two_input_arrays(x0, x1=None): + if x1 is None: + arr = np.asarray(x0).reshape(-1) + if arr.size >= 2: + x_0, x_1 = arr[0], arr[1] + elif arr.size == 1: + x_0, x_1 = arr[0], arr[0] + else: + x_0, x_1 = -np.inf, -np.inf + else: + x_0, x_1 = x0, x1 + + if isinstance(x_0, (float, int, np.integer)): + length0 = 1 + else: + length0 = np.asarray(x_0).size + + if isinstance(x_1, (float, int, np.integer)): + length1 = 1 + else: + length1 = np.asarray(x_1).size + + if length0 == length1: + pass + elif (length0 == 1) and (length1 > 1): + x_0 = [x_0] * length1 + length0 = length1 + elif (length0 > 1) and (length1 == 1): + x_1 = [x_1] * length0 + length1 = length0 + else: + raise ValueError( + "Either one variable is a float, int or of size 1, or two variables should have the same " + "size." + ) + + return ( + np.asarray(x_0).reshape(-1).astype(float), + np.asarray(x_1).reshape(-1).astype(float), + ) + + def _transform_coordinates(x_0, x_1): + mask_oob = (x_0 < arg_0_min) | (x_0 > arg_0_max) | (x_1 < arg_1_min) | (x_1 > arg_1_max) + if allow_extrapolation: + _x_0 = _clip_with_fraction(x_0, arg_0_min, arg_0_max, max_extrapolation_fraction) + _x_1 = _clip_with_fraction(x_1, arg_1_min, arg_1_max, max_extrapolation_fraction) + else: + _x_0 = np.clip(x_0, arg_0_min, arg_0_max) + _x_1 = np.clip(x_1, arg_1_min, arg_1_max) + + if independent[0] in ["Teff", "age"]: + _x_0 = np.log10(np.maximum(_x_0, arg_0_floor)) + if independent[1] in ["Teff", "age"]: + _x_1 = np.log10(np.maximum(_x_1, arg_1_floor)) + + return _x_0, _x_1, mask_oob if interpolator.lower() == "ct": # Interpolate with the scipy CloughTocher2DInterpolator @@ -450,60 +589,30 @@ def atmosphere_interpolator(_x): model[dependent], **_kwargs_for_CT, ) + _rbf_extrapolator = None + if allow_extrapolation: + _rbf_extrapolator = RBFInterpolator( + np.stack((arg_0, arg_1), -1), + model[dependent], + **_kwargs_for_RBF, + ) def atmosphere_interpolator(x0, x1=None): - # Support scalar/array inputs for both coordinates, with simple broadcasting - if x1 is None: - arr = np.asarray(x0).reshape(-1) - if arr.size >= 2: - x_0, x_1 = arr[0], arr[1] - elif arr.size == 1: - x_0, x_1 = arr[0], arr[0] - else: - x_0, x_1 = -np.inf, -np.inf - else: - x_0, x_1 = x0, x1 - - if isinstance(x_0, (float, int, np.integer)): - length0 = 1 - else: - length0 = np.asarray(x_0).size - - if isinstance(x_1, (float, int, np.integer)): - length1 = 1 - else: - length1 = np.asarray(x_1).size - - if length0 == length1: - pass - elif (length0 == 1) and (length1 > 1): - x_0 = [x_0] * length1 - length0 = length1 - elif (length0 > 1) and (length1 == 1): - x_1 = [x_1] * length0 - length1 = length0 - else: - raise ValueError( - "Either one variable is a float, int or of size 1, or two variables should have the same" - "size." - ) - - _x_0 = np.asarray(x_0).reshape(-1).astype(float) - _x_1 = np.asarray(x_1).reshape(-1).astype(float) - - # mark out-of-range inputs to avoid excessive extrapolation - mask_oob = (_x_0 < arg_0_min) | (_x_0 > arg_0_max) | (_x_1 < arg_1_min) | (_x_1 > arg_1_max) - - if independent[0] in ["Teff", "age"]: - _x_0 = np.log10(_x_0) - - if independent[1] in ["Teff", "age"]: - _x_1 = np.log10(_x_1) - + x_0, x_1 = _prepare_two_input_arrays(x0, x1) + _x_0, _x_1, mask_oob = _transform_coordinates(x_0, x_1) out = _atmosphere_interpolator(_x_0, _x_1) out = np.asarray(out).reshape(-1) - if not allow_extrapolation and np.any(mask_oob): + + if allow_extrapolation: + bad = ~np.isfinite(out) + if np.any(bad): + out[bad] = np.asarray(_rbf_extrapolator(np.column_stack((_x_0[bad], _x_1[bad])))).reshape( + -1 + ) + out = _sanitize_extrapolated(out, mask_oob) + elif np.any(mask_oob): out[mask_oob] = -np.inf + return out elif interpolator.lower() == "rbf": @@ -514,63 +623,14 @@ def atmosphere_interpolator(x0, x1=None): **_kwargs_for_RBF, ) - def atmosphere_interpolator(*x): - # Accept (x0, x1) or single array-like; use first two values, duplicate if only one - if len(x) == 2: - x_0, x_1 = x - elif len(x) == 1: - arr = np.asarray(x[0]).reshape(-1) - if arr.size >= 2: - x_0, x_1 = arr[0], arr[1] - elif arr.size == 1: - x_0, x_1 = arr[0], arr[0] - else: - x_0, x_1 = -np.inf, -np.inf - else: - x_0, x_1 = -np.inf, -np.inf - - if isinstance(x_0, (float, int, np.integer)): - length0 = 1 - else: - length0 = np.asarray(x_0).size - - if isinstance(x_1, (float, int, np.integer)): - length1 = 1 - else: - length1 = np.asarray(x_1).size - - if length0 == length1: - pass - - elif (length0 == 1) & (length1 > 1): - x_0 = [x_0] * length1 - length0 = length1 - - elif (length0 > 1) & (length1 == 1): - x_1 = [x_1] * length0 - length1 = length0 - - else: - raise ValueError( - "Either one variable is a float, int or of size 1, or two variables should have the same " - "size." - ) - - _x_0 = np.asarray(x_0).reshape(-1).astype(float) - _x_1 = np.asarray(x_1).reshape(-1).astype(float) - - # mark out-of-range inputs to avoid excessive extrapolation - mask_oob = (_x_0 < arg_0_min) | (_x_0 > arg_0_max) | (_x_1 < arg_1_min) | (_x_1 > arg_1_max) - - if independent[0] in ["Teff", "age"]: - _x_0 = np.log10(_x_0) - - if independent[1] in ["Teff", "age"]: - _x_1 = np.log10(_x_1) - + def atmosphere_interpolator(x0, x1=None): + x_0, x_1 = _prepare_two_input_arrays(x0, x1) + _x_0, _x_1, mask_oob = _transform_coordinates(x_0, x_1) out = _atmosphere_interpolator(np.column_stack((_x_0, _x_1))) out = np.asarray(out).reshape(-1) - if np.any(mask_oob): + if allow_extrapolation: + out = _sanitize_extrapolated(out, mask_oob) + elif np.any(mask_oob): out[mask_oob] = -np.inf return out diff --git a/src/WDPhotTools/cooling_model_reader.py b/src/WDPhotTools/cooling_model_reader.py index 01582b0..499c336 100644 --- a/src/WDPhotTools/cooling_model_reader.py +++ b/src/WDPhotTools/cooling_model_reader.py @@ -1723,6 +1723,10 @@ def compute_cooling_age_interpolator( Keyword argument for the interpolator. See `scipy.interpolate.RBFInterpolator`. kwargs_for_CT: dict (Default: {'fill_value': -np.inf, 'tol': 1e-10, 'maxiter': 100000}) Keyword argument for the interpolator. See `scipy.interpolate.CloughTocher2DInterpolator`. + allow_extrapolation: bool (Default: False) + If True, allow smooth extrapolation up to 50% outside each interpolation axis span. Extrapolated cooling + ages are sanitised to finite positive values and extrapolated cooling rates are sanitised to finite + non-positive values. """ @@ -1844,98 +1848,175 @@ def compute_cooling_age_interpolator( } _kwargs_for_RBF.update(**kwargs_for_RBF) + max_extrapolation_fraction = 0.5 + lum_grid = np.log10(self.luminosity) + lum_min = float(np.nanmin(lum_grid)) + lum_max = float(np.nanmax(lum_grid)) + mass_min = float(np.nanmin(self.mass)) + mass_max = float(np.nanmax(self.mass)) + age_floor = ( + float(np.nanmin(self.age[self.age > 0.0])) if np.any(self.age > 0.0) else float(np.finfo(float).tiny) + ) + finite_age = self.age[np.isfinite(self.age)] + age_fallback = float(np.nanmedian(finite_age)) if finite_age.size else age_floor + + def _clip_with_fraction(values, value_min, value_max, fraction): + span = max(float(value_max - value_min), float(np.finfo(float).eps)) + return np.clip(values, value_min - span * fraction, value_max + span * fraction) + + def _prepare_inputs(x_0, x_1): + _x_0 = np.asarray(x_0, dtype=float).reshape(-1) + _x_1 = np.asarray(x_1, dtype=float).reshape(-1) + + if (_x_0.size == 1) and (_x_1.size > 1): + _x_0 = np.repeat(_x_0, _x_1.size) + elif (_x_1.size == 1) and (_x_0.size > 1): + _x_1 = np.repeat(_x_1, _x_0.size) + elif _x_0.size != _x_1.size: + raise ValueError("x_0 and x_1 should have the same size or one should have size 1.") + + mask_oob = (_x_0 < lum_min) | (_x_0 > lum_max) | (_x_1 < mass_min) | (_x_1 > mass_max) + + if allow_extrapolation: + _x_0 = _clip_with_fraction(_x_0, lum_min, lum_max, max_extrapolation_fraction) + _x_1 = _clip_with_fraction(_x_1, mass_min, mass_max, max_extrapolation_fraction) + else: + _x_0 = np.clip(_x_0, lum_min, lum_max) + _x_1 = np.clip(_x_1, mass_min, mass_max) + + return _x_0, _x_1, mask_oob + + def _sanitize_cooling_age(values): + out = np.asarray(values, dtype=float).reshape(-1) + out[~np.isfinite(out)] = age_fallback + out = np.maximum(out, age_floor) + return out + + def _sanitize_extrapolated(values, mask_oob, sanitizer): + out = np.asarray(values, dtype=float).reshape(-1) + mask_oob = np.asarray(mask_oob, dtype=bool).reshape(-1) + sanitize_mask = mask_oob | (~np.isfinite(out)) + if np.any(sanitize_mask): + out[sanitize_mask] = sanitizer(out[sanitize_mask]) + return out + if interpolator.lower() == "ct": - # Interpolate with the scipy CloughTocher2DInterpolator - self.cooling_interpolator = CloughTocher2DInterpolator( - (np.log10(self.luminosity), self.mass), + _cooling_interpolator = CloughTocher2DInterpolator( + (lum_grid, self.mass), self.age, **_kwargs_for_CT, ) + _cooling_rbf_extrapolator = None + if allow_extrapolation: + _cooling_rbf_extrapolator = RBFInterpolator( + np.stack((lum_grid, self.mass), -1), + self.age, + **_kwargs_for_RBF, + ) + + def cooling_interpolator(x_0, x_1): + _x_0, _x_1, mask_oob = _prepare_inputs(x_0, x_1) + out = np.asarray(_cooling_interpolator(_x_0, _x_1), dtype=float).reshape(-1) + if allow_extrapolation: + bad = ~np.isfinite(out) + if np.any(bad): + out[bad] = np.asarray( + _cooling_rbf_extrapolator(np.column_stack((_x_0[bad], _x_1[bad]))), + dtype=float, + ).reshape(-1) + out = _sanitize_extrapolated(out, mask_oob, _sanitize_cooling_age) + elif np.any(mask_oob): + out[mask_oob] = -np.inf + return out + + self.cooling_interpolator = cooling_interpolator elif interpolator.lower() == "rbf": - # Interpolate with the scipy RBFInterpolator _cooling_interpolator = RBFInterpolator( - np.stack((np.log10(self.luminosity), self.mass), -1), + np.stack((lum_grid, self.mass), -1), self.age, **_kwargs_for_RBF, ) - lum_min = np.nanmin(np.log10(self.luminosity)) - lum_max = np.nanmax(np.log10(self.luminosity)) - mass_min = np.nanmin(self.mass) - mass_max = np.nanmax(self.mass) - def cooling_interpolator(x_0, x_1): - _x_0 = np.array(x_0) - _x_1 = np.array(x_1) - - if (_x_0.size == 1) & (_x_1.size > 1): - _x_0 = np.repeat(_x_0, _x_1.size) - - if (_x_1.size == 1) & (_x_0.size > 1): - _x_1 = np.repeat(_x_1, _x_0.size) - - if not allow_extrapolation: - _x_0[_x_0 < lum_min] = lum_min - _x_0[_x_0 > lum_max] = lum_max - _x_1[_x_1 < mass_min] = mass_min - _x_1[_x_1 > mass_max] = mass_max - - length0 = _x_0.size - - return _cooling_interpolator(np.array([_x_0, _x_1], dtype="object").T.reshape(length0, 2)) + _x_0, _x_1, mask_oob = _prepare_inputs(x_0, x_1) + out = np.asarray( + _cooling_interpolator(np.column_stack((_x_0, _x_1))), + dtype=float, + ).reshape(-1) + if allow_extrapolation: + out = _sanitize_extrapolated(out, mask_oob, _sanitize_cooling_age) + elif np.any(mask_oob): + out[mask_oob] = -np.inf + return out self.cooling_interpolator = cooling_interpolator else: raise ValueError(f"Interpolator should be CT or RBF, {interpolator} is given.") - self.dLdt = self._itp2d_gradient(self.cooling_interpolator, np.log10(self.luminosity), self.mass) - + self.dLdt = self._itp2d_gradient(self.cooling_interpolator, lum_grid, self.mass) finite_mask = np.isfinite(self.dLdt) + finite_rate = self.dLdt[finite_mask] + finite_rate_nonpositive = finite_rate[finite_rate <= 0.0] + rate_fallback = float(np.nanmedian(finite_rate_nonpositive)) if finite_rate_nonpositive.size else 0.0 + + def _sanitize_cooling_rate(values): + out = np.asarray(values, dtype=float).reshape(-1) + out[~np.isfinite(out)] = rate_fallback + out[out > 0.0] = 0.0 + return out + if interpolator.lower() == "ct": - self.cooling_rate_interpolator = CloughTocher2DInterpolator( - ( - np.log10(self.luminosity)[finite_mask], - self.mass[finite_mask], - ), + _cooling_rate_interpolator = CloughTocher2DInterpolator( + (lum_grid[finite_mask], self.mass[finite_mask]), self.dLdt[finite_mask], **_kwargs_for_CT, ) + _cooling_rate_rbf_extrapolator = None + if allow_extrapolation: + _cooling_rate_rbf_extrapolator = RBFInterpolator( + np.stack((lum_grid[finite_mask], self.mass[finite_mask]), -1), + self.dLdt[finite_mask], + **_kwargs_for_RBF, + ) + + def cooling_rate_interpolator(x_0, x_1): + _x_0, _x_1, mask_oob = _prepare_inputs(x_0, x_1) + out = np.asarray(_cooling_rate_interpolator(_x_0, _x_1), dtype=float).reshape(-1) + if allow_extrapolation: + bad = ~np.isfinite(out) + if np.any(bad): + out[bad] = np.asarray( + _cooling_rate_rbf_extrapolator(np.column_stack((_x_0[bad], _x_1[bad]))), + dtype=float, + ).reshape(-1) + out = _sanitize_extrapolated(out, mask_oob, _sanitize_cooling_rate) + elif np.any(mask_oob): + out[mask_oob] = -np.inf + return out + + self.cooling_rate_interpolator = cooling_rate_interpolator elif interpolator.lower() == "rbf": - # Interpolate with the scipy RBFInterpolator _cooling_rate_interpolator = RBFInterpolator( - np.stack((np.log10(self.luminosity)[finite_mask], self.mass[finite_mask]), -1), + np.stack((lum_grid[finite_mask], self.mass[finite_mask]), -1), self.dLdt[finite_mask], **_kwargs_for_RBF, ) - lum_min = np.nanmin(np.log10(self.luminosity)) - lum_max = np.nanmax(np.log10(self.luminosity)) - mass_min = np.nanmin(self.mass) - mass_max = np.nanmax(self.mass) - def cooling_rate_interpolator(x_0, x_1): - _x_0 = np.asarray(x_0) - _x_1 = np.asarray(x_1) - - if (_x_0.size == 1) & (_x_1.size > 1): - _x_0 = np.repeat(_x_0, _x_1.size) - - if (_x_1.size == 1) & (_x_0.size > 1): - _x_0 = np.repeat(_x_1, _x_0.size) - - if not allow_extrapolation: - _x_0[_x_0 < lum_min] = lum_min - _x_0[_x_0 > lum_max] = lum_max - _x_1[_x_1 < mass_min] = mass_min - _x_1[_x_1 > mass_max] = mass_max - - length0 = _x_0.size - - return _cooling_rate_interpolator(np.asarray([_x_0, _x_1], dtype="object").T.reshape(length0, 2)) + _x_0, _x_1, mask_oob = _prepare_inputs(x_0, x_1) + out = np.asarray( + _cooling_rate_interpolator(np.column_stack((_x_0, _x_1))), + dtype=float, + ).reshape(-1) + if allow_extrapolation: + out = _sanitize_extrapolated(out, mask_oob, _sanitize_cooling_rate) + elif np.any(mask_oob): + out[mask_oob] = -np.inf + return out self.cooling_rate_interpolator = cooling_rate_interpolator diff --git a/src/WDPhotTools/diff2_functions_least_square.py b/src/WDPhotTools/diff2_functions_least_square.py index 387a9f3..b8b22b1 100644 --- a/src/WDPhotTools/diff2_functions_least_square.py +++ b/src/WDPhotTools/diff2_functions_least_square.py @@ -2,6 +2,40 @@ from .extinction import get_extinction_fraction +_MAG_DISTANCE_FACTOR = 2.17147241 +_MAG_TO_FRAC_FLUX_VAR = 0.8483036976765438 + + +def _compute_residual_terms( + obs, + errors, + model_mag, + distance, + distance_err, + photometry_space, +): + """Compute chi-square terms in either magnitude or relative flux space.""" + if photometry_space == "magnitude": + if distance_err is None: + e2 = errors**2.0 + else: + # 5 / ln(10) converts fractional distance error to magnitude error. + # (ln(10) / 2.5)^2 converts magnitude variance to fractional flux variance. + e2 = (errors**2.0 + (distance_err / distance * _MAG_DISTANCE_FACTOR) ** 2.0) * _MAG_TO_FRAC_FLUX_VAR + d2 = ((10.0 ** ((obs - model_mag) / 2.5) - 1.0) ** 2.0) / e2 + return d2, e2 + + if photometry_space == "flux": + model_flux = 10.0 ** (-0.4 * model_mag) + e2 = errors**2.0 + if distance_err is not None: + # Relative flux follows d^-2, so sigma_f = 2 * sigma_d / d * f. + e2 = e2 + (2.0 * distance_err / distance * model_flux) ** 2.0 + d2 = ((obs - model_flux) ** 2.0) / e2 + return d2, e2 + + raise ValueError("Unknown photometry_space. Please choose from 'magnitude' and 'flux'.") + def diff2( _x, @@ -11,6 +45,7 @@ def diff2( distance_err, interpolator_filter, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for scipy.optimize.least_squares). @@ -22,14 +57,15 @@ def diff2( for interp in interpolator_filter: mag.append(interp(_x[:2])) - mag = np.asarray(mag).reshape(-1) - - # 5 / ln(10) = 2.17147241 converts fractional distance error to magnitude error - # (ln(10) / 2.5)^2 = 0.8483036976765438 converts magnitude variance to fractional flux variance - e2 = (errors**2.0 + (distance_err / distance * 2.17147241) ** 2.0) * 0.8483036976765438 - - d2 = ((10.0 ** ((obs - mag - 5.0 * np.log10(distance) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 - + model_mag = np.asarray(mag).reshape(-1) + 5.0 * np.log10(distance) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=distance, + distance_err=distance_err, + photometry_space=photometry_space, + ) # Ensure finite residuals d2 = np.where(np.isfinite(d2), d2, np.float64(1e30)) if return_err: @@ -39,7 +75,14 @@ def diff2( return d2 -def diff2_distance(_x, obs, errors, interpolator_filter, return_err): +def diff2_distance( + _x, + obs, + errors, + interpolator_filter, + return_err, + photometry_space="magnitude", +): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for scipy.optimize.least_squares). @@ -58,9 +101,15 @@ def diff2_distance(_x, obs, errors, interpolator_filter, return_err): for interp in interpolator_filter: mag.append(interp(_x[:-1])) - mag = np.asarray(mag).reshape(-1) - e2 = errors**2.0 - d2 = ((10.0 ** ((obs - mag - 5.0 * np.log10(_x[-1]) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + 5.0 * np.log10(_x[-1]) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=_x[-1], + distance_err=None, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -91,6 +140,7 @@ def diff2_distance_red_interpolated( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided. @@ -116,9 +166,15 @@ def diff2_distance_red_interpolated( extinction_fraction = get_extinction_fraction(_x[-1], ra, dec, zmin, zmax) av = np.array([i(rv) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - e2 = errors**2.0 - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(_x[-1]) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(_x[-1]) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=_x[-1], + distance_err=None, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -149,6 +205,7 @@ def diff2_distance_red_interpolated_fixed_logg( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided. @@ -174,10 +231,15 @@ def diff2_distance_red_interpolated_fixed_logg( extinction_fraction = get_extinction_fraction(_x[-1], ra, dec, zmin, zmax) av = np.array([i(rv) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - e2 = errors**2.0 - - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(_x[-1]) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(_x[-1]) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=_x[-1], + distance_err=None, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -210,6 +272,7 @@ def diff2_distance_red_filter( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when @@ -238,9 +301,15 @@ def diff2_distance_red_filter( teff = float(np.asarray(interpolator_teff(_x[:2])).reshape(-1)[0]) logg = _x[logg_pos] av = np.array([i([logg, teff, rv]) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - e2 = errors**2.0 - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(_x[-1]) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(_x[-1]) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=_x[-1], + distance_err=None, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -273,6 +342,7 @@ def diff2_distance_red_filter_fixed_logg( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided. @@ -299,10 +369,15 @@ def diff2_distance_red_filter_fixed_logg( teff = float(np.asarray(interpolator_teff(_x[:-1])).reshape(-1)[0]) av = np.array([i([logg, teff, rv]) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - e2 = errors**2.0 - - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(_x[-1]) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(_x[-1]) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=_x[-1], + distance_err=None, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -335,6 +410,7 @@ def diff2_red_interpolated( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value. @@ -353,13 +429,15 @@ def diff2_red_interpolated( extinction_fraction = get_extinction_fraction(distance, ra, dec, zmin, zmax) av = np.array([i(rv) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - - # 5 / ln(10) = 2.17147241 converts fractional distance error to mag error - # (ln(10) / 2.5)^2 = 0.8483036976765438 converts magnitude variance to fractional flux variance - # of residuals squared - e2 = (errors**2.0 + (distance_err / distance * 2.17147241) ** 2.0) * 0.8483036976765438 - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(distance) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(distance) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=distance, + distance_err=distance_err, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -394,6 +472,7 @@ def diff2_red_filter( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for scipy.optimize.least_square). @@ -422,13 +501,15 @@ def diff2_red_filter( logg = _x[logg_pos] av = np.array([i([logg, teff, rv]) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - - # 5 / ln(10) = 2.17147241 converts fractional distance error to mag error - # (ln(10) / 2.5)^2 = 0.8483036976765438 converts magnitude variance to fractional flux variance - # of residuals squared - e2 = (errors**2.0 + (distance_err / distance * 2.17147241) ** 2.0) * 0.8483036976765438 - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(distance) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(distance) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=distance, + distance_err=distance_err, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: @@ -463,6 +544,7 @@ def diff2_red_filter_fixed_logg( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for scipy.optimize.least_square). @@ -482,13 +564,15 @@ def diff2_red_filter_fixed_logg( teff = float(np.asarray(interpolator_teff(_x)).reshape(-1)[0]) av = np.array([i([logg, teff, rv]) for i in reddening_vector]).reshape(-1) * ebv * extinction_fraction - mag = np.asarray(mag).reshape(-1) - - # 5 / ln(10) = 2.17147241 converts fractional distance error to mag error - # (ln(10) / 2.5)^2 = 0.8483036976765438 converts magnitude variance to fractional flux variance - # of residuals squared - e2 = (errors**2.0 + (distance_err / distance * 2.17147241) ** 2.0) * 0.8483036976765438 - d2 = ((10.0 ** ((obs - av - mag - 5.0 * np.log10(distance) + 5.0) / 2.5) - 1.0) ** 2.0) / e2 + model_mag = np.asarray(mag).reshape(-1) + av + 5.0 * np.log10(distance) - 5.0 + d2, e2 = _compute_residual_terms( + obs=obs, + errors=errors, + model_mag=model_mag, + distance=distance, + distance_err=distance_err, + photometry_space=photometry_space, + ) if np.isfinite(d2).all(): if return_err: diff --git a/src/WDPhotTools/diff2_functions_minimize.py b/src/WDPhotTools/diff2_functions_minimize.py index 090ed62..d8e9b5c 100644 --- a/src/WDPhotTools/diff2_functions_minimize.py +++ b/src/WDPhotTools/diff2_functions_minimize.py @@ -13,13 +13,31 @@ ) -def diff2_summed(_x, obs, errors, distance, distance_err, interpolator_filter, return_err): +def diff2_summed( + _x, + obs, + errors, + distance, + distance_err, + interpolator_filter, + return_err, + photometry_space="magnitude", +): """ Internal method for computing the ch2-squared value (for scipy.optimize.minimize). """ - d2, e2 = diff2(_x, obs, errors, distance, distance_err, interpolator_filter, True) + d2, e2 = diff2( + _x, + obs, + errors, + distance, + distance_err, + interpolator_filter, + True, + photometry_space=photometry_space, + ) if return_err: return np.sum(d2), 1.0 / np.sum(1.0 / e2) @@ -46,6 +64,7 @@ def diff2_red_filter_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for scipy.optimize.minimize). @@ -69,6 +88,7 @@ def diff2_red_filter_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: @@ -96,6 +116,7 @@ def diff2_red_filter_fixed_logg_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for scipy.optimize.minimize). @@ -119,6 +140,7 @@ def diff2_red_filter_fixed_logg_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: @@ -144,6 +166,7 @@ def diff2_red_interpolated_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for @@ -167,6 +190,7 @@ def diff2_red_interpolated_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: return np.sum(d2), 1.0 / np.sum(1.0 / e2) @@ -191,6 +215,7 @@ def diff2_distance_red_filter_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for @@ -214,6 +239,7 @@ def diff2_distance_red_filter_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: return np.sum(d2), 1.0 / np.sum(1.0 / e2) @@ -238,6 +264,7 @@ def diff2_distance_red_filter_fixed_logg_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for @@ -261,6 +288,7 @@ def diff2_distance_red_filter_fixed_logg_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: return np.sum(d2), 1.0 / np.sum(1.0 / e2) @@ -269,14 +297,28 @@ def diff2_distance_red_filter_fixed_logg_summed( return np.sum(d2) -def diff2_distance_summed(_x, obs, errors, interpolator_filter, return_error): +def diff2_distance_summed( + _x, + obs, + errors, + interpolator_filter, + return_error, + photometry_space="magnitude", +): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for scipy.optimize.minimize). """ - d2, e2 = diff2_distance(_x, obs, errors, interpolator_filter, True) + d2, e2 = diff2_distance( + _x, + obs, + errors, + interpolator_filter, + True, + photometry_space=photometry_space, + ) if return_error: return np.sum(d2), 1.0 / np.sum(1.0 / e2) @@ -299,6 +341,7 @@ def diff2_distance_red_interpolated_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for @@ -320,6 +363,7 @@ def diff2_distance_red_interpolated_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: @@ -343,6 +387,7 @@ def diff2_distance_red_interpolated_fixed_logg_summed( zmin, zmax, return_err, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for @@ -364,6 +409,7 @@ def diff2_distance_red_interpolated_fixed_logg_summed( zmin, zmax, True, + photometry_space=photometry_space, ) if return_err: diff --git a/src/WDPhotTools/fitter.py b/src/WDPhotTools/fitter.py index bd6ed76..2a3763e 100644 --- a/src/WDPhotTools/fitter.py +++ b/src/WDPhotTools/fitter.py @@ -4,6 +4,7 @@ """Core of the WD photometry fitter""" import copy +import logging import os import time from functools import partial @@ -55,6 +56,14 @@ plt.rc("font", size=18) plt.rc("legend", fontsize=12) +logger = logging.getLogger(__name__) +if not logger.handlers: + _handler = logging.StreamHandler() + _handler.setFormatter(logging.Formatter("%(asctime)s %(name)s %(levelname)s: %(message)s")) + logger.addHandler(_handler) +logger.setLevel(logging.INFO) +logger.propagate = False + class WDfitter(AtmosphereModelReader): """ @@ -70,7 +79,7 @@ def __init__(self): # Only used if minimize or least_squares are the fitting method self.results = {"H": {}, "He": {}} self.best_fit_params = {"H": {}, "He": {}} - self.best_fit_mag = {"H": [], "He": []} + self.best_fit_photometry = {"H": [], "He": []} # Only used if emcee is the fitting method self.sampler = {"H": [], "He": []} self.samples = {"H": [], "He": []} @@ -195,8 +204,9 @@ def fit( self, atmosphere=["H", "He"], filters=["G3", "G3_BP", "G3_RP"], - mags=[None, None, None], - mag_errors=[1.0, 1.0, 1.0], + photometry=None, + photometry_errors=None, + photometry_space="magnitude", allow_none=False, distance=None, distance_err=None, @@ -227,12 +237,13 @@ def fit( allow_extrapolation=False, ): """ - The method to execute a photometric fit. Pure hydrogen and helium atmospheres fitting are supported. See - `atmosphere_model_reader` for more information. Set allow_none to True so that `mags` can be provided in None - to Default non-detection, it is not used in the fit but it allows the fitter to be reused over a large dataset - where non-detections occur occasionally. In practice, one can add the full list of filters and set None for all - the non-detections, however this is highly inefficent in memory usage: most of the interpolated grid is not - used, and masking takes time. + The method to execute a photometric fit. Pure hydrogen and helium + atmospheres fitting are supported. See `atmosphere_model_reader` for + more information. Set allow_none to True so that `photometry` can be + provided with None values for non-detections. In practice, one can add + the full list of filters and set None for all the non-detections, + however this is highly inefficient in memory usage: most of the + interpolated grid is not used, and masking takes time. Parameters ---------- @@ -240,15 +251,17 @@ def fit( Choose to fit with pure hydrogen atmosphere model and/or pure helium atmosphere model. filters: list/array of str (Default: ['G3', 'G3_BP', 'G3_RP']) Choose the filters to be fitted with. - mags: list/array of float (Default: [None, None, None]) - The magnitudes in the chosen filters, in their respective magnitude system. None can be provided as - non-detection, it does not contribute to the fitting. - mag_errors: list/array of float (Default: [1., 1., 1.]) - The uncertainties in the magnitudes provided. + photometry: list/array of float (Default: None) + Observed photometry in the selected `photometry_space`. + This is required. + photometry_errors: list/array of float (Default: None) + Uncertainties of `photometry`. This is required. + photometry_space: str (Default: 'magnitude') + Choose from 'magnitude' and 'flux'. allow_none: bool (Default: False) - Set to True to detect None in the `mags` list to create a mask, this check requires extra run-time. + Set to True to detect None in the input photometry list to create a mask, this check requires extra run-time. distance: float (Default: None) - The distance to the source, in parsec. Set to None if the distance is to be fitted simultanenous. Provide + The distance to the source, in parsec. Set to None if the distance is to be fitted simultaneously. Provide an initial guess in the `initial_guess`, or it will be initialised at 10.0 pc. distance_err: float (Default: None) The uncertainty of the distance. @@ -274,7 +287,7 @@ def fit( be initialise as 50.0 pc if not provided. logg: float (Default: 8.0) Only used if 'logg' is not included in the `independent` argument. - atmosphere_interpolator: str (Default: 'RBF') + atmosphere_interpolator: str (Default: 'CT') Choose between 'RBF' and 'CT'. reuse_interpolator: bool (Default: False) Set to use the existing interpolated grid, it should be set to True if the same collection of data is @@ -290,10 +303,10 @@ def fit( Number of steps is discarded as burn-in (emcee method only). progress: bool (Default: True) Show the progress of the emcee sampling (emcee method only). - refine: cool (Default: True) + refine: bool (Default: False) Set to True to refine the minimum with `scipy.optimize.minimize`. - refine_bounds: str (Default: [5, 95]) - The bounds of the minimizer are definited by the percentiles ofthe samples. + refine_bounds: list/array of float (Default: [5, 95]) + The bounds of the minimizer are defined by the percentiles of the samples. prior: function (Default: log_dummy_prior) The prior function for the emcee method. Default is a dummy prior that always return 0. kwargs_for_RBF: dict (Default: {}) @@ -303,7 +316,7 @@ def fit( kwargs_for_minimize: dict (Default: {'method': 'Powell', 'options': {'xtol': 0.001}}) Keyword argument for the minimizer, see `scipy.optimize.minimize`. kwargs_for_least_squares: dict (Default: {}) - keywprd argument for the minimizer, see `scipy.optimize.least_squares`. + Keyword argument for the minimizer, see `scipy.optimize.least_squares`. kwargs_for_emcee: dict (Default: {}) Keyword argument for the emcee walker. @@ -347,25 +360,58 @@ def fit( if isinstance(initial_guess, np.ndarray): initial_guess = list(initial_guess.reshape(-1)) + if photometry_space not in ["magnitude", "flux"]: + raise ValueError("Unknown photometry_space. Please choose from 'magnitude' and 'flux'.") + if (photometry is None) or (photometry_errors is None): + raise ValueError("Please provide photometry and photometry_errors for the selected photometry_space.") + + logger.info( + "Starting fit: method=%s space=%s atmospheres=%s filters=%d allow_none=%s", + method, + photometry_space, + atmosphere, + len(filters), + allow_none, + ) if isinstance(distance, (float, int, np.floating)): if not isinstance(distance_err, (float, int, np.floating)): distance_err = np.sqrt(distance) + logger.info( + "distance_err not provided; using sqrt(distance)=%.6g", + distance_err, + ) - if distance is None: - if len(initial_guess) == len(independent): - initial_guess = initial_guess + [50.0] + if distance is None and len(initial_guess) == len(independent): + initial_guess = initial_guess + [50.0] + logger.info("Distance is free parameter; extending initial_guess with default distance=50.0 pc.") - # Mask the data and interpolator if set to detect None - if allow_none: - mask = np.array([m is not None for m in mags], dtype=bool) - mags = np.asarray(mags, dtype=float)[mask] - mag_errors = np.asarray(mag_errors, dtype=float)[mask] - filters = np.asarray(filters, dtype=object)[mask] + filters = np.asarray(filters, dtype=object) + if len(photometry) != len(photometry_errors): + raise ValueError("photometry and photometry_errors must have the same length.") + if len(photometry) != len(filters): + raise ValueError("filters and photometry must have the same length.") + # Mask the data and interpolator if set to detect None. + if allow_none: + mask = np.asarray([m is not None for m in photometry], dtype=bool) + photometry = np.asarray(photometry, dtype=float)[mask] + photometry_errors = np.asarray(photometry_errors, dtype=float)[mask] + filters = filters[mask] + logger.info( + "Applied allow_none mask; using %d/%d photometry points.", + photometry.size, + mask.size, + ) + if photometry.size == 0: + raise ValueError("No valid photometry remains after allow_none masking.") else: - mags = np.array(mags, dtype=float) - mag_errors = np.array(mag_errors, dtype=float) - filters = np.array(filters) + photometry = np.asarray(photometry, dtype=float) + photometry_errors = np.asarray(photometry_errors, dtype=float) + + if photometry.size != photometry_errors.size or photometry.size != filters.size: + raise ValueError( + "Length mismatch after preprocessing: filters, photometry, and photometry_errors must match." + ) if ( ((rv >= 0.0) and (self.reddening_vector is None)) @@ -395,10 +441,18 @@ def fit( diffs[~np.isfinite(diffs)] = np.inf nearest_idx = int(np.argmin(diffs)) ig[idx] = float(grid_vals[nearest_idx]) - print( - "Because extrapolation is not allowed in the initial guess(es) are outside the grid, initial_guess " - f"is updated from {initial_guess} to {ig}." - ) + if not np.allclose( + np.asarray(ig, dtype=float), + np.asarray(initial_guess, dtype=float), + rtol=0.0, + atol=0.0, + equal_nan=True, + ): + logger.warning( + "Adjusted initial_guess to nearest grid point because allow_extrapolation=False: %s -> %s", + initial_guess, + ig, + ) initial_guess = ig # Reuse the interpolator if instructed or possible @@ -408,9 +462,15 @@ def fit( & (self.interpolator[atmosphere[0]] != []) & (len(self.interpolator[atmosphere[0]]) == (len(filters) + 4)) ): - pass + logger.info("Reusing existing interpolators for %d filters.", len(filters)) else: + logger.info( + "Building interpolators: interpolator=%s filters=%d atmospheres=%s", + atmosphere_interpolator, + len(filters), + atmosphere, + ) self.interpolator = {"H": {}, "He": {}} for j in atmosphere: @@ -432,8 +492,9 @@ def fit( self.fitting_params = { "atmosphere": atmosphere, "filters": filters, - "mags": mags, - "mag_errors": mag_errors, + "photometry_space": photometry_space, + "photometry": photometry, + "photometry_errors": photometry_errors, "distance": distance, "distance_err": distance_err, "independent": independent, @@ -464,10 +525,57 @@ def fit( logg_pos_arr = np.where(np.array(self.fitting_params["independent"]) == "logg")[0] logg_pos = int(logg_pos_arr[0]) if logg_pos_arr.size > 0 else None + def _bind_photometry_space(functions): + return {name: partial(func, photometry_space=photometry_space) for name, func in functions.items()} + + objective_minimize = _bind_photometry_space( + { + "main": diff2_summed, + "red_filter": diff2_red_filter_summed, + "red_filter_fixed_logg": diff2_red_filter_fixed_logg_summed, + "red_interpolated": diff2_red_interpolated_summed, + "distance": diff2_distance_summed, + "distance_red_filter": diff2_distance_red_filter_summed, + "distance_red_filter_fixed_logg": diff2_distance_red_filter_fixed_logg_summed, + "distance_red_interpolated": diff2_distance_red_interpolated_summed, + "distance_red_interpolated_fixed_logg": diff2_distance_red_interpolated_fixed_logg_summed, + } + ) + + objective_least_squares = _bind_photometry_space( + { + "main": diff2, + "red_filter": diff2_red_filter, + "red_filter_fixed_logg": diff2_red_filter_fixed_logg, + "red_interpolated": diff2_red_interpolated, + "distance": diff2_distance, + "distance_red_filter": diff2_distance_red_filter, + "distance_red_filter_fixed_logg": diff2_distance_red_filter_fixed_logg, + "distance_red_interpolated": diff2_distance_red_interpolated, + "distance_red_interpolated_fixed_logg": diff2_distance_red_interpolated_fixed_logg, + } + ) + + objective_emcee = _bind_photometry_space( + { + "main": log_likelihood, + "red_filter": log_likelihood_red_filter, + "red_filter_fixed_logg": log_likelihood_red_filter_fixed_logg, + "red_interpolated": log_likelihood_red_interpolated, + "distance": log_likelihood_distance, + "distance_red_filter": log_likelihood_distance_red_filter, + "distance_red_filter_fixed_logg": log_likelihood_distance_red_filter_fixed_logg, + "distance_red_interpolated": log_likelihood_distance_red_interpolated, + "distance_red_interpolated_fixed_logg": log_likelihood_distance_red_interpolated_fixed_logg, + } + ) + # If using the scipy.optimize.minimize() if method == "minimize": + logger.info("Running scipy.optimize.minimize for atmospheres=%s", atmosphere) # Iterative through the list of atmospheres for j in atmosphere: + logger.info("Starting minimize fit for atmosphere=%s", j) if extinction_convolved: interpolator_teff = self.interpolator[j]["Teff"] @@ -477,11 +585,11 @@ def fit( if ebv <= 0.0: # with or without logg takes the same for, it is handled in the interpolator self.results[j] = optimize.minimize( - diff2_distance_summed, + objective_minimize["distance"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], False, ), @@ -492,11 +600,11 @@ def fit( if not extinction_convolved: if "logg" in independent: self.results[j] = optimize.minimize( - diff2_distance_red_interpolated_summed, + objective_minimize["distance_red_interpolated"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], rv, self.extinction_mode, @@ -513,11 +621,11 @@ def fit( else: self.results[j] = optimize.minimize( - diff2_distance_red_interpolated_fixed_logg_summed, + objective_minimize["distance_red_interpolated_fixed_logg"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], rv, self.extinction_mode, @@ -534,11 +642,11 @@ def fit( else: if "logg" in independent: self.results[j] = optimize.minimize( - diff2_distance_red_filter_summed, + objective_minimize["distance_red_filter"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], self.interpolator[j]["Teff"], logg_pos, @@ -557,11 +665,11 @@ def fit( else: self.results[j] = optimize.minimize( - diff2_distance_red_filter_fixed_logg_summed, + objective_minimize["distance_red_filter_fixed_logg"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], self.interpolator[j]["Teff"], logg, @@ -582,11 +690,11 @@ def fit( else: if ebv <= 0.0: self.results[j] = optimize.minimize( - diff2_summed, + objective_minimize["main"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -598,11 +706,11 @@ def fit( else: if not extinction_convolved: self.results[j] = optimize.minimize( - diff2_red_interpolated_summed, + objective_minimize["red_interpolated"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -624,11 +732,11 @@ def fit( _arr = np.where(np.array(self.fitting_params["independent"]) == "logg")[0] logg_pos = int(_arr[0]) if _arr.size > 0 else None self.results[j] = optimize.minimize( - diff2_red_filter_summed, + objective_minimize["red_filter"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -649,11 +757,11 @@ def fit( else: self.results[j] = optimize.minimize( - diff2_red_filter_fixed_logg_summed, + objective_minimize["red_filter_fixed_logg"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -708,8 +816,10 @@ def fit( # If using scipy.optimize.least_squares elif method == "least_squares": + logger.info("Running scipy.optimize.least_squares for atmospheres=%s", atmosphere) # Iterative through the list of atmospheres for j in atmosphere: + logger.info("Starting least_squares fit for atmosphere=%s", j) if extinction_convolved: interpolator_teff = self.interpolator[j]["Teff"] @@ -731,11 +841,11 @@ def fit( else (lb[0], ub[0]) ) self.results[j] = optimize.least_squares( - diff2_distance, + objective_least_squares["distance"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], False, ), @@ -747,11 +857,11 @@ def fit( if not extinction_convolved: if "logg" in independent: self.results[j] = optimize.least_squares( - diff2_distance_red_interpolated, + objective_least_squares["distance_red_interpolated"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], rv, self.extinction_mode, @@ -768,11 +878,11 @@ def fit( else: self.results[j] = optimize.least_squares( - diff2_distance_red_interpolated_fixed_logg, + objective_least_squares["distance_red_interpolated_fixed_logg"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], rv, self.extinction_mode, @@ -789,11 +899,11 @@ def fit( else: if "logg" in independent: self.results[j] = optimize.least_squares( - diff2_distance_red_filter, + objective_least_squares["distance_red_filter"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], self.interpolator[j]["Teff"], logg_pos, @@ -812,11 +922,11 @@ def fit( else: self.results[j] = optimize.least_squares( - diff2_distance_red_filter_fixed_logg, + objective_least_squares["distance_red_filter_fixed_logg"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], self.interpolator[j]["Teff"], logg, @@ -844,11 +954,11 @@ def fit( ub = lohi[:, 1] - eps bounds = (lb.tolist(), ub.tolist()) self.results[j] = optimize.least_squares( - diff2, + objective_least_squares["main"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -865,11 +975,11 @@ def fit( # as the extinction from SFD12 table 6 has no dependency on # temperature and logg self.results[j] = optimize.least_squares( - diff2_red_interpolated, + objective_least_squares["red_interpolated"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -888,11 +998,11 @@ def fit( else: if "logg" in independent: self.results[j] = optimize.least_squares( - diff2_red_filter, + objective_least_squares["red_filter"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -913,11 +1023,11 @@ def fit( else: self.results[j] = optimize.least_squares( - diff2_red_filter_fixed_logg, + objective_least_squares["red_filter_fixed_logg"], initial_guess, args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -975,6 +1085,7 @@ def fit( # If using emcee elif method == "emcee": + logger.info("Running emcee for atmospheres=%s", atmosphere) _initial_guess = np.array(initial_guess) ndim = len(_initial_guess) nwalkers = int(nwalkers) @@ -982,6 +1093,7 @@ def fit( # Iterative through the list of atmospheres for j in atmosphere: + logger.info("Starting emcee fit for atmosphere=%s", j) if extinction_convolved: interpolator_teff = self.interpolator[j]["Teff"] @@ -992,10 +1104,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_distance, + objective_emcee["distance"], args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], prior, ), @@ -1008,10 +1120,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_distance_red_interpolated, + objective_emcee["distance_red_interpolated"], args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], rv, self.extinction_mode, @@ -1030,10 +1142,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_distance_red_interpolated_fixed_logg, + objective_emcee["distance_red_interpolated_fixed_logg"], args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], rv, self.extinction_mode, @@ -1052,10 +1164,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_distance_red_filter, + objective_emcee["distance_red_filter"], args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], interpolator_teff, logg_pos, @@ -1076,10 +1188,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_distance_red_filter_fixed_logg, + objective_emcee["distance_red_filter_fixed_logg"], args=( - mags, - mag_errors, + photometry, + photometry_errors, [self.interpolator[j][i] for i in filters], interpolator_teff, logg, @@ -1102,10 +1214,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood, + objective_emcee["main"], args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -1119,10 +1231,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_red_interpolated, + objective_emcee["red_interpolated"], args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -1144,10 +1256,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_red_filter, + objective_emcee["red_filter"], args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -1173,10 +1285,10 @@ def fit( self.sampler[j] = emcee.EnsembleSampler( nwalkers, ndim, - log_likelihood_red_filter_fixed_logg, + objective_emcee["red_filter_fixed_logg"], args=( - mags, - mag_errors, + photometry, + photometry_errors, distance, distance_err, [self.interpolator[j][i] for i in filters], @@ -1197,6 +1309,14 @@ def fit( self.sampler[j].run_mcmc(pos, nsteps, progress=progress) self.samples[j] = self.sampler[j].get_chain(discard=nburns, flat=True) + logger.info( + "Completed emcee sampling for atmosphere=%s: nwalkers=%d nsteps=%d nburns=%d nsamples=%d", + j, + nwalkers, + nsteps, + nburns, + self.samples[j].shape[0], + ) # Save the best fit results if len(independent) == 1: @@ -1220,7 +1340,11 @@ def fit( kwargs = copy.deepcopy(_kwargs_for_minimize) kwargs["bounds"] = np.percentile(self.samples[j], refine_bounds, axis=0).T - print("Refining") + logger.info( + "Refining emcee solution with minimize for atmosphere=%s bounds=%s", + j, + refine_bounds, + ) _initial_guess = np.percentile(self.samples[j], 50.0, axis=0) @@ -1230,8 +1354,9 @@ def fit( # intial_guess when distance has to be found self.fit( filters=filters, - mags=mags, - mag_errors=mag_errors, + photometry=photometry, + photometry_errors=photometry_errors, + photometry_space=photometry_space, allow_none=allow_none, atmosphere=atmosphere, logg=logg, @@ -1251,8 +1376,9 @@ def fit( else: self.fit( filters=filters, - mags=mags, - mag_errors=mag_errors, + photometry=photometry, + photometry_errors=photometry_errors, + photometry_space=photometry_space, allow_none=allow_none, atmosphere=atmosphere, logg=logg, @@ -1296,16 +1422,13 @@ def fit( else: raise ValueError("Unknown method. Please choose from minimize, least_squares and emcee.") - # Save the pivot wavelength and magnitude for each filter + # Save the pivot wavelength for each filter. self.pivot_wavelengths = [] for i in self.fitting_params["filters"]: self.pivot_wavelengths.append(self.column_wavelengths[i]) for j in atmosphere: - self.best_fit_mag[j] = [] - - for i in self.fitting_params["filters"]: - self.best_fit_mag[j].append(self.best_fit_params[j][i]) + model_mag = [self.best_fit_params[j][i] for i in self.fitting_params["filters"]] for name in ["Teff", "mass", "Mbol", "age"]: if len(independent) == 1: @@ -1375,6 +1498,25 @@ def fit( for i, _f in enumerate(self.fitting_params["filters"]): self.best_fit_params[j]["Av_" + _f] = Av[i] + reddening = np.array( + [self.best_fit_params[j]["Av_" + f] for f in self.fitting_params["filters"]], + dtype=np.float64, + ) + apparent_mag = np.array(model_mag, dtype=np.float64) + self.best_fit_params[j]["dist_mod"] + reddening + if photometry_space == "flux": + self.best_fit_photometry[j] = (10.0 ** (-0.4 * apparent_mag)).tolist() + else: + self.best_fit_photometry[j] = apparent_mag.tolist() + + logger.info( + "Completed fit atmosphere=%s: Teff=%.6g logg=%.6g distance=%.6g chi2=%.6g", + j, + self.best_fit_params[j].get("Teff", np.nan), + self.best_fit_params[j].get("logg", np.nan), + self.best_fit_params[j].get("distance", np.nan), + self.best_fit_params[j].get("chi2", np.nan), + ) + def show_corner_plot( self, figsize=(8, 8), @@ -1502,7 +1644,7 @@ def show_best_fit( return_fig=True, ): """ - Generate a figure with the given and fitted photometry. + Generate a figure with the given and fitted photometry in magnitude or flux space. Parameters ---------- @@ -1543,12 +1685,16 @@ def show_best_fit( fig = plt.figure(figsize=figsize) _ax = fig.gca() + photometry_space = self.fitting_params.get("photometry_space", "magnitude") + + # Plot the photometry provided. + observed = self.fitting_params["photometry"] + observed_err = self.fitting_params["photometry_errors"] - # Plot the photometry provided _ax.errorbar( self.pivot_wavelengths, - self.fitting_params["mags"], - yerr=self.fitting_params["mag_errors"], + observed, + yerr=observed_err, linestyle="None", capsize=3, fmt="s", @@ -1561,11 +1707,10 @@ def show_best_fit( if self.best_fit_params[k] == {}: continue - reddening = [self.best_fit_params[k]["Av_" + f] for f in self.fitting_params["filters"]] - + best_fit = np.array(self.best_fit_photometry[k], dtype=np.float64) _ax.scatter( np.array(self.pivot_wavelengths), - np.array(self.best_fit_mag[k]) + self.best_fit_params[k]["dist_mod"] + np.array(reddening), + best_fit, label=f"Best-fit {k}", color=color[j], zorder=15, @@ -1573,11 +1718,14 @@ def show_best_fit( # Other decorative stuff _ax.legend() - _ax.invert_yaxis() _ax.grid() _ax.set_xlabel("Wavelength / A") - _ax.set_ylabel("Magnitude / mag") + if photometry_space == "flux": + _ax.set_ylabel("Relative flux") + else: + _ax.invert_yaxis() + _ax.set_ylabel("Magnitude / mag") # Configure the title if title is None: diff --git a/src/WDPhotTools/likelihood_functions.py b/src/WDPhotTools/likelihood_functions.py index 66a53b9..0b14b4c 100644 --- a/src/WDPhotTools/likelihood_functions.py +++ b/src/WDPhotTools/likelihood_functions.py @@ -30,13 +30,23 @@ def log_likelihood( distance_err, interpolator_filter, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). """ - d2, e2 = diff2(_x, obs, errors, distance, distance_err, interpolator_filter, True) + d2, e2 = diff2( + _x, + obs, + errors, + distance, + distance_err, + interpolator_filter, + True, + photometry_space=photometry_space, + ) p = prior(*_x) if np.isfinite(d2).all(): @@ -46,14 +56,28 @@ def log_likelihood( return -np.inf -def log_likelihood_distance(_x, obs, errors, interpolator_filter, prior): +def log_likelihood_distance( + _x, + obs, + errors, + interpolator_filter, + prior, + photometry_space="magnitude", +): """ Internal method for computing the ch2-squared value in cases when the distance is not provided (for emcee). """ - d2, e2 = diff2_distance(_x, obs, errors, interpolator_filter, True) + d2, e2 = diff2_distance( + _x, + obs, + errors, + interpolator_filter, + True, + photometry_space=photometry_space, + ) p = prior(*_x) if np.isfinite(d2).all(): @@ -79,6 +103,7 @@ def log_likelihood_distance_red_filter( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -101,6 +126,7 @@ def log_likelihood_distance_red_filter( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) @@ -127,6 +153,7 @@ def log_likelihood_distance_red_filter_fixed_logg( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -149,6 +176,7 @@ def log_likelihood_distance_red_filter_fixed_logg( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) @@ -173,6 +201,7 @@ def log_likelihood_distance_red_interpolated( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -193,6 +222,7 @@ def log_likelihood_distance_red_interpolated( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) @@ -217,6 +247,7 @@ def log_likelihood_distance_red_interpolated_fixed_logg( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -237,6 +268,7 @@ def log_likelihood_distance_red_interpolated_fixed_logg( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) @@ -265,6 +297,7 @@ def log_likelihood_red_filter( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -289,6 +322,7 @@ def log_likelihood_red_filter( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) @@ -317,6 +351,7 @@ def log_likelihood_red_filter_fixed_logg( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -341,6 +376,7 @@ def log_likelihood_red_filter_fixed_logg( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) @@ -367,6 +403,7 @@ def log_likelihood_red_interpolated( z_min, z_max, prior, + photometry_space="magnitude", ): """ Internal method for computing the ch2-squared value (for emcee). @@ -389,6 +426,7 @@ def log_likelihood_red_interpolated( z_min, z_max, True, + photometry_space=photometry_space, ) p = prior(*_x) diff --git a/test/test_diff2_functions_edges.py b/test/test_diff2_functions_edges.py index dbf02a2..bda02c6 100644 --- a/test/test_diff2_functions_edges.py +++ b/test/test_diff2_functions_edges.py @@ -22,6 +22,25 @@ def test_diff2_basic_shapes(): assert d2.shape == obs.shape and e2.shape == obs.shape +def test_diff2_flux_mode_shapes(): + interps = [_const_interp(10.0) for _ in range(3)] + model_flux = np.ones(3) * 10.0 ** (-0.4 * (10.0 + 5.0 * np.log10(10.0) - 5.0)) + obs = model_flux.copy() + err = np.ones(3) * 1e-4 + d2, e2 = diff2( + np.array([0.0]), + obs, + err, + 10.0, + 0.1, + interps, + True, + photometry_space="flux", + ) + assert d2.shape == obs.shape and e2.shape == obs.shape + assert np.isfinite(d2).all() and np.isfinite(e2).all() + + def test_diff2_distance_red_filter_invalid_distance_returns_inf(): interps = [_const_interp(10.0) for _ in range(2)] teff_itp = _const_interp(10000.0) diff --git a/test/test_fitter.py b/test/test_fitter.py index 3398855..d73d18d 100644 --- a/test/test_fitter.py +++ b/test/test_fitter.py @@ -77,8 +77,8 @@ def test_fitting_teff(mock_show): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", @@ -107,8 +107,8 @@ def test_fitting_teff_with_none(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, None], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, None], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], allow_none=True, atmosphere="H", logg=7.5, @@ -121,8 +121,8 @@ def test_fitting_teff_with_none(): assert np.isclose(ftr.results["H"].x, np.array([13000.0]), rtol=2.5e-02, atol=2.5e-02).all() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, 10.350], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 0.1], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, 10.350], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 0.1], allow_none=True, atmosphere="H", logg=7.5, @@ -135,8 +135,8 @@ def test_fitting_teff_with_none(): assert np.isclose(ftr.results["H"].x, np.array([13000.0]), rtol=2.5e-02, atol=2.5e-02).all() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, None], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, None], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], allow_none=True, atmosphere="H", logg=7.5, @@ -159,8 +159,8 @@ def test_fitting_logg_and_mbol(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", distance=10.0, @@ -192,8 +192,8 @@ def test_fitting_logg_teff_distance(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], @@ -219,8 +219,8 @@ def test_fitting_logg_teff_distance_nelder_mead(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], @@ -246,8 +246,8 @@ def test_fitting_teff_red(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", @@ -278,8 +278,8 @@ def test_fitting_logg_and_teff_red(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", distance=10.0, @@ -314,8 +314,8 @@ def test_fitting_logg_teff_distance_red(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], @@ -344,8 +344,8 @@ def test_fitting_logg_teff_distance_red_best_fit_plot_colour(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], @@ -385,8 +385,8 @@ def test_fitting_teff_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], logg=7.5, independent=["Teff"], distance=10.0, @@ -415,8 +415,8 @@ def test_fitting_teff_with_none_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, None], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, None], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], allow_none=True, atmosphere="H", logg=7.5, @@ -430,8 +430,8 @@ def test_fitting_teff_with_none_lsq(): assert np.isclose(ftr.results["H"].x, np.array([13000.0]), rtol=2.5e-02, atol=2.5e-02).all() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, 10.350], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 0.1], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, 10.350], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 0.1], allow_none=True, atmosphere="H", logg=7.5, @@ -444,8 +444,8 @@ def test_fitting_teff_with_none_lsq(): assert np.isclose(ftr.results["H"].x, np.array([13000.0]), rtol=2.5e-02, atol=2.5e-02).all() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, None], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, None], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], allow_none=True, atmosphere="H", logg=7.5, @@ -469,8 +469,8 @@ def test_fitting_logg_and_teff_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], distance=10.0, distance_err=0.1, @@ -503,8 +503,8 @@ def test_fitting_logg_teff_distance_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="CT", @@ -531,8 +531,8 @@ def test_fitting_logg_teff_distance_nelder_mead_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], initial_guess=[13000.0, 7.5, 10.0], method="least_squares", @@ -558,8 +558,8 @@ def test_fitting_teff_red_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], logg=7.5, independent=["Teff"], distance=10.0, @@ -591,8 +591,8 @@ def test_fitting_logg_and_teff_red_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], distance=10.0, distance_err=0.1, @@ -628,8 +628,8 @@ def test_fitting_logg_teff_distance_red_lsq(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], initial_guess=[13000.0, 7.5, 10.0], method="least_squares", @@ -663,12 +663,16 @@ def test_fitting_teff_emcee(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", method="emcee", + nwalkers=12, + nsteps=120, + nburns=20, + progress=False, distance=10.0, distance_err=0.1, initial_guess=[13000.0], @@ -705,14 +709,18 @@ def test_fitting_teff_with_none_emcee(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, None], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, None], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], allow_none=True, atmosphere="H", logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", method="emcee", + nwalkers=12, + nsteps=120, + nburns=20, + progress=False, distance=10.0, distance_err=0.1, refine=True, @@ -728,14 +736,18 @@ def test_fitting_teff_with_none_emcee(): ).all() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, 10.350], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 0.1], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, 10.350], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 0.1], allow_none=True, atmosphere="H", logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", method="emcee", + nwalkers=12, + nsteps=120, + nburns=20, + progress=False, distance=10.0, distance_err=0.1, refine=True, @@ -750,14 +762,18 @@ def test_fitting_teff_with_none_emcee(): ).all() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV", "U"], - mags=[10.882, 10.853, 10.946, 11.301, 11.183, None], - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], + photometry=[10.882, 10.853, 10.946, 11.301, 11.183, None], + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02, 10.0], allow_none=True, atmosphere="H", logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", method="emcee", + nwalkers=12, + nsteps=120, + nburns=20, + progress=False, distance=10.0, distance_err=0.1, refine=True, @@ -782,12 +798,16 @@ def test_fitting_teff_red_emcee(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], logg=7.5, independent=["Teff"], atmosphere_interpolator="CT", method="emcee", + nwalkers=12, + nsteps=120, + nburns=20, + progress=False, distance=10.0, distance_err=0.1, initial_guess=[13000.0], @@ -827,8 +847,8 @@ def test_chi2_minimization_red_interpolated(): ftr = WDfitter() ftr.fit( filters=["G3", "G3_BP", "G3_RP", "FUV", "NUV"], - mags=mags + extinction_interpolated, - mag_errors=[0.02, 0.02, 0.02, 0.02, 0.02], + photometry=mags + extinction_interpolated, + photometry_errors=[0.02, 0.02, 0.02, 0.02, 0.02], independent=["Teff", "logg"], atmosphere_interpolator="CT", method="least_squares", diff --git a/test/test_fitter_ct.py b/test/test_fitter_ct.py index d0a9556..6673320 100644 --- a/test/test_fitter_ct.py +++ b/test/test_fitter_ct.py @@ -70,8 +70,8 @@ def test_minimize_teff_logg(): ftr1 = WDfitter() ftr1.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="CT", @@ -103,8 +103,8 @@ def test_minimize_teff(): ftr2 = WDfitter() ftr2.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="CT", @@ -137,8 +137,8 @@ def test_minimize_teff_reddening(): ftr3 = WDfitter() ftr3.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="CT", @@ -173,8 +173,8 @@ def test_minimize_teff_reddening_interpolated(): ftr4 = WDfitter() ftr4.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", extinction_convolved=False, @@ -210,8 +210,8 @@ def test_minimize_teff_logg_reddening(): ftr5 = WDfitter() ftr5.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="CT", @@ -245,8 +245,8 @@ def test_minimize_teff_logg_reddening_interpolated(): ftr6 = WDfitter() ftr6.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", extinction_convolved=False, @@ -281,8 +281,8 @@ def test_minimize_teff_logg_distance(): ftr7 = WDfitter() ftr7.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="CT", @@ -312,8 +312,8 @@ def test_minimize_teff_distance(): ftr8 = WDfitter() ftr8.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="CT", @@ -344,8 +344,8 @@ def test_minimize_teff_distance_reddening(): ftr9 = WDfitter() ftr9.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="CT", @@ -378,8 +378,8 @@ def test_minimize_teff_distance_reddening_interpolated(): ftr10 = WDfitter() ftr10.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", extinction_convolved=False, @@ -413,8 +413,8 @@ def test_minimize_teff_logg_distance_reddening(): ftr11 = WDfitter() ftr11.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="CT", @@ -446,8 +446,8 @@ def test_minimize_teff_logg_distance_reddening_interpolated(): ftr12 = WDfitter() ftr12.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", extinction_convolved=False, @@ -483,8 +483,8 @@ def test_lsq_teff_logg(): ftr13 = WDfitter() ftr13.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="CT", @@ -516,8 +516,8 @@ def test_lsq_teff(): ftr14 = WDfitter() ftr14.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="CT", @@ -550,8 +550,8 @@ def test_lsq_teff_reddening(): ftr15 = WDfitter() ftr15.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="CT", @@ -586,8 +586,8 @@ def test_lsq_teff_reddening_interpolated(): ftr16 = WDfitter() ftr16.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", extinction_convolved=False, @@ -623,8 +623,8 @@ def test_lsq_teff_logg_reddening(): ftr17 = WDfitter() ftr17.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="CT", @@ -658,8 +658,8 @@ def test_lsq_teff_logg_reddening_interpolated(): ftr18 = WDfitter() ftr18.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", extinction_convolved=False, @@ -694,8 +694,8 @@ def test_lsq_teff_logg_distance(): ftr19 = WDfitter() ftr19.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="CT", @@ -725,8 +725,8 @@ def test_lsq_teff_distance(): ftr20 = WDfitter() ftr20.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="CT", @@ -757,8 +757,8 @@ def test_lsq_teff_distance_reddening(): ftr21 = WDfitter() ftr21.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="CT", @@ -791,8 +791,8 @@ def test_lsq_teff_distance_reddening_interpolated(): ftr22 = WDfitter() ftr22.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", extinction_convolved=False, @@ -826,8 +826,8 @@ def test_lsq_teff_logg_distance_reddening(): ftr23 = WDfitter() ftr23.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="CT", @@ -859,8 +859,8 @@ def test_lsq_teff_logg_distance_reddening_interpolated(): ftr24 = WDfitter() ftr24.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", extinction_convolved=False, @@ -899,16 +899,16 @@ def test_emcee_teff_logg(): ftr25.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, ) @@ -937,16 +937,16 @@ def test_emcee_teff(): ftr26.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -976,16 +976,16 @@ def test_emcee_teff_reddening(): ftr27.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -1017,17 +1017,17 @@ def test_emcee_teff_reddening_interpolated(): ftr28.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="CT", extinction_convolved=False, initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -1059,16 +1059,16 @@ def test_emcee_teff_logg_reddening(): ftr29.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, rv=RV, @@ -1099,17 +1099,17 @@ def test_emcee_teff_logg_reddening_interpolated(): ftr30.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", extinction_convolved=False, atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, rv=RV, @@ -1140,16 +1140,16 @@ def test_emcee_teff_logg_distance(): ftr31.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags, mags_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags, mags_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=50, - nsteps=5000, - nburns=500, + nwalkers=16, + nsteps=180, + nburns=30, ) assert np.isclose( ftr31.best_fit_params["H"]["Teff"], @@ -1176,16 +1176,16 @@ def test_emcee_teff_distance(): ftr32.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags, mags_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags, mags_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, ) assert np.isclose( @@ -1213,16 +1213,16 @@ def test_emcee_teff_distance_reddening(): ftr33.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, rv=RV, ebv=EBV, @@ -1252,22 +1252,22 @@ def test_emcee_teff_distance_reddening_interpolated(): ftr34.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate( + photometry=np.concatenate( ( mags + extinction_interpolated, mags_grizyJHK + extinction_grizyJHK_interpolated, ) ), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", extinction_convolved=False, atmosphere_interpolator="CT", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, rv=RV, ebv=EBV, @@ -1297,16 +1297,16 @@ def test_emcee_teff_logg_distance_reddening(): ftr35.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=50, - nsteps=5000, - nburns=500, + nwalkers=16, + nsteps=180, + nburns=30, rv=RV, ebv=EBV, ) @@ -1335,22 +1335,22 @@ def test_emcee_teff_logg_distance_reddening_interpolated(): ftr36.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate( + photometry=np.concatenate( ( mags + extinction_interpolated, mags_grizyJHK + extinction_grizyJHK_interpolated, ) ), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", extinction_convolved=False, atmosphere_interpolator="CT", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, rv=RV, ebv=EBV, ) diff --git a/test/test_fitter_extrpolation.py b/test/test_fitter_extrpolation.py index 37ca7ad..0d3fc51 100644 --- a/test/test_fitter_extrpolation.py +++ b/test/test_fitter_extrpolation.py @@ -20,8 +20,8 @@ def test_fitter_extrapolation_rbf(allow_extrapolation): f.fit( atmosphere=["H"], filters=filters, - mags=mags, - mag_errors=mag_errors, + photometry=mags, + photometry_errors=mag_errors, distance=10.0, distance_err=1.0, independent=["Mbol", "logg"], @@ -50,8 +50,8 @@ def test_fitter_extrapolation_ct(allow_extrapolation): f.fit( atmosphere=["H"], filters=filters, - mags=mags, - mag_errors=mag_errors, + photometry=mags, + photometry_errors=mag_errors, distance=10.0, distance_err=1.0, independent=["Mbol", "logg"], @@ -83,8 +83,8 @@ def test_fitter_extrapolation_rbf_initial_guess_outisde_grid(allow_extrapolation f.fit( atmosphere=["H"], filters=filters, - mags=mags, - mag_errors=mag_errors, + photometry=mags, + photometry_errors=mag_errors, distance=10.0, distance_err=1.0, independent=["Mbol", "logg"], @@ -115,8 +115,8 @@ def test_fitter_extrapolation_ct_initial_guess_outisde_grid(allow_extrapolation) f.fit( atmosphere=["H"], filters=filters, - mags=mags, - mag_errors=mag_errors, + photometry=mags, + photometry_errors=mag_errors, distance=10.0, distance_err=1.0, independent=["Mbol", "logg"], diff --git a/test/test_fitter_flux.py b/test/test_fitter_flux.py new file mode 100644 index 0000000..dab5212 --- /dev/null +++ b/test/test_fitter_flux.py @@ -0,0 +1,260 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""Tests for WD fitter in relative flux space.""" + +import numpy as np +import pytest + +from WDPhotTools.fitter import WDfitter + +FIVE_FILTERS = np.array(["G3", "G3_BP", "G3_RP", "FUV", "NUV"]) +MAGS = np.array([10.882, 10.853, 10.946, 11.301, 11.183], dtype=float) +MAG_ERRORS = np.ones(FIVE_FILTERS.size, dtype=float) * 0.02 +RELATIVE_FLUX = 10.0 ** (-0.4 * MAGS) +RELATIVE_FLUX_ERRORS = RELATIVE_FLUX * np.log(10.0) / 2.5 * MAG_ERRORS + + +def _run_fit(space, **kwargs): + ftr = WDfitter() + if space == "magnitude": + ftr.fit( + atmosphere="H", + filters=FIVE_FILTERS, + photometry=MAGS, + photometry_errors=MAG_ERRORS, + photometry_space=space, + **kwargs, + ) + else: + ftr.fit( + atmosphere="H", + filters=FIVE_FILTERS, + photometry=RELATIVE_FLUX, + photometry_errors=RELATIVE_FLUX_ERRORS, + photometry_space=space, + **kwargs, + ) + return ftr + + +def test_flux_minimize_teff(): + ftr = _run_fit( + "flux", + independent=["Teff"], + method="minimize", + atmosphere_interpolator="CT", + initial_guess=[13000.0], + logg=7.5, + distance=10.0, + distance_err=0.1, + ) + assert np.isclose( + ftr.best_fit_params["H"]["Teff"], + 13000.0, + rtol=1e-01, + atol=1e-01, + ) + + +def test_flux_least_squares_teff_logg(): + ftr = _run_fit( + "flux", + independent=["Teff", "logg"], + method="least_squares", + atmosphere_interpolator="CT", + initial_guess=[13000.0, 7.5], + distance=10.0, + distance_err=0.1, + ) + assert np.isclose( + ftr.best_fit_params["H"]["Teff"], + 13000.0, + rtol=1e-01, + atol=1e-01, + ) + assert np.isclose( + ftr.best_fit_params["H"]["logg"], + 7.5, + rtol=1e-01, + atol=1e-01, + ) + + +def test_flux_allow_none_and_plot_axis(): + photometry = [RELATIVE_FLUX[0], RELATIVE_FLUX[1], RELATIVE_FLUX[2], None, RELATIVE_FLUX[4]] + photometry_errors = [ + RELATIVE_FLUX_ERRORS[0], + RELATIVE_FLUX_ERRORS[1], + RELATIVE_FLUX_ERRORS[2], + 1.0, + RELATIVE_FLUX_ERRORS[4], + ] + ftr = WDfitter() + ftr.fit( + atmosphere="H", + filters=FIVE_FILTERS, + photometry=photometry, + photometry_errors=photometry_errors, + photometry_space="flux", + allow_none=True, + independent=["Teff"], + method="minimize", + atmosphere_interpolator="CT", + initial_guess=[13000.0], + logg=7.5, + distance=10.0, + distance_err=0.1, + ) + fig = ftr.show_best_fit(display=False, savefig=False, return_fig=True) + assert fig.gca().get_ylabel() == "Relative flux" + + +def test_flux_emcee_smoke(): + np.random.seed(42) + ftr = _run_fit( + "flux", + independent=["Teff"], + method="emcee", + atmosphere_interpolator="CT", + initial_guess=[13000.0], + logg=7.5, + distance=10.0, + distance_err=0.1, + nwalkers=20, + nsteps=80, + nburns=20, + progress=False, + ) + assert np.isfinite(ftr.best_fit_params["H"]["Teff"]) + + +def test_flux_input_validation(): + with pytest.raises(ValueError): + WDfitter().fit( + filters=FIVE_FILTERS, + photometry=RELATIVE_FLUX, + photometry_errors=RELATIVE_FLUX_ERRORS, + photometry_space="unknown", + independent=["Teff"], + initial_guess=[13000.0], + distance=10.0, + distance_err=0.1, + ) + + with pytest.raises(ValueError): + WDfitter().fit( + filters=FIVE_FILTERS, + photometry_space="flux", + independent=["Teff"], + initial_guess=[13000.0], + distance=10.0, + distance_err=0.1, + ) + + +def test_canonical_storage_fields(): + ftr = _run_fit( + "magnitude", + independent=["Teff"], + method="minimize", + atmosphere_interpolator="CT", + initial_guess=[13000.0], + logg=7.5, + distance=10.0, + distance_err=0.1, + ) + assert "photometry" in ftr.fitting_params + assert "photometry_errors" in ftr.fitting_params + assert "mags" not in ftr.fitting_params + assert "mag_errors" not in ftr.fitting_params + assert "fluxes" not in ftr.fitting_params + assert "flux_errors" not in ftr.fitting_params + assert len(ftr.best_fit_photometry["H"]) == FIVE_FILTERS.size + + +def test_legacy_keyword_arguments_removed(): + with pytest.raises(TypeError): + WDfitter().fit( + filters=FIVE_FILTERS, + mags=MAGS, + photometry_errors=MAG_ERRORS, + independent=["Teff"], + initial_guess=[13000.0], + distance=10.0, + distance_err=0.1, + ) + + +PARITY_CASES = [ + { + "name": "ct_minimize_teff", + "kwargs": { + "independent": ["Teff"], + "method": "minimize", + "atmosphere_interpolator": "CT", + "initial_guess": [13000.0], + "logg": 7.5, + "distance": 10.0, + "distance_err": 0.1, + }, + }, + { + "name": "ct_lsq_teff_logg", + "kwargs": { + "independent": ["Teff", "logg"], + "method": "least_squares", + "atmosphere_interpolator": "CT", + "initial_guess": [13000.0, 7.5], + "distance": 10.0, + "distance_err": 0.1, + }, + }, + { + "name": "ct_lsq_teff_reddening", + "kwargs": { + "independent": ["Teff"], + "method": "least_squares", + "atmosphere_interpolator": "CT", + "initial_guess": [13000.0], + "logg": 7.5, + "distance": 10.0, + "distance_err": 0.1, + "rv": 3.1, + "ebv": 0.123, + }, + }, + { + "name": "rbf_minimize_teff_logg_distance", + "kwargs": { + "independent": ["Teff", "logg"], + "method": "minimize", + "atmosphere_interpolator": "RBF", + "initial_guess": [13000.0, 7.5, 10.0], + "distance": None, + "distance_err": None, + }, + }, +] + + +@pytest.mark.parametrize("case", PARITY_CASES, ids=[c["name"] for c in PARITY_CASES]) +def test_flux_vs_magnitude_parameter_parity(case): + fit_mag = _run_fit("magnitude", **case["kwargs"]) + fit_flux = _run_fit("flux", **case["kwargs"]) + + mag = fit_mag.best_fit_params["H"] + flux = fit_flux.best_fit_params["H"] + + assert np.isclose(mag["Teff"], flux["Teff"], rtol=5e-03, atol=0.0) + assert np.isclose(mag["logg"], flux["logg"], atol=5e-03, rtol=0.0) + assert np.isclose(mag["distance"], flux["distance"], rtol=5e-03, atol=0.0) + assert np.isclose(mag["Mbol"], flux["Mbol"], atol=3e-02, rtol=0.0) + assert np.isfinite(mag["chi2"]) and np.isfinite(flux["chi2"]) + + # Avoid unstable relative comparison when chi2 values are effectively zero. + if max(abs(mag["chi2"]), abs(flux["chi2"])) > 1e-12: + assert abs(np.log10(abs(mag["chi2"])) - np.log10(abs(flux["chi2"]))) <= 1.0 + + assert len(fit_mag.best_fit_photometry["H"]) == FIVE_FILTERS.size + assert len(fit_flux.best_fit_photometry["H"]) == FIVE_FILTERS.size diff --git a/test/test_fitter_mass_ct.py b/test/test_fitter_mass_ct.py index 1b67f8c..5ddc9da 100644 --- a/test/test_fitter_mass_ct.py +++ b/test/test_fitter_mass_ct.py @@ -70,8 +70,8 @@ def test_minimize_teff_mass_reddening(): ftr5 = WDfitter() ftr5.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="minimize", atmosphere_interpolator="CT", @@ -105,8 +105,8 @@ def test_minimize_teff_mass_reddening_interpolated(): ftr6 = WDfitter() ftr6.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="minimize", extinction_convolved=False, @@ -141,8 +141,8 @@ def test_minimize_teff_mass_distance(): ftr7 = WDfitter() ftr7.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="minimize", atmosphere_interpolator="CT", @@ -172,8 +172,8 @@ def test_minimize_teff_mass_distance_reddening(): ftr11 = WDfitter() ftr11.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="minimize", atmosphere_interpolator="CT", @@ -205,8 +205,8 @@ def test_minimize_teff_mass_distance_reddening_interpolated(): ftr12 = WDfitter() ftr12.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="minimize", extinction_convolved=False, @@ -242,8 +242,8 @@ def test_lsq_teff_mass(): ftr13 = WDfitter() ftr13.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="least_squares", atmosphere_interpolator="CT", @@ -275,8 +275,8 @@ def test_lsq_teff_mass_reddening(): ftr17 = WDfitter() ftr17.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="least_squares", atmosphere_interpolator="CT", @@ -310,8 +310,8 @@ def test_lsq_teff_mass_reddening_interpolated(): ftr18 = WDfitter() ftr18.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="least_squares", extinction_convolved=False, @@ -346,8 +346,8 @@ def test_lsq_teff_mass_distance(): ftr19 = WDfitter() ftr19.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="least_squares", atmosphere_interpolator="CT", @@ -377,8 +377,8 @@ def test_lsq_teff_mass_distance_reddening(): ftr23 = WDfitter() ftr23.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="least_squares", atmosphere_interpolator="CT", @@ -410,8 +410,8 @@ def test_lsq_teff_mass_distance_reddening_interpolated(): ftr24 = WDfitter() ftr24.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "mass"], method="least_squares", extinction_convolved=False, diff --git a/test/test_fitter_rbf.py b/test/test_fitter_rbf.py index 3408fbf..da10544 100644 --- a/test/test_fitter_rbf.py +++ b/test/test_fitter_rbf.py @@ -70,8 +70,8 @@ def test_minimize_teff_logg(): ftr1 = WDfitter() ftr1.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="RBF", @@ -93,8 +93,8 @@ def test_minimize_teff(): ftr2 = WDfitter() ftr2.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="RBF", @@ -117,8 +117,8 @@ def test_minimize_teff_reddening(): ftr3 = WDfitter() ftr3.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="RBF", @@ -143,8 +143,8 @@ def test_minimize_teff_reddening_interpolated(): ftr4 = WDfitter() ftr4.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", extinction_convolved=False, @@ -170,8 +170,8 @@ def test_minimize_teff_logg_reddening(): ftr5 = WDfitter() ftr5.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="RBF", @@ -195,8 +195,8 @@ def test_minimize_teff_logg_reddening_interpolated(): ftr6 = WDfitter() ftr6.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", extinction_convolved=False, @@ -221,8 +221,8 @@ def test_minimize_teff_logg_distance(): ftr7 = WDfitter() ftr7.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="RBF", @@ -242,8 +242,8 @@ def test_minimize_teff_distance(): ftr8 = WDfitter() ftr8.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="RBF", @@ -264,8 +264,8 @@ def test_minimize_teff_distance_reddening(): ftr9 = WDfitter() ftr9.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", atmosphere_interpolator="RBF", @@ -288,8 +288,8 @@ def test_minimize_teff_distance_reddening_interpolated(): ftr10 = WDfitter() ftr10.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="minimize", extinction_convolved=False, @@ -313,8 +313,8 @@ def test_minimize_teff_logg_distance_reddening(): ftr11 = WDfitter() ftr11.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", atmosphere_interpolator="RBF", @@ -336,8 +336,8 @@ def test_minimize_teff_logg_distance_reddening_interpolated(): ftr12 = WDfitter() ftr12.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="minimize", extinction_convolved=False, @@ -366,8 +366,8 @@ def test_lsq_teff_logg(): ftr13 = WDfitter() ftr13.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="RBF", @@ -389,8 +389,8 @@ def test_lsq_teff(): ftr14 = WDfitter() ftr14.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="RBF", @@ -413,8 +413,8 @@ def test_lsq_teff_reddening(): ftr15 = WDfitter() ftr15.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="RBF", @@ -439,8 +439,8 @@ def test_lsq_teff_reddening_interpolated(): ftr16 = WDfitter() ftr16.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", extinction_convolved=False, @@ -466,8 +466,8 @@ def test_lsq_teff_logg_reddening(): ftr17 = WDfitter() ftr17.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="RBF", @@ -491,8 +491,8 @@ def test_lsq_teff_logg_reddening_interpolated(): ftr18 = WDfitter() ftr18.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", extinction_convolved=False, @@ -517,8 +517,8 @@ def test_lsq_teff_logg_distance(): ftr19 = WDfitter() ftr19.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="RBF", @@ -538,8 +538,8 @@ def test_lsq_teff_distance(): ftr20 = WDfitter() ftr20.fit( filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="RBF", @@ -560,8 +560,8 @@ def test_lsq_teff_distance_reddening(): ftr21 = WDfitter() ftr21.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", atmosphere_interpolator="RBF", @@ -584,8 +584,8 @@ def test_lsq_teff_distance_reddening_interpolated(): ftr22 = WDfitter() ftr22.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="least_squares", extinction_convolved=False, @@ -609,8 +609,8 @@ def test_lsq_teff_logg_distance_reddening(): ftr23 = WDfitter() ftr23.fit( filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", atmosphere_interpolator="RBF", @@ -632,8 +632,8 @@ def test_lsq_teff_logg_distance_reddening_interpolated(): ftr24 = WDfitter() ftr24.fit( filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="least_squares", extinction_convolved=False, @@ -660,16 +660,16 @@ def test_emcee_teff_logg(): ftr25.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, ) @@ -688,16 +688,16 @@ def test_emcee_teff(): ftr26.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -717,16 +717,16 @@ def test_emcee_teff_reddening(): ftr27.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -748,17 +748,17 @@ def test_emcee_teff_reddening_interpolated(): ftr28.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", extinction_convolved=False, initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -780,16 +780,16 @@ def test_emcee_teff_logg_reddening(): ftr29.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, rv=RV, @@ -810,17 +810,17 @@ def test_emcee_teff_logg_reddening_interpolated(): ftr30.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", extinction_convolved=False, atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, rv=RV, @@ -841,16 +841,16 @@ def test_emcee_teff_logg_distance(): ftr31.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags, mags_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags, mags_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, ) results31 = ftr31.best_fit_params["H"]["Teff"] assert np.isclose(results31, 13000.0, rtol=1e-01, atol=1e-01) @@ -867,16 +867,16 @@ def test_emcee_teff_distance(): ftr32.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags, mags_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags, mags_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, ) results32 = ftr32.best_fit_params["H"]["Teff"] @@ -894,16 +894,16 @@ def test_emcee_teff_distance_reddening(): ftr33.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, rv=RV, ebv=EBV, @@ -923,22 +923,22 @@ def test_emcee_teff_distance_reddening_interpolated(): ftr34.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate( + photometry=np.concatenate( ( mags + extinction_interpolated, mags_grizyJHK + extinction_grizyJHK_interpolated, ) ), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", extinction_convolved=False, atmosphere_interpolator="RBF", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, rv=RV, ebv=EBV, @@ -958,16 +958,16 @@ def test_emcee_teff_logg_distance_reddening(): ftr35.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, rv=RV, ebv=EBV, ) @@ -986,22 +986,22 @@ def test_emcee_teff_logg_distance_reddening_interpolated(): ftr36.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate( + photometry=np.concatenate( ( mags + extinction_interpolated, mags_grizyJHK + extinction_grizyJHK_interpolated, ) ), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", extinction_convolved=False, atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, rv=RV, ebv=EBV, ) diff --git a/test/test_fitter_rbf_emcee_prior.py b/test/test_fitter_rbf_emcee_prior.py index d5c34b9..753d3ef 100644 --- a/test/test_fitter_rbf_emcee_prior.py +++ b/test/test_fitter_rbf_emcee_prior.py @@ -111,16 +111,16 @@ def test_emcee_teff_logg(): ftr01.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, prior=prior_teff_logg, @@ -150,16 +150,16 @@ def test_emcee_teff(): ftr26.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -190,16 +190,16 @@ def test_emcee_teff_reddening(): ftr27.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -232,17 +232,17 @@ def test_emcee_teff_reddening_interpolated(): ftr28.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", extinction_convolved=False, initial_guess=[13000.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, distance=10.0, distance_err=0.1, @@ -275,16 +275,16 @@ def test_emcee_teff_logg_reddening(): ftr29.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, rv=RV, @@ -316,17 +316,17 @@ def test_emcee_teff_logg_reddening_interpolated(): ftr30.fit( atmosphere="H", filters=five_filters_name_list, - mags=mags + extinction_interpolated, - mag_errors=np.ones(five_filters_name_list.size) * 0.02, + photometry=mags + extinction_interpolated, + photometry_errors=np.ones(five_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", extinction_convolved=False, atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, distance=10.0, distance_err=0.1, rv=RV, @@ -358,16 +358,16 @@ def test_emcee_teff_logg_distance(): ftr31.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags, mags_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags, mags_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, prior=prior_teff_logg_distance, ) assert np.isclose( @@ -395,16 +395,16 @@ def test_emcee_teff_distance(): ftr32.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags, mags_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags, mags_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, prior=prior_teff_distance, ) @@ -433,16 +433,16 @@ def test_emcee_teff_distance_reddening(): ftr33.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, rv=RV, ebv=EBV, @@ -473,22 +473,22 @@ def test_emcee_teff_distance_reddening_interpolated(): ftr34.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate( + photometry=np.concatenate( ( mags + extinction_interpolated, mags_grizyJHK + extinction_grizyJHK_interpolated, ) ), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff"], method="emcee", extinction_convolved=False, atmosphere_interpolator="RBF", initial_guess=[13000.0, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, logg=7.5, rv=RV, ebv=EBV, @@ -519,16 +519,16 @@ def test_emcee_teff_logg_distance_reddening(): ftr35.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry=np.concatenate((mags + extinction, mags_grizyJHK + extinction_grizyJHK)), + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, rv=RV, ebv=EBV, prior=prior_teff_logg_distance, @@ -558,22 +558,22 @@ def test_emcee_teff_logg_distance_reddening_interpolated(): ftr36.fit( atmosphere="H", filters=thirteen_filters_name_list, - mags=np.concatenate( + photometry=np.concatenate( ( mags + extinction_interpolated, mags_grizyJHK + extinction_grizyJHK_interpolated, ) ), - mag_errors=np.ones(thirteen_filters_name_list.size) * 0.02, + photometry_errors=np.ones(thirteen_filters_name_list.size) * 0.02, independent=["Teff", "logg"], method="emcee", extinction_convolved=False, atmosphere_interpolator="RBF", initial_guess=[13000.0, 7.5, 10.0], refine=False, - nwalkers=20, - nsteps=1000, - nburns=100, + nwalkers=12, + nsteps=120, + nburns=20, rv=RV, ebv=EBV, prior=prior_teff_logg_distance, diff --git a/test/test_fitter_regression_main.py b/test/test_fitter_regression_main.py new file mode 100644 index 0000000..3ac9bf6 --- /dev/null +++ b/test/test_fitter_regression_main.py @@ -0,0 +1,133 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""Regression checks against pinned origin/main fitter results. + +Baseline values were generated from origin/main at commit 59f833d. +""" + +import numpy as np +import pytest + +from WDPhotTools.fitter import WDfitter + +FILTERS = np.array(["G3", "G3_BP", "G3_RP", "FUV", "NUV"]) +MAGS = np.array([10.882, 10.853, 10.946, 11.301, 11.183], dtype=float) +MAG_ERRORS = np.ones(FILTERS.size, dtype=float) * 0.02 + + +REGRESSION_CASES = [ + { + "name": "ct_min_teff_fixed_dist", + "kwargs": { + "atmosphere": "H", + "filters": FILTERS, + "photometry": MAGS, + "photometry_errors": MAG_ERRORS, + "photometry_space": "magnitude", + "logg": 7.5, + "independent": ["Teff"], + "atmosphere_interpolator": "CT", + "method": "minimize", + "distance": 10.0, + "distance_err": 0.1, + "initial_guess": [13000.0], + }, + "expected": { + "Teff": 13000.0, + "logg": 7.5, + "distance": 10.0, + "Mbol": 9.962, + "chi2": 0.0, + }, + }, + { + "name": "ct_lsq_teff_logg_fixed_dist", + "kwargs": { + "atmosphere": "H", + "filters": FILTERS, + "photometry": MAGS, + "photometry_errors": MAG_ERRORS, + "photometry_space": "magnitude", + "independent": ["Teff", "logg"], + "atmosphere_interpolator": "CT", + "method": "least_squares", + "distance": 10.0, + "distance_err": 0.1, + "initial_guess": [13000.0, 7.5], + }, + "expected": { + "Teff": 13000.0, + "logg": 7.5, + "distance": 10.0, + "Mbol": 9.962, + "chi2": 0.0, + }, + }, + { + "name": "rbf_min_teff_logg_dist", + "kwargs": { + "atmosphere": "H", + "filters": FILTERS, + "photometry": MAGS, + "photometry_errors": MAG_ERRORS, + "photometry_space": "magnitude", + "independent": ["Teff", "logg"], + "atmosphere_interpolator": "RBF", + "method": "minimize", + "initial_guess": [13000.0, 7.5, 10.0], + }, + "expected": { + "Teff": 13000.000000000116, + "logg": 7.5, + "distance": 10.0, + "Mbol": 9.96199999999996, + "chi2": 2.4194140247106092e-21, + }, + }, + { + "name": "ct_lsq_teff_reddening", + "kwargs": { + "atmosphere": "H", + "filters": FILTERS, + "photometry": MAGS, + "photometry_errors": MAG_ERRORS, + "photometry_space": "magnitude", + "logg": 7.5, + "independent": ["Teff"], + "atmosphere_interpolator": "CT", + "method": "least_squares", + "distance": 10.0, + "distance_err": 0.1, + "initial_guess": [13000.0], + "rv": 3.1, + "ebv": 0.123, + }, + "expected": { + "Teff": 15562.152043762782, + "logg": 7.5, + "distance": 10.0, + "Mbol": 9.152973607611566, + "chi2": 44.26036599996003, + }, + }, +] + + +@pytest.mark.parametrize("case", REGRESSION_CASES, ids=[c["name"] for c in REGRESSION_CASES]) +def test_regression_against_origin_main(case): + ftr = WDfitter() + ftr.fit(**case["kwargs"]) + params = ftr.best_fit_params["H"] + expected = case["expected"] + + assert np.isclose(params["Teff"], expected["Teff"], rtol=1e-02, atol=0.0) + assert np.isclose(params["logg"], expected["logg"], atol=1e-02, rtol=0.0) + assert np.isclose(params["distance"], expected["distance"], rtol=1e-02, atol=0.0) + assert np.isclose(params["Mbol"], expected["Mbol"], atol=1e-02, rtol=0.0) + + assert np.isfinite(params["chi2"]) + if expected["chi2"] <= 1e-20: + assert abs(params["chi2"]) <= 1e-6 + else: + assert abs(np.log10(abs(params["chi2"])) - np.log10(abs(expected["chi2"]))) <= 1.0 diff --git a/test/test_interpolator_extrapolation_success.py b/test/test_interpolator_extrapolation_success.py new file mode 100644 index 0000000..d9fa432 --- /dev/null +++ b/test/test_interpolator_extrapolation_success.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""Extrapolation quality checks for atmosphere and cooling interpolators.""" + +import inspect + +import numpy as np +import pytest + +from WDPhotTools.atmosphere_model_reader import AtmosphereModelReader +from WDPhotTools.cooling_model_reader import CoolingModelReader + + +def _sample_with_fraction(rng, value_min, value_max, fraction, size): + span = value_max - value_min + return rng.uniform(value_min - span * fraction, value_max + span * fraction, size=size) + + +def _estimate_atmosphere_success_rates(interpolator, sample_size=200, seed=12345): + rng = np.random.default_rng(seed) + atm = AtmosphereModelReader() + itp = atm.interp_am( + dependent="Teff", + atmosphere="H", + independent=["logg", "Mbol"], + interpolator=interpolator, + allow_extrapolation=True, + ) + + logg_grid = np.asarray(atm.model_da["logg"], dtype=float) + mbol_grid = np.asarray(atm.model_da["Mbol"], dtype=float) + logg_min, logg_max = float(np.nanmin(logg_grid)), float(np.nanmax(logg_grid)) + mbol_min, mbol_max = float(np.nanmin(mbol_grid)), float(np.nanmax(mbol_grid)) + + rates = {} + for fraction in (0.1, 0.2, 0.3, 0.4, 0.5): + logg = _sample_with_fraction(rng, logg_min, logg_max, fraction, sample_size) + mbol = _sample_with_fraction(rng, mbol_min, mbol_max, fraction, sample_size) + extrapolated = (logg < logg_min) | (logg > logg_max) | (mbol < mbol_min) | (mbol > mbol_max) + teff = np.asarray(itp(logg, mbol), dtype=float).reshape(-1) + valid = np.isfinite(teff) & (teff > 0.0) + rates[fraction] = float(np.mean(valid[extrapolated])) if np.any(extrapolated) else 1.0 + + return rates + + +def _estimate_cooling_success_rates(interpolator, sample_size=200, seed=12345): + rng = np.random.default_rng(seed) + cmr = CoolingModelReader() + cmr.compute_cooling_age_interpolator( + interpolator=interpolator, + allow_extrapolation=True, + ) + + logl_grid = np.log10(np.asarray(cmr.luminosity, dtype=float)) + mass_grid = np.asarray(cmr.mass, dtype=float) + logl_min, logl_max = float(np.nanmin(logl_grid)), float(np.nanmax(logl_grid)) + mass_min, mass_max = float(np.nanmin(mass_grid)), float(np.nanmax(mass_grid)) + + rates = {} + for fraction in (0.1, 0.2, 0.3, 0.4, 0.5): + logl = _sample_with_fraction(rng, logl_min, logl_max, fraction, sample_size) + mass = _sample_with_fraction(rng, mass_min, mass_max, fraction, sample_size) + extrapolated = (logl < logl_min) | (logl > logl_max) | (mass < mass_min) | (mass > mass_max) + cooling_age = np.asarray(cmr.cooling_interpolator(logl, mass), dtype=float).reshape(-1) + cooling_rate = np.asarray(cmr.cooling_rate_interpolator(logl, mass), dtype=float).reshape(-1) + valid = np.isfinite(cooling_age) & np.isfinite(cooling_rate) & (cooling_age > 0.0) & (cooling_rate <= 0.0) + rates[fraction] = float(np.mean(valid[extrapolated])) if np.any(extrapolated) else 1.0 + + return rates + + +def test_extrapolation_defaults_remain_disabled(): + """Ensure extrapolation defaults remain disabled for public APIs.""" + interp_am_sig = inspect.signature(AtmosphereModelReader.interp_am) + compute_sig = inspect.signature(CoolingModelReader.compute_cooling_age_interpolator) + assert interp_am_sig.parameters["allow_extrapolation"].default is False + assert compute_sig.parameters["allow_extrapolation"].default is False + + +@pytest.mark.parametrize("interpolator", ["CT", "RBF"]) +def test_atmosphere_extrapolation_success_rate(interpolator): + """10% atmosphere extrapolation must return finite and physically valid Teff.""" + rates = _estimate_atmosphere_success_rates(interpolator) + assert rates[0.1] >= 0.99 + assert rates[0.5] >= 0.95 + + +@pytest.mark.parametrize("interpolator", ["CT", "RBF"]) +def test_cooling_extrapolation_success_rate(interpolator): + """10% cooling extrapolation must return finite positive age and finite non-positive dL/dt.""" + rates = _estimate_cooling_success_rates(interpolator) + assert rates[0.1] >= 0.99 + assert rates[0.5] >= 0.95 + + +@pytest.mark.parametrize("interpolator", ["CT", "RBF"]) +def test_atmosphere_boundary_behaviour_unchanged(interpolator): + """Boundary and out-of-bounds behaviour stays unchanged when extrapolation is disabled.""" + atm = AtmosphereModelReader() + interp_no_extrap = atm.interp_am( + dependent="Teff", + atmosphere="H", + independent=["logg", "Mbol"], + interpolator=interpolator, + allow_extrapolation=False, + ) + interp_with_extrap = atm.interp_am( + dependent="Teff", + atmosphere="H", + independent=["logg", "Mbol"], + interpolator=interpolator, + allow_extrapolation=True, + ) + + logg_grid = np.asarray(atm.model_da["logg"], dtype=float) + mbol_grid = np.asarray(atm.model_da["Mbol"], dtype=float) + logg_min, logg_max = float(np.nanmin(logg_grid)), float(np.nanmax(logg_grid)) + mbol_min, mbol_max = float(np.nanmin(mbol_grid)), float(np.nanmax(mbol_grid)) + + boundary_logg = np.array([logg_min, logg_min, logg_max, logg_max], dtype=float) + boundary_mbol = np.array([mbol_min, mbol_max, mbol_min, mbol_max], dtype=float) + no_extrap_boundary = np.asarray(interp_no_extrap(boundary_logg, boundary_mbol), dtype=float).reshape(-1) + with_extrap_boundary = np.asarray(interp_with_extrap(boundary_logg, boundary_mbol), dtype=float).reshape(-1) + finite = np.isfinite(no_extrap_boundary) & np.isfinite(with_extrap_boundary) + assert np.allclose(no_extrap_boundary[finite], with_extrap_boundary[finite], rtol=1e-10, atol=1e-6) + + out_of_bounds_logg = np.array([logg_min - 0.1, logg_max + 0.1], dtype=float) + out_of_bounds_mbol = np.array([mbol_min - 0.1, mbol_max + 0.1], dtype=float) + no_extrap_oob = np.asarray(interp_no_extrap(out_of_bounds_logg, out_of_bounds_mbol), dtype=float).reshape(-1) + assert np.all(np.isneginf(no_extrap_oob)) + + +@pytest.mark.parametrize("interpolator", ["CT", "RBF"]) +def test_cooling_boundary_behaviour_unchanged(interpolator): + """Cooling interpolator boundary behaviour stays unchanged when extrapolation is disabled.""" + cmr_no_extrap = CoolingModelReader() + cmr_no_extrap.compute_cooling_age_interpolator(interpolator=interpolator, allow_extrapolation=False) + cmr_with_extrap = CoolingModelReader() + cmr_with_extrap.compute_cooling_age_interpolator(interpolator=interpolator, allow_extrapolation=True) + + lum_grid = np.log10(np.asarray(cmr_no_extrap.luminosity, dtype=float)) + mass_grid = np.asarray(cmr_no_extrap.mass, dtype=float) + lum_min, lum_max = float(np.nanmin(lum_grid)), float(np.nanmax(lum_grid)) + mass_min, mass_max = float(np.nanmin(mass_grid)), float(np.nanmax(mass_grid)) + + boundary_lum = np.array([lum_min, lum_min, lum_max, lum_max], dtype=float) + boundary_mass = np.array([mass_min, mass_max, mass_min, mass_max], dtype=float) + no_extrap_boundary = np.asarray( + cmr_no_extrap.cooling_interpolator(boundary_lum, boundary_mass), + dtype=float, + ).reshape(-1) + with_extrap_boundary = np.asarray( + cmr_with_extrap.cooling_interpolator(boundary_lum, boundary_mass), + dtype=float, + ).reshape(-1) + finite = np.isfinite(no_extrap_boundary) & np.isfinite(with_extrap_boundary) + assert np.allclose(no_extrap_boundary[finite], with_extrap_boundary[finite], rtol=1e-10, atol=1e-6) + + out_of_bounds_lum = np.array([lum_min - 0.05, lum_max + 0.05], dtype=float) + out_of_bounds_mass = np.array([mass_min - 0.05, mass_max + 0.05], dtype=float) + no_extrap_oob = np.asarray( + cmr_no_extrap.cooling_interpolator(out_of_bounds_lum, out_of_bounds_mass), + dtype=float, + ).reshape(-1) + assert np.all(np.isneginf(no_extrap_oob)) diff --git a/test/test_plotter.py b/test/test_plotter.py index 96d08c2..6193ca2 100644 --- a/test/test_plotter.py +++ b/test/test_plotter.py @@ -60,9 +60,7 @@ def test_plot_atmosphere_model(mock_show): @patch("matplotlib.pyplot.show") def test_plot_atmosphere_model_different_filters(mock_show): - plotter.plot_atmosphere_model( - x="B-V", y="U", invert_yaxis=True, display=True - ) + plotter.plot_atmosphere_model(x="B-V", y="U", invert_yaxis=True, display=True) @patch("matplotlib.pyplot.show") @@ -164,7 +162,9 @@ def test_plot_atmosphere_models_none_folder_savefig(mock_show): _folder = os.getcwd() for e in ["png"]: _filename = "test_plot_atmosphere_models_none_folder_savefig" + "." + e - assert os.path.isfile(os.path.join(_folder, _filename)) + _filepath = os.path.join(_folder, _filename) + assert os.path.isfile(_filepath) + os.remove(_filepath) # Testing plot_cooling_model savefig diff --git a/test/test_readme_figures.py b/test/test_readme_figures.py new file mode 100644 index 0000000..61d46bc --- /dev/null +++ b/test/test_readme_figures.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""Smoke tests for figure-generation snippets in README.""" + +from pathlib import Path + +import matplotlib +import numpy as np +from matplotlib import pyplot as plt + +from WDPhotTools import plotter +from WDPhotTools.atmosphere_model_reader import AtmosphereModelReader + +matplotlib.use("Agg") + + +def test_readme_plotter_atmosphere_figure(tmp_path): + output = Path(tmp_path) + plotter.plot_atmosphere_model( + invert_yaxis=True, + display=False, + savefig=True, + folder=str(output), + filename="DA_cooling_tracks_from_plotter", + ext="png", + ) + assert (output / "DA_cooling_tracks_from_plotter.png").is_file() + + +def test_readme_custom_atmosphere_figure(tmp_path): + atm = AtmosphereModelReader() + + g_band = atm.interp_am() + bp_band = atm.interp_am(dependent="G3_BP") + rp_band = atm.interp_am(dependent="G3_RP") + + logg = np.arange(7.0, 9.5, 0.5) + mbol = np.arange(0.0, 20.0, 0.1) + + fig = plt.figure(figsize=(8, 8)) + ax = fig.gca() + for value in logg: + ax.plot( + bp_band(value, mbol) - rp_band(value, mbol), + g_band(value, mbol), + ) + ax.set_ylim(20.0, 6.0) + ax.grid() + ax.set_xlabel("(BP - RP) / mag") + ax.set_ylabel("G / mag") + ax.set_title("DA Cooling tracks") + fig.tight_layout() + + output = Path(tmp_path) / "DA_cooling_tracks.png" + fig.savefig(output) + plt.close(fig) + assert output.is_file() + + +def test_readme_plotter_cooling_figure(tmp_path): + output = Path(tmp_path) + plotter.plot_cooling_model( + mass=[0.2, 0.4, 0.6, 0.8, 1.0], + display=False, + savefig=True, + folder=str(output), + filename="DA_cooling_model_from_plotter", + ext="png", + ) + assert (output / "DA_cooling_model_from_plotter.png").is_file()