20978266610

Committed 14 Jan 2026 12:54AM UTC coverage: 95.519% (-0.6%) from 96.096%

Build # 20978266610

Build Type

push

github

Committed by

web-flow

Commit Message

v0.0.13 (#49)

* handling numpy 2 (#47)

* Compat: use numpy generic types for v1/v2 (np.floating/np.integer) in src; no stylistic changes

* Allow NumPy 2 in install: relax to >=1.23,<3; update README to state v1+v2 support

* updated array handling

* removed test output figure

* updated setup.cfg to include numpy2

* use asarray

* use np.where

* use where instead of argwhere; wrapped interpolator outputs with np.asarray(...).itme() or reshape(-1)[0].

* pack values in correct array shape for atmosphere interpolator

* Tests/coverage improvements (#48)

* tests: increase coverage and add edge-case tests across modules

* added more tests

* final fixes to the minimisation function and the interpolation functions. updated example scripts and tests.

* updated test

* avoid extrapolation in both ct and rbf interpolator for the atmosphere model

* fixed extrapolation issue

Run Details

299 of 337 new or added lines in 11 files covered. (88.72%)

4 existing lines in 3 files now uncovered.

3389 of 3548 relevant lines covered (95.52%)

1.91 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

90.64

/src/WDPhotTools/atmosphere_model_reader.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""Handling the formatting of different atmosphere models"""

import os

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.interpolate import RBFInterpolator


class AtmosphereModelReader(object):
    """Handling the formatting of different atmosphere models"""

    def __init__(self):
        super(AtmosphereModelReader, self).__init__()

        self.this_file = os.path.dirname(os.path.abspath(__file__))

        self.model_list = {
            "montreal_co_da_20": "Bedard et al. 2020 CO DA",
            "montreal_co_db_20": "Bedard et al. 2020 CO DB",
            "lpcode_he_da_07": "Panei et al. 2007 He DA",
            "lpcode_co_da_07": "Panei et al. 2007 CO DA",
            "lpcode_he_da_09": "Althaus et al. 2009 He DA",
            "lpcode_co_da_10_z001": "Renedo et al. 2010 CO DA Z=0.01",
            "lpcode_co_da_10_z0001": "Renedo et al. 2010 CO DA Z=0.001",
            "lpcode_co_da_15_z00003": "Althaus et al. 2015 DA Z=0.00003",
            "lpcode_co_da_15_z0001": "Althaus et al. 2015 DA Z=0.0001",
            "lpcode_co_da_15_z0005": "Althaus et al. 2015 DA Z=0.0005",
            "lpcode_co_db_17_z00005": "Althaus et al. 2017 DB Y=0.4",
            "lpcode_co_db_17_z0001": "Althaus et al. 2017 DB Y=0.4",
            "lpcode_co_db_17": "Camisassa et al. 2017 DB",
            "lpcode_one_da_07": "Althaus et al. 2007 ONe DA",
            "lpcode_one_da_19": "Camisassa et al. 2019 ONe DA",
            "lpcode_one_db_19": "Camisassa et al. 2019 ONe DB",
            "lpcode_da_22": "Althaus et al. 2013 He DA, Camisassa et al. 2016 CO DA,  Camisassa et al. 2019 ONe DA",
            "lpcode_db_22": "Camisassa et al. 2017 CO DB, " + "Camisassa et al. 2019 ONe DB",
        }

        # DA atmosphere
        filepath_da = os.path.join(
            os.path.dirname(os.path.abspath(__file__)),
            "wd_photometry/Table_DA_13012021.txt",
        )

        # DB atmosphere
        filepath_db = os.path.join(
            os.path.dirname(os.path.abspath(__file__)),
            "wd_photometry/Table_DB_13012021.txt",
        )

        # Prepare the array column dtype
        self.column_key = np.array(
            [
                "Teff",
                "logg",
                "mass",
                "Mbol",
                "BC",
                "U",
                "B",
                "V",
                "R",
                "I",
                "J",
                "H",
                "Ks",
                "Y_mko",
                "J_mko",
                "H_mko",
                "K_mko",
                "W1",
                "W2",
                "W3",
                "W4",
                "S36",
                "S45",
                "S58",
                "S80",
                "u_sdss",
                "g_sdss",
                "r_sdss",
                "i_sdss",
                "z_sdss",
                "g_ps1",
                "r_ps1",
                "i_ps1",
                "z_ps1",
                "y_ps1",
                "G2",
                "G2_BP",
                "G2_RP",
                "G3",
                "G3_BP",
                "G3_RP",
                "FUV",
                "NUV",
                "age",
            ]
        )
        self.column_key_formatted = np.array(
            [
                r"T$_{\mathrm{eff}}$",
                "log(g)",
                "Mass",
                r"M$_{\mathrm{bol}}$",
                "BC",
                r"$U$",
                r"$B$",
                r"$V$",
                r"$R$",
                r"$I$",
                r"$J$",
                r"$H$",
                r"$K_{\mathrm{s}}$",
                r"$Y_{\mathrm{MKO}}$",
                r"$J_{\mathrm{MKO}}$",
                r"$H_{\mathrm{MKO}}$",
                r"$K_{\mathrm{MKO}}$",
                r"$W_{1}$",
                r"$W_{2}$",
                r"$W_{3}$",
                r"$W_{4}$",
                r"$S_{36}$",
                r"$S_{45}$",
                r"$S_{58}$",
                r"$S_{80}$",
                r"u$_{\mathrm{SDSS}}$",
                r"$g_{\mathrm{SDSS}}$",
                r"$r_{\mathrm{SDSS}}$",
                r"$i_{\mathrm{SDSS}}$",
                r"$z_{\mathrm{SDSS}}$",
                r"$g_{\mathrm{PS1}}$",
                r"$r_{\mathrm{PS1}}$",
                r"$i_{\mathrm{PS1}}$",
                r"$z_{\mathrm{PS1}}$",
                r"$y_{\mathrm{PS1}}$",
                r"$G_{\mathrm{DR2}}$",
                r"$G_{\mathrm{BP, DR2}}$",
                r"$G_{\mathrm{RP, DR2}}$",
                r"$G_{\mathrm{DR3}}$",
                r"$G_{\mathrm{BP, DR3}}$",
                r"$G_{\mathrm{RP, DR3}}$",
                "FUV",
                "NUV",
                "Age",
            ]
        )
        self.column_key_unit = np.array(
            [
                "K",
                r"(cm/s$^2$)",
                r"M$_\odot$",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "mag",
                "yr",
            ]
        )
        self.column_key_wavelength = np.array(
            [
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                3585.0,
                4371.0,
                5478.0,
                6504.0,
                8020.0,
                12350.0,
                16460.0,
                21600.0,
                10310.0,
                12500.0,
                16360.0,
                22060.0,
                33682.0,
                46179.0,
                120717.0,
                221944.0,
                35378.0,
                44780.0,
                56962.0,
                77978.0,
                3557.0,
                4702.0,
                6175.0,
                7491.0,
                8946.0,
                4849.0,
                6201.0,
                7535.0,
                8674.0,
                9628.0,
                6229.0,
                5037.0,
                7752.0,
                6218.0,
                5110.0,
                7769.0,
                1535.0,
                2301.0,
                0.0,
            ]
        )

        self.column_names = {}
        self.column_units = {}
        self.column_wavelengths = {}
        for i, j, k, _l in zip(
            self.column_key,
            self.column_key_formatted,
            self.column_key_unit,
            self.column_key_wavelength,
        ):
            self.column_names[i] = j
            self.column_units[i] = k
            self.column_wavelengths[i] = _l

        self.column_type = np.array(([np.float64] * len(self.column_key)))
        self.dtype = list(zip(self.column_key, self.column_type))

        # Load the synthetic photometry file in a recarray
        self.model_da = np.loadtxt(filepath_da, skiprows=2, dtype=self.dtype)
        self.model_db = np.loadtxt(filepath_db, skiprows=2, dtype=self.dtype)

        self.model_da["age"][self.model_da["age"] <= 1.0] += 1.0
        self.model_db["age"][self.model_db["age"] <= 1.0] += 1.0

    def list_atmosphere_parameters(self):
        """
        Print the formatted list of parameters available from the atmophere
        models.

        """

        for i, j in zip(self.column_names.items(), self.column_units.items()):
            print(f"Parameter: {i[1]}, Column Name: {i[0]}, Unit: {j[1]}")

    def interp_am(
        self,
        dependent="G3",
        atmosphere="H",
        independent=["logg", "Mbol"],
        logg=8.0,
        interpolator="CT",
        kwargs_for_RBF={},
        kwargs_for_CT={},
        allow_extrapolation=False,
    ):
        """
        This function interpolates the grid of synthetic photometry and a few
        other physical properties as a function of 2 independent variables,
        the Default choices are 'logg' and 'Mbol'.

        Parameters
        ----------
        dependent: str (Default: 'G3')
            The value to be interpolated over. Choose from:
            'Teff', 'logg', 'mass', 'Mbol', 'BC', 'U', 'B', 'V', 'R', 'I', 'J', 'H', 'Ks', 'Y_mko', 'J_mko', 'H_mko',
            'K_mko', 'W1', 'W2', 'W3', 'W4', 'S36', 'S45', 'S58', 'S80', 'u_sdss', 'g_sdss', 'r_sdss', 'i_sdss',
            'z_sdss', 'g_ps1', 'r_ps1', 'i_ps1', 'z_ps1', 'y_ps1', 'G2', 'G2_BP', 'G2_RP', 'G3', 'G3_BP', 'G3_RP',
            'FUV', 'NUV', 'age'.
        atmosphere: str (Default: 'H')
            The atmosphere type, 'H' or 'He'.
        independent: list (Default: ['logg', 'Mbol'])
            The parameters to be interpolated over for dependent.
        logg: float (Default: 8.0)
            Only used if independent is of length 1.
        interpolator: str (Default: 'RBF')
            Choose between 'RBF' and 'CT'.
        kwargs_for_RBF: dict (Default: {"neighbors": None, "smoothing": 0.0, "kernel": "thin_plate_spline",
            "epsilon": None, "degree": None,})
            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`.

        Returns
        -------
            A callable function of CloughTocher2DInterpolator.

        """
        _kwargs_for_RBF = {
            "neighbors": None,
            "smoothing": 0.0,
            "kernel": "thin_plate_spline",
            "epsilon": None,
            "degree": None,
        }
        _kwargs_for_RBF.update(**kwargs_for_RBF)

        _kwargs_for_CT = {
            "fill_value": np.inf,
            "tol": 1e-10,
            "maxiter": 100000,
            "rescale": True,
        }
        _kwargs_for_CT.update(**kwargs_for_CT)

        # DA atmosphere
        if atmosphere.lower() in ["h", "hydrogen", "da"]:
            model = self.model_da

        # DB atmosphere
        elif atmosphere.lower() in ["he", "helium", "db"]:
            model = self.model_db

        else:
            raise ValueError(
                'Please choose from "h", "hydrogen", "da", "he", "helium" or "db" as the atmophere type, you have '
                "provided {}.format(atmosphere.lower())"
            )

        independent = np.asarray(independent, dtype=object).reshape(-1)

        independent_list = self.column_key
        independent_list_lower_cases = np.char.lower(independent_list)

        # If only performing a 1D interpolation, the logg has to be assumed.
        if len(independent) == 1:
            if independent[0].lower() in independent_list_lower_cases:
                independent = np.array(("logg", independent[0]))

            else:
                raise ValueError(
                    "When ony interpolating in 1-dimension, the independent variable has to be one of: Teff, mass, "
                    "Mbol, or age."
                )

            _independent_arg_0 = np.where(independent[0].lower() == independent_list_lower_cases)[0][0]
            _independent_arg_1 = np.where(independent[1].lower() == independent_list_lower_cases)[0][0]

            independent = np.array([independent_list[_independent_arg_0], independent_list[_independent_arg_1]])

            arg_0 = model[independent[0]]
            arg_1 = model[independent[1]]

            arg_1_min = np.nanmin(arg_1)
            arg_1_max = np.nanmax(arg_1)

            if independent[1] in ["Teff", "age"]:
                arg_1 = np.log10(arg_1)

            if interpolator.lower() == "ct":
                # Interpolate with the scipy CloughTocher2DInterpolator
                _atmosphere_interpolator = CloughTocher2DInterpolator(
                    (arg_0, arg_1),
                    model[dependent],
                    **_kwargs_for_CT,
                )

                def atmosphere_interpolator(_x):
                    if independent[1] in ["Teff", "age"]:
                        _x = np.log10(_x)

                    return _atmosphere_interpolator(logg, _x)

            elif interpolator.lower() == "rbf":
                # Interpolate with the scipy RBFInterpolator
                _atmosphere_interpolator = RBFInterpolator(
                    np.stack((arg_0, arg_1), -1),
                    model[dependent],
                    **_kwargs_for_RBF,
                )

                def atmosphere_interpolator(_x):
                    _x_arr = np.asarray(_x).reshape(-1).astype(float)
                    length = _x_arr.size
                    _logg = np.full(length, logg, dtype=float)

                    if not allow_extrapolation:
                        _x_arr[_x_arr < arg_1_min] = arg_1_min
                        _x_arr[_x_arr > arg_1_max] = arg_1_max

                    if independent[1] in ["Teff", "age"]:
                        _x_arr = np.log10(_x_arr)

                    return _atmosphere_interpolator(np.column_stack((_logg, _x_arr)))

            else:
                raise ValueError("Interpolator should be CT or RBF, {interpolator} is given.")

        # If a 2D grid is to be interpolated, normally is the logg and another
        # parameter
        elif len(independent) == 2:
            _independent_arg_0 = np.where(independent[0].lower() == independent_list_lower_cases)[0][0]
            _independent_arg_1 = np.where(independent[1].lower() == independent_list_lower_cases)[0][0]

            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_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)

            if independent[0] in ["Teff", "age"]:
                arg_0 = np.log10(arg_0)

            if independent[1] in ["Teff", "age"]:
                arg_1 = np.log10(arg_1)

            if interpolator.lower() == "ct":
                # Interpolate with the scipy CloughTocher2DInterpolator
                _atmosphere_interpolator = CloughTocher2DInterpolator(
                    (arg_0, arg_1),
                    model[dependent],
                    **_kwargs_for_CT,
                )

                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)

                    out = _atmosphere_interpolator(_x_0, _x_1)
                    out = np.asarray(out).reshape(-1)
                    if not allow_extrapolation and np.any(mask_oob):
                        out[mask_oob] = -np.inf
                    return out

            elif interpolator.lower() == "rbf":
                # Interpolate with the scipy RBFInterpolator
                _atmosphere_interpolator = RBFInterpolator(
                    np.stack((arg_0, arg_1), -1),
                    model[dependent],
                    **_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)

                    out = _atmosphere_interpolator(np.column_stack((_x_0, _x_1)))
                    out = np.asarray(out).reshape(-1)
                    if np.any(mask_oob):
                        out[mask_oob] = -np.inf
                    return out

            else:
                raise ValueError("This should never happen.")

        else:
            raise TypeError("Please provide ONE varaible name as a string or list, or TWO varaible names in a list.")

        return atmosphere_interpolator

1	#!/usr/bin/env python3
2	# -- coding: utf-8 --
3
4	"""Handling the formatting of different atmosphere models"""
5
6	import os	2✔
7
8	import numpy as np	2✔
9	from scipy.interpolate import CloughTocher2DInterpolator	2✔
10	from scipy.interpolate import RBFInterpolator	2✔
11
12
13	class AtmosphereModelReader(object):	2✔
14	"""Handling the formatting of different atmosphere models"""
15
16	def __init__(self):	2✔
17	super(AtmosphereModelReader, self).__init__()	2✔
18
19	self.this_file = os.path.dirname(os.path.abspath(__file__))	2✔
20
21	self.model_list = {	2✔
22	"montreal_co_da_20": "Bedard et al. 2020 CO DA",
23	"montreal_co_db_20": "Bedard et al. 2020 CO DB",
24	"lpcode_he_da_07": "Panei et al. 2007 He DA",
25	"lpcode_co_da_07": "Panei et al. 2007 CO DA",
26	"lpcode_he_da_09": "Althaus et al. 2009 He DA",
27	"lpcode_co_da_10_z001": "Renedo et al. 2010 CO DA Z=0.01",
28	"lpcode_co_da_10_z0001": "Renedo et al. 2010 CO DA Z=0.001",
29	"lpcode_co_da_15_z00003": "Althaus et al. 2015 DA Z=0.00003",
30	"lpcode_co_da_15_z0001": "Althaus et al. 2015 DA Z=0.0001",
31	"lpcode_co_da_15_z0005": "Althaus et al. 2015 DA Z=0.0005",
32	"lpcode_co_db_17_z00005": "Althaus et al. 2017 DB Y=0.4",
33	"lpcode_co_db_17_z0001": "Althaus et al. 2017 DB Y=0.4",
34	"lpcode_co_db_17": "Camisassa et al. 2017 DB",
35	"lpcode_one_da_07": "Althaus et al. 2007 ONe DA",
36	"lpcode_one_da_19": "Camisassa et al. 2019 ONe DA",
37	"lpcode_one_db_19": "Camisassa et al. 2019 ONe DB",
38	"lpcode_da_22": "Althaus et al. 2013 He DA, Camisassa et al. 2016 CO DA, Camisassa et al. 2019 ONe DA",
39	"lpcode_db_22": "Camisassa et al. 2017 CO DB, " + "Camisassa et al. 2019 ONe DB",
40	}
41
42	# DA atmosphere
43	filepath_da = os.path.join(	2✔
44	os.path.dirname(os.path.abspath(__file__)),
45	"wd_photometry/Table_DA_13012021.txt",
46	)
47
48	# DB atmosphere
49	filepath_db = os.path.join(	2✔
50	os.path.dirname(os.path.abspath(__file__)),
51	"wd_photometry/Table_DB_13012021.txt",
52	)
53
54	# Prepare the array column dtype
55	self.column_key = np.array(	2✔
56	[
57	"Teff",
58	"logg",
59	"mass",
60	"Mbol",
61	"BC",
62	"U",
63	"B",
64	"V",
65	"R",
66	"I",
67	"J",
68	"H",
69	"Ks",
70	"Y_mko",
71	"J_mko",
72	"H_mko",
73	"K_mko",
74	"W1",
75	"W2",
76	"W3",
77	"W4",
78	"S36",
79	"S45",
80	"S58",
81	"S80",
82	"u_sdss",
83	"g_sdss",
84	"r_sdss",
85	"i_sdss",
86	"z_sdss",
87	"g_ps1",
88	"r_ps1",
89	"i_ps1",
90	"z_ps1",
91	"y_ps1",
92	"G2",
93	"G2_BP",
94	"G2_RP",
95	"G3",
96	"G3_BP",
97	"G3_RP",
98	"FUV",
99	"NUV",
100	"age",
101	]
102	)
103	self.column_key_formatted = np.array(	2✔
104	[
105	r"T$_{\mathrm{eff}}$",
106	"log(g)",
107	"Mass",
108	r"M$_{\mathrm{bol}}$",
109	"BC",
110	r"$U$",
111	r"$B$",
112	r"$V$",
113	r"$R$",
114	r"$I$",
115	r"$J$",
116	r"$H$",
117	r"$K_{\mathrm{s}}$",
118	r"$Y_{\mathrm{MKO}}$",
119	r"$J_{\mathrm{MKO}}$",
120	r"$H_{\mathrm{MKO}}$",
121	r"$K_{\mathrm{MKO}}$",
122	r"$W_{1}$",
123	r"$W_{2}$",
124	r"$W_{3}$",
125	r"$W_{4}$",
126	r"$S_{36}$",
127	r"$S_{45}$",
128	r"$S_{58}$",
129	r"$S_{80}$",
130	r"u$_{\mathrm{SDSS}}$",
131	r"$g_{\mathrm{SDSS}}$",
132	r"$r_{\mathrm{SDSS}}$",
133	r"$i_{\mathrm{SDSS}}$",
134	r"$z_{\mathrm{SDSS}}$",
135	r"$g_{\mathrm{PS1}}$",
136	r"$r_{\mathrm{PS1}}$",
137	r"$i_{\mathrm{PS1}}$",
138	r"$z_{\mathrm{PS1}}$",
139	r"$y_{\mathrm{PS1}}$",
140	r"$G_{\mathrm{DR2}}$",
141	r"$G_{\mathrm{BP, DR2}}$",
142	r"$G_{\mathrm{RP, DR2}}$",
143	r"$G_{\mathrm{DR3}}$",
144	r"$G_{\mathrm{BP, DR3}}$",
145	r"$G_{\mathrm{RP, DR3}}$",
146	"FUV",
147	"NUV",
148	"Age",
149	]
150	)
151	self.column_key_unit = np.array(	2✔
152	[
153	"K",
154	r"(cm/s$^2$)",
155	r"M$_\odot$",
156	"mag",
157	"mag",
158	"mag",
159	"mag",
160	"mag",
161	"mag",
162	"mag",
163	"mag",
164	"mag",
165	"mag",
166	"mag",
167	"mag",
168	"mag",
169	"mag",
170	"mag",
171	"mag",
172	"mag",
173	"mag",
174	"mag",
175	"mag",
176	"mag",
177	"mag",
178	"mag",
179	"mag",
180	"mag",
181	"mag",
182	"mag",
183	"mag",
184	"mag",
185	"mag",
186	"mag",
187	"mag",
188	"mag",
189	"mag",
190	"mag",
191	"mag",
192	"mag",
193	"mag",
194	"mag",
195	"mag",
196	"yr",
197	]
198	)
199	self.column_key_wavelength = np.array(	2✔
200	[
201	0.0,
202	0.0,
203	0.0,
204	0.0,
205	0.0,
206	3585.0,
207	4371.0,
208	5478.0,
209	6504.0,
210	8020.0,
211	12350.0,
212	16460.0,
213	21600.0,
214	10310.0,
215	12500.0,
216	16360.0,
217	22060.0,
218	33682.0,
219	46179.0,
220	120717.0,
221	221944.0,
222	35378.0,
223	44780.0,
224	56962.0,
225	77978.0,
226	3557.0,
227	4702.0,
228	6175.0,
229	7491.0,
230	8946.0,
231	4849.0,
232	6201.0,
233	7535.0,
234	8674.0,
235	9628.0,
236	6229.0,
237	5037.0,
238	7752.0,
239	6218.0,
240	5110.0,
241	7769.0,
242	1535.0,
243	2301.0,
244	0.0,
245	]
246	)
247
248	self.column_names = {}	2✔
249	self.column_units = {}	2✔
250	self.column_wavelengths = {}	2✔
251	for i, j, k, _l in zip(	2✔
252	self.column_key,
253	self.column_key_formatted,
254	self.column_key_unit,
255	self.column_key_wavelength,
256	):
257	self.column_names[i] = j	2✔
258	self.column_units[i] = k	2✔
259	self.column_wavelengths[i] = _l	2✔
260
261	self.column_type = np.array(([np.float64] * len(self.column_key)))	2✔
262	self.dtype = list(zip(self.column_key, self.column_type))	2✔
263
264	# Load the synthetic photometry file in a recarray
265	self.model_da = np.loadtxt(filepath_da, skiprows=2, dtype=self.dtype)	2✔
266	self.model_db = np.loadtxt(filepath_db, skiprows=2, dtype=self.dtype)	2✔
267
268	self.model_da["age"][self.model_da["age"] <= 1.0] += 1.0	2✔
269	self.model_db["age"][self.model_db["age"] <= 1.0] += 1.0	2✔
270
271	def list_atmosphere_parameters(self):	2✔
272	"""
273	Print the formatted list of parameters available from the atmophere
274	models.
275
276	"""
277
278	for i, j in zip(self.column_names.items(), self.column_units.items()):	2✔
279	print(f"Parameter: {i[1]}, Column Name: {i[0]}, Unit: {j[1]}")	2✔
280
281	def interp_am(	2✔
282	self,
283	dependent="G3",
284	atmosphere="H",
285	independent=["logg", "Mbol"],
286	logg=8.0,
287	interpolator="CT",
288	kwargs_for_RBF={},
289	kwargs_for_CT={},
290	allow_extrapolation=False,
291	):
292	"""
293	This function interpolates the grid of synthetic photometry and a few
294	other physical properties as a function of 2 independent variables,
295	the Default choices are 'logg' and 'Mbol'.
296
297	Parameters
298	----------
299	dependent: str (Default: 'G3')
300	The value to be interpolated over. Choose from:
301	'Teff', 'logg', 'mass', 'Mbol', 'BC', 'U', 'B', 'V', 'R', 'I', 'J', 'H', 'Ks', 'Y_mko', 'J_mko', 'H_mko',
302	'K_mko', 'W1', 'W2', 'W3', 'W4', 'S36', 'S45', 'S58', 'S80', 'u_sdss', 'g_sdss', 'r_sdss', 'i_sdss',
303	'z_sdss', 'g_ps1', 'r_ps1', 'i_ps1', 'z_ps1', 'y_ps1', 'G2', 'G2_BP', 'G2_RP', 'G3', 'G3_BP', 'G3_RP',
304	'FUV', 'NUV', 'age'.
305	atmosphere: str (Default: 'H')
306	The atmosphere type, 'H' or 'He'.
307	independent: list (Default: ['logg', 'Mbol'])
308	The parameters to be interpolated over for dependent.
309	logg: float (Default: 8.0)
310	Only used if independent is of length 1.
311	interpolator: str (Default: 'RBF')
312	Choose between 'RBF' and 'CT'.
313	kwargs_for_RBF: dict (Default: {"neighbors": None, "smoothing": 0.0, "kernel": "thin_plate_spline",
314	"epsilon": None, "degree": None,})
315	Keyword argument for the interpolator. See `scipy.interpolate.RBFInterpolator`.
316	kwargs_for_CT: dict (Default: {'fill_value': -np.inf, 'tol': 1e-10, 'maxiter': 100000})
317	Keyword argument for the interpolator. See `scipy.interpolate.CloughTocher2DInterpolator`.
318
319	Returns
320	-------
321	A callable function of CloughTocher2DInterpolator.
322
323	"""
324	_kwargs_for_RBF = {	2✔
325	"neighbors": None,
326	"smoothing": 0.0,
327	"kernel": "thin_plate_spline",
328	"epsilon": None,
329	"degree": None,
330	}
331	_kwargs_for_RBF.update(**kwargs_for_RBF)	2✔
332
333	_kwargs_for_CT = {	2✔
334	"fill_value": np.inf,
335	"tol": 1e-10,
336	"maxiter": 100000,
337	"rescale": True,
338	}
339	_kwargs_for_CT.update(**kwargs_for_CT)	2✔
340
341	# DA atmosphere
342	if atmosphere.lower() in ["h", "hydrogen", "da"]:	2✔
343	model = self.model_da	2✔
344
345	# DB atmosphere
346	elif atmosphere.lower() in ["he", "helium", "db"]:	2✔
347	model = self.model_db	2✔
348
349	else:
350	raise ValueError(	2✔
351	'Please choose from "h", "hydrogen", "da", "he", "helium" or "db" as the atmophere type, you have '
352	"provided {}.format(atmosphere.lower())"
353	)
354
355	independent = np.asarray(independent, dtype=object).reshape(-1)	2✔
356
357	independent_list = self.column_key	2✔
358	independent_list_lower_cases = np.char.lower(independent_list)	2✔
359
360	# If only performing a 1D interpolation, the logg has to be assumed.
361	if len(independent) == 1:	2✔
362	if independent[0].lower() in independent_list_lower_cases:	2✔
363	independent = np.array(("logg", independent[0]))	2✔
364
365	else:
366	raise ValueError(	×
367	"When ony interpolating in 1-dimension, the independent variable has to be one of: Teff, mass, "
368	"Mbol, or age."
369	)
370
371	_independent_arg_0 = np.where(independent[0].lower() == independent_list_lower_cases)[0][0]	2✔
372	_independent_arg_1 = np.where(independent[1].lower() == independent_list_lower_cases)[0][0]	2✔
373
374	independent = np.array([independent_list[_independent_arg_0], independent_list[_independent_arg_1]])	2✔
375
376	arg_0 = model[independent[0]]	2✔
377	arg_1 = model[independent[1]]	2✔
378
379	arg_1_min = np.nanmin(arg_1)	2✔
380	arg_1_max = np.nanmax(arg_1)	2✔
381
382	if independent[1] in ["Teff", "age"]:	2✔
383	arg_1 = np.log10(arg_1)	2✔
384
385	if interpolator.lower() == "ct":	2✔
386	# Interpolate with the scipy CloughTocher2DInterpolator
387	_atmosphere_interpolator = CloughTocher2DInterpolator(	2✔
388	(arg_0, arg_1),
389	model[dependent],
390	**_kwargs_for_CT,
391	)
392
393	def atmosphere_interpolator(_x):	2✔
394	if independent[1] in ["Teff", "age"]:	2✔
395	_x = np.log10(_x)	2✔
396
397	return _atmosphere_interpolator(logg, _x)	2✔
398
399	elif interpolator.lower() == "rbf":	2✔
400	# Interpolate with the scipy RBFInterpolator
401	_atmosphere_interpolator = RBFInterpolator(	2✔
402	np.stack((arg_0, arg_1), -1),
403	model[dependent],
404	**_kwargs_for_RBF,
405	)
406
407	def atmosphere_interpolator(_x):	2✔
408	_x_arr = np.asarray(_x).reshape(-1).astype(float)	2✔
409	length = _x_arr.size	2✔
410	_logg = np.full(length, logg, dtype=float)	2✔
411
412	if not allow_extrapolation:	2✔
413	_x_arr[_x_arr < arg_1_min] = arg_1_min	2✔
414	_x_arr[_x_arr > arg_1_max] = arg_1_max	2✔
415
416	if independent[1] in ["Teff", "age"]:	2✔
417	_x_arr = np.log10(_x_arr)	2✔
418
419	return _atmosphere_interpolator(np.column_stack((_logg, _x_arr)))	2✔
420
421	else:
422	raise ValueError("Interpolator should be CT or RBF, {interpolator} is given.")	2✔
423
424	# If a 2D grid is to be interpolated, normally is the logg and another
425	# parameter
426	elif len(independent) == 2:	2✔
427	_independent_arg_0 = np.where(independent[0].lower() == independent_list_lower_cases)[0][0]	2✔
428	_independent_arg_1 = np.where(independent[1].lower() == independent_list_lower_cases)[0][0]	2✔
429
430	independent = np.array([independent_list[_independent_arg_0], independent_list[_independent_arg_1]])	2✔
431
432	arg_0 = model[independent[0]]	2✔
433	arg_1 = model[independent[1]]	2✔
434
435	arg_0_min = np.nanmin(arg_0)	2✔
436	arg_0_max = np.nanmax(arg_0)	2✔
437	arg_1_min = np.nanmin(arg_1)	2✔
438	arg_1_max = np.nanmax(arg_1)	2✔
439
440	if independent[0] in ["Teff", "age"]:	2✔
441	arg_0 = np.log10(arg_0)	2✔
442
443	if independent[1] in ["Teff", "age"]:	2✔
444	arg_1 = np.log10(arg_1)	2✔
445
446	if interpolator.lower() == "ct":	2✔
447	# Interpolate with the scipy CloughTocher2DInterpolator
448	_atmosphere_interpolator = CloughTocher2DInterpolator(	2✔
449	(arg_0, arg_1),
450	model[dependent],
451	**_kwargs_for_CT,
452	)
453
454	def atmosphere_interpolator(x0, x1=None):	2✔
455	# Support scalar/array inputs for both coordinates, with simple broadcasting
456	if x1 is None:	2✔
457	arr = np.asarray(x0).reshape(-1)	2✔
458	if arr.size >= 2:	2✔
459	x_0, x_1 = arr[0], arr[1]	2✔
NEW 460	elif arr.size == 1:	×
NEW 461	x_0, x_1 = arr[0], arr[0]	×
462	else:
NEW 463	x_0, x_1 = -np.inf, -np.inf	×
464	else:
465	x_0, x_1 = x0, x1	2✔
466
467	if isinstance(x_0, (float, int, np.integer)):	2✔
468	length0 = 1	2✔
469	else:
470	length0 = np.asarray(x_0).size	2✔
471
472	if isinstance(x_1, (float, int, np.integer)):	2✔
473	length1 = 1	2✔
474	else:
475	length1 = np.asarray(x_1).size	2✔
476
477	if length0 == length1:	2✔
478	pass	2✔
479	elif (length0 == 1) and (length1 > 1):	2✔
480	x_0 = [x_0] * length1	2✔
481	length0 = length1	2✔
NEW 482	elif (length0 > 1) and (length1 == 1):	×
NEW 483	x_1 = [x_1] * length0	×
NEW 484	length1 = length0	×
485	else:
NEW 486	raise ValueError(	×
487	"Either one variable is a float, int or of size 1, or two variables should have the same"
488	"size."
489	)
490
491	_x_0 = np.asarray(x_0).reshape(-1).astype(float)	2✔
492	_x_1 = np.asarray(x_1).reshape(-1).astype(float)	2✔
493
494	# mark out-of-range inputs to avoid excessive extrapolation
495	mask_oob = (_x_0 < arg_0_min) \| (_x_0 > arg_0_max) \| (_x_1 < arg_1_min) \| (_x_1 > arg_1_max)	2✔
496
497	if independent[0] in ["Teff", "age"]:	2✔
498	_x_0 = np.log10(_x_0)	2✔
499
500	if independent[1] in ["Teff", "age"]:	2✔
501	_x_1 = np.log10(_x_1)	2✔
502
503	out = _atmosphere_interpolator(_x_0, _x_1)	2✔
504	out = np.asarray(out).reshape(-1)	2✔
505	if not allow_extrapolation and np.any(mask_oob):	2✔
506	out[mask_oob] = -np.inf	2✔
507	return out	2✔
508
509	elif interpolator.lower() == "rbf":	2✔
510	# Interpolate with the scipy RBFInterpolator
511	_atmosphere_interpolator = RBFInterpolator(	2✔
512	np.stack((arg_0, arg_1), -1),
513	model[dependent],
514	**_kwargs_for_RBF,
515	)
516
517	def atmosphere_interpolator(*x):	2✔
518	# Accept (x0, x1) or single array-like; use first two values, duplicate if only one
519	if len(x) == 2:	2✔
520	x_0, x_1 = x	2✔
521	elif len(x) == 1:	2✔
522	arr = np.asarray(x[0]).reshape(-1)	2✔
523	if arr.size >= 2:	2✔
524	x_0, x_1 = arr[0], arr[1]	2✔
NEW 525	elif arr.size == 1:	×
NEW 526	x_0, x_1 = arr[0], arr[0]	×
527	else:
NEW 528	x_0, x_1 = -np.inf, -np.inf	×
529	else:
NEW 530	x_0, x_1 = -np.inf, -np.inf	×
531
532	if isinstance(x_0, (float, int, np.integer)):	2✔
533	length0 = 1	2✔
534	else:
535	length0 = np.asarray(x_0).size	2✔
536
537	if isinstance(x_1, (float, int, np.integer)):	2✔
538	length1 = 1	2✔
539	else:
540	length1 = np.asarray(x_1).size	2✔
541
542	if length0 == length1:	2✔
543	pass	2✔
544
545	elif (length0 == 1) & (length1 > 1):	2✔
546	x_0 = [x_0] * length1	2✔
547	length0 = length1	2✔
548
549	elif (length0 > 1) & (length1 == 1):	×
550	x_1 = [x_1] * length0	×
551	length1 = length0	×
552
553	else:
554	raise ValueError(	×
555	"Either one variable is a float, int or of size 1, or two variables should have the same "
556	"size."
557	)
558
559	_x_0 = np.asarray(x_0).reshape(-1).astype(float)	2✔
560	_x_1 = np.asarray(x_1).reshape(-1).astype(float)	2✔
561
562	# mark out-of-range inputs to avoid excessive extrapolation
563	mask_oob = (_x_0 < arg_0_min) \| (_x_0 > arg_0_max) \| (_x_1 < arg_1_min) \| (_x_1 > arg_1_max)	2✔
564
565	if independent[0] in ["Teff", "age"]:	2✔
566	_x_0 = np.log10(_x_0)	2✔
567
568	if independent[1] in ["Teff", "age"]:	2✔
569	_x_1 = np.log10(_x_1)	2✔
570
571	out = _atmosphere_interpolator(np.column_stack((_x_0, _x_1)))	2✔
572	out = np.asarray(out).reshape(-1)	2✔
573	if np.any(mask_oob):	2✔
574	out[mask_oob] = -np.inf	2✔
575	return out	2✔
576
577	else:
578	raise ValueError("This should never happen.")	2✔
579
580	else:
581	raise TypeError("Please provide ONE varaible name as a string or list, or TWO varaible names in a list.")	2✔
582
583	return atmosphere_interpolator	2✔

cylammarco / WDPhotTools / 20978266610

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous