Skip to content

API Reference

This section provides a detailed reference for the code in the Star Log-extended eMulator package.

SLM

augment_data_multiple_columns(X)

Augment the data matrix X with nonlinear terms for multiple variables.

Parameters:

Name Type Description Default
X ndarray

The data matrix where each row is a variable, and each column is a snapshot in time.

 required

Returns:

Name Type Description
augmented_X ndarray

The augmented data matrix with quadratic and cross-product terms appended below the original rows.
Source code in src/slmemulator/SLM.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def augment_data_multiple_columns(X):
    r"""
    Augment the data matrix X with nonlinear terms for multiple variables.

    Parameters:
        X (np.ndarray): The data matrix where each row is a variable, and each
            column is a snapshot in time.

    Returns:
        augmented_X (np.ndarray): The augmented data matrix with quadratic and
            cross-product terms appended below the original rows.
    """
    n_variables, n_snapshots = X.shape

    pairs = list(combinations_with_replacement(range(n_variables), 2))
    augmented_X = np.empty((n_variables + len(pairs), n_snapshots), dtype=X.dtype)
    augmented_X[:n_variables, :] = X

    for row_idx, (i, j) in enumerate(pairs, start=n_variables):
        augmented_X[row_idx, :] = X[i, :] * X[j, :]

    return augmented_X

SLM(X, dt, error_threshold=0.0001, max_r=None)

Dynamic Mode Decomposition of the augmented data. Automatically determines the number of modes (r) based on an error threshold.

Parameters:

Name Type Description Default
X ndarray

The data matrix where each row is a variable, and each column is a snapshot in time. Expected to be log-transformed where appropriate.

 required
dt float

The time difference of linear DMDs.

 required
error_threshold float

(Optional) The maximum allowed absolute difference between the original data and the DMD reconstruction. Defaults to 1e-4.

 0.0001
max_r int

(Optional) The maximum number of modes to consider. If None, goes up to the maximum possible rank (min(X.shape)).

 None

Returns:

Name Type Description
tuple

(Phi, omega, lambda_vals, b, Xdmd, S, r_optimal) where Phi are

the DMD modes (truncated to the original variables), omega the

continuous-time eigenvalues, lambda_vals the discrete-time

eigenvalues, b the mode amplitudes, Xdmd the reconstruction of the

original (non-augmented) variables, S the singular values of the

augmented snapshot matrix, and r_optimal the selected rank.
Source code in src/slmemulator/SLM.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
def SLM(X, dt, error_threshold=1e-4, max_r=None):
    r"""
    Dynamic Mode Decomposition of the augmented data.
    Automatically determines the number of modes (r) based on an error threshold.

    Parameters:
        X (np.ndarray): The data matrix where each row is a variable, and each
            column is a snapshot in time. Expected to be log-transformed where
            appropriate.
        dt (float): The time difference of linear DMDs.
        error_threshold (float): (Optional) The maximum allowed absolute
            difference between the original data and the DMD reconstruction.
            Defaults to 1e-4.
        max_r (int): (Optional) The maximum number of modes to consider.
            If None, goes up to the maximum possible rank (min(X.shape)).

    Returns:
        tuple: (Phi, omega, lambda_vals, b, Xdmd, S, r_optimal) where Phi are
        the DMD modes (truncated to the original variables), omega the
        continuous-time eigenvalues, lambda_vals the discrete-time
        eigenvalues, b the mode amplitudes, Xdmd the reconstruction of the
        original (non-augmented) variables, S the singular values of the
        augmented snapshot matrix, and r_optimal the selected rank.
    """
    n = X.shape[0]  # Original number of variables before augmentation

    X_augmented = augment_data_multiple_columns(X)
    X1 = X_augmented[:, :-1]  # All columns except the last
    X2 = X_augmented[:, 1:]  # All columns except the first

    # Compute SVD of X1 once
    U_full, S_full, Vt_full = np.linalg.svd(X1, full_matrices=False)

    max_possible_r = min(X1.shape)
    max_r_to_check = max_possible_r if max_r is None else min(max_r, max_possible_r)

    # Reconstruction times (all snapshots)
    t = np.arange(X1.shape[1] + 1) * dt
    x1 = X1[:, 0]

    r_optimal = 1
    min_error = np.inf
    best = None

    # Find the smallest r whose reconstruction meets the threshold; otherwise
    # keep the r with the smallest error.
    for r_current in range(1, max_r_to_check + 1):
        U_r = U_full[:, :r_current]
        S_r_inv = np.diag(1.0 / S_full[:r_current])
        V_r = Vt_full[:r_current, :]

        Atilde = U_r.T @ X2 @ V_r.T @ S_r_inv
        lambda_vals, W_r = np.linalg.eig(Atilde)  # discrete-time eigenvalues

        Phi = X2 @ V_r.T @ S_r_inv @ W_r  # DMD modes
        omega = np.log(lambda_vals) / dt  # continuous-time eigenvalues

        # DMD mode amplitudes and reconstruction
        b = np.linalg.lstsq(Phi, x1, rcond=None)[0]
        time_dynamics = b[:, np.newaxis] * np.exp(omega[:, np.newaxis] * t)
        Xdmd = (Phi @ time_dynamics)[:n, :]  # truncate to original variables

        current_error = np.max(np.abs(X - Xdmd))

        if current_error <= error_threshold or current_error < min_error:
            r_optimal = r_current
            min_error = current_error
            best = (Phi[:n, :], omega, lambda_vals, b, Xdmd)

        if current_error <= error_threshold:
            break

    print(f"Optimal 'r' determined: {r_optimal} (Max absolute error = {min_error:.6f})")

    best_Phi, best_omega, best_lambda_vals, best_b, best_Xdmd = best
    return best_Phi, best_omega, best_lambda_vals, best_b, best_Xdmd, S_full, r_optimal

solve_tov(fileName, tidal=False, parametric=False, mseos=True)

Solves the TOV equation and returns radius, central pressure and mass.

Parameters:

Name Type Description Default
fileName str

Filename containing the EOS in the format nb (fm^-3), E (MeV), P (MeV/fm^3). For non-parametric runs this is the name of a file in the packaged EOS_Data; for parametric runs it is looked up in the generated EOS_files directory.

 required
tidal bool

Also compute the tidal Love number k2. Default False.

 False
parametric bool

Whether the EOS file comes from a parametric run.

 False
mseos bool

For parametric runs, MSEOS (True) or Quarkyonia (False).

 True

Returns:

Name Type Description
dataArray ndarray

Data array containing radii, central pressure and mass (includes tidal deformability k_2 if tidal is True).
Source code in src/slmemulator/SLM.py 
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def solve_tov(fileName, tidal=False, parametric=False, mseos=True):
    r"""
    Solves the TOV equation and returns radius, central pressure and mass.

    Parameters:
        fileName (str): Filename containing the EOS in the format nb (fm^-3),
            E (MeV), P (MeV/fm^3). For non-parametric runs this is the name of
            a file in the packaged EOS_Data; for parametric runs it is looked
            up in the generated EOS_files directory.
        tidal (bool): Also compute the tidal Love number k2. Default False.
        parametric (bool): Whether the EOS file comes from a parametric run.
        mseos (bool): For parametric runs, MSEOS (True) or Quarkyonia (False).

    Returns:
        dataArray (np.ndarray): Data array containing radii, central pressure
            and mass (includes tidal deformability k_2 if tidal is True).
    """
    paths = get_paths()

    if not parametric:
        # EOS file shipped with the package
        eos_file_path = paths["package_eos_data_dir"] / fileName
        if not eos_file_path.is_file():
            raise FileNotFoundError(
                f"Internal EOS file '{fileName}' not found in package data."
            )
        tov_path_target = paths["user_tov_data_dir"]
    elif mseos:
        eos_file_path = paths["mseos_path_specific"] / fileName
        tov_path_target = paths["mseos_tov_path_specific"]
    else:
        eos_file_path = paths["qeos_path_specific"] / fileName
        tov_path_target = paths["qeos_tov_path_specific"]

    tov = TOV(str(eos_file_path), tidal=tidal)
    tov_path_target.mkdir(parents=True, exist_ok=True)

    tov.tov_routine(verbose=False, write_to_file=False)
    print("R of 1.4 solar mass star: ", tov.canonical_NS_radius())

    dataArray = [
        tov.total_radius.flatten(),
        tov.total_pres_central.flatten(),
        tov.total_mass.flatten(),
    ]
    if tidal:
        dataArray.append(tov.k2.flatten())
    dataArray = np.asarray(dataArray, dtype=np.float64)

    # Output name: MR_<params>.txt for parametric files, MR_<name>_TOV.txt otherwise
    name_parts = Path(fileName).name.removesuffix(".txt").split("_")
    if len(name_parts) > 2:
        output_file_name = "MR_" + "_".join(name_parts[1:]) + ".txt"
    else:
        output_file_name = "_".join(["MR", name_parts[0], "TOV"]) + ".txt"

    np.savetxt(tov_path_target / output_file_name, dataArray.T, fmt="%1.8e")
    return dataArray

cleanData

clean_directory(directory: str | None = None) -> None

Recursively cleans a specified directory by removing common project artifacts and specific, code-generated subdirectories.

The function targets temporary files (by extension) and removes specific directories generated during modeling, plotting, and data processing.

Parameters:

Name Type Description Default
directory str

The path to the directory to clean. If :obj:None, the function The path to the directory to clean. If :obj:None, the function

None
defaults to cleaning the **current working directory** (

func:os.getcwd).

 required

Returns:

Name Type Description
None None

The function modifies the filesystem but does not return a value.
Source code in src/slmemulator/cleanData.py 
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def clean_directory(directory: str | None = None) -> None:
    """
    Recursively cleans a specified directory by removing common project artifacts
    and specific, code-generated subdirectories.

    The function targets temporary files (by extension) and removes specific
    directories generated during modeling, plotting, and data processing.

    Parameters:
        directory (str, optional): The path to the directory to clean. If :obj:`None`, the function         The path to the directory to clean. If :obj:`None`, the function 
        defaults to cleaning the **current working directory** (:func:`os.getcwd`).

    Returns:
        None: The function modifies the filesystem but does not return a value.

    """
    if directory is None:
        directory = os.getcwd()

    # Ensure the target directory exists and is a directory
    target_dir_path = Path(directory).resolve()
    if not target_dir_path.is_dir():
        print(f"Error: Directory '{directory}' not found or is not a directory.")
        return

    print(f"Cleaning directory: {target_dir_path}")

    # List of additional folders created by the code that you want to remove recursively.
    # These are now relative to the 'directory' argument.
    additional_folders_to_clean_names = [
        DEFAULT_EOS_FILES_SUBDIR_NAME,
        DEFAULT_RESULTS_SUBDIR_NAME,
        DEFAULT_TOV_DATA_SUBDIR_NAME,
        DEFAULT_TEST_DATA_SUBDIR_NAME,
        DEFAULT_TRAIN_DATA_SUBDIR_NAME,
        DEFAULT_VAL_DATA_SUBDIR_NAME,
        DEFAULT_PLOTS_SUBDIR_NAME,
    ]

    # Convert these names to full paths within the target directory
    additional_folders_to_clean_paths = [
        target_dir_path / name for name in additional_folders_to_clean_names
    ]

    for root, dirs, files in os.walk(
        target_dir_path, topdown=False
    ):  # Traverse from bottom to top
        # Clean files matching patterns in cleanup_targets
        for file in files:
            for target_ext in [
                t for t in cleanup_targets if t.startswith(".")
            ]:  # Only check extensions
                if file.endswith(target_ext):
                    file_path = os.path.join(root, file)
                    print(f"Removing file: {file_path}")
                    os.remove(file_path)

        # Clean directories matching patterns in cleanup_targets (by name)
        # or directories matching full paths in additional_folders_to_clean_paths
        for dir_name in dirs:
            dir_full_path = Path(root) / dir_name  # Use Path for comparison
            if (
                dir_name in cleanup_targets
                or dir_full_path in additional_folders_to_clean_paths
            ):
                print(f"Removing directory: {dir_full_path}")
                shutil.rmtree(dir_full_path)

config

Project path configuration for slmemulator.

get_paths(output_base_dir: Path | None = None, eos_name: str = 'MSEOS', is_parametric_run: bool = True, include_slm_paths: bool = True) -> dict[str, Path]

Generates and returns a dictionary of resolved project paths, dynamically structuring subdirectories based on the Equation of State (EOS) name and run configuration.

The function provides paths for input data, model binaries, general output, and specific subdirectories for results, plots, and test data related to SLM (Star Log-extended eMulator) or pSLM (parametric SLM) runs.

Parameters:

Name Type Description Default
output_base_dir Path

The root directory where all generated project outputs (results, plots, test data) will be stored. If None, the repository root is used. Defaults to None.

 None
eos_name str

The name of the Equation of State (e.g., "MSEOS", "QEOS", "APR"). This name dictates the specific subdirectory created for the current run within the results, plots, and test directories. Defaults to "MSEOS".

 'MSEOS'
is_parametric_run bool

Flag indicating if the current modeling run is using the parametric SLM (pSLM) approach. If True, the output paths will include 'pSLM'; otherwise 'SLM'. Defaults to True.

 True
include_slm_paths bool

If True, the dictionary includes the run-specific directories (current_slm_results_dir, current_slm_plots_dir, current_slm_tests_dir). Defaults to True.

 True

Returns:

Type Description
dict[str, Path]

dict[str, pathlib.Path]: A dictionary containing all relevant path
dict[str, Path]

configurations.
Source code in src/slmemulator/config.py 
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def get_paths(
    output_base_dir: Path | None = None,
    eos_name: str = "MSEOS",
    is_parametric_run: bool = True,
    include_slm_paths: bool = True,
) -> dict[str, Path]:
    """
    Generates and returns a dictionary of resolved project paths, dynamically
    structuring subdirectories based on the Equation of State (EOS) name and
    run configuration.

    The function provides paths for input data, model binaries, general output,
    and specific subdirectories for results, plots, and test data related to
    SLM (Star Log-extended eMulator) or pSLM (parametric SLM) runs.

    Parameters:
        output_base_dir (pathlib.Path, optional):
            The root directory where all generated project outputs (results,
            plots, test data) will be stored. If ``None``, the repository root
            is used. Defaults to ``None``.
        eos_name (str, optional):
            The name of the Equation of State (e.g., "MSEOS", "QEOS", "APR").
            This name dictates the specific subdirectory created for the
            current run within the results, plots, and test directories.
            Defaults to "MSEOS".
        is_parametric_run (bool, optional):
            Flag indicating if the current modeling run is using the parametric
            SLM (pSLM) approach. If ``True``, the output paths will include
            'pSLM'; otherwise 'SLM'. Defaults to ``True``.
        include_slm_paths (bool, optional):
            If ``True``, the dictionary includes the run-specific directories
            (``current_slm_results_dir``, ``current_slm_plots_dir``,
            ``current_slm_tests_dir``). Defaults to ``True``.

    Returns:
        dict[str, pathlib.Path]: A dictionary containing all relevant path
        configurations.
    """
    output_data_base = output_base_dir or PROJECT_ROOT
    src_dir = PROJECT_ROOT / DEFAULT_SRC_SUBDIR_NAME

    # Strip a known extension and upper-case the EOS name for directory naming
    clean_eos_name = eos_name
    for ext in _EXTENSIONS_TO_STRIP:
        if clean_eos_name.lower().endswith(ext):
            clean_eos_name = clean_eos_name[: -len(ext)]
            break
    eos_folder_name = clean_eos_name.upper()

    slm_subdir = (
        DEFAULT_PSLM_SUBDIR_NAME if is_parametric_run else DEFAULT_SLM_SUBDIR_NAME
    )

    # Concrete Path is fine here: the package is installed on the filesystem
    # (packaged EOS data could not be read through np.loadtxt otherwise).
    package_root = Path(str(resources.files("slmemulator")))

    paths: dict[str, Path] = {
        "project_root": PROJECT_ROOT,
        "src_dir": src_dir,
        # Internal package resources (read-only)
        "package_eos_codes_dir": package_root.joinpath(DEFAULT_EOS_CODES_SUBDIR_NAME),
        "package_eos_data_dir": package_root.joinpath(DEFAULT_EOS_DATA_SUBDIR_NAME),
        # General output directories
        "plots_dir": output_data_base / DEFAULT_PLOTS_SUBDIR_NAME,
        "results_dir": output_data_base / DEFAULT_RESULTS_SUBDIR_NAME,
        "docs_dir": output_data_base / DEFAULT_DOCS_SUBDIR_NAME,
        "tests_dir": output_data_base / DEFAULT_TESTS_SUBDIR_NAME,
        "tutorials_dir": output_data_base / DEFAULT_TUTORIALS_SUBDIR_NAME,
        # User-managed/project-level EOS and TOV data
        "user_eos_data_dir": output_data_base / DEFAULT_EOS_DATA_SUBDIR_NAME,
        "test_data_dir": output_data_base / DEFAULT_TEST_DATA_SUBDIR_NAME,
        "train_data_dir": output_data_base / DEFAULT_TRAIN_DATA_SUBDIR_NAME,
        "val_data_dir": output_data_base / DEFAULT_VAL_DATA_SUBDIR_NAME,
        "generated_eos_files_dir": output_data_base / DEFAULT_EOS_FILES_SUBDIR_NAME,
        # Current dynamic EOS/TOV/generated paths
        "current_eos_input_dir": output_data_base
        / DEFAULT_EOS_FILES_SUBDIR_NAME
        / eos_folder_name,
        "current_tov_data_dir": output_data_base
        / DEFAULT_TOV_DATA_SUBDIR_NAME
        / eos_folder_name,
        "current_train_data_dir": output_data_base
        / DEFAULT_TRAIN_DATA_SUBDIR_NAME
        / eos_folder_name,
        "current_val_data_dir": output_data_base
        / DEFAULT_VAL_DATA_SUBDIR_NAME
        / eos_folder_name,
        "current_test_data_dir": output_data_base
        / DEFAULT_TEST_DATA_SUBDIR_NAME
        / eos_folder_name,
        # Kept old names for minimal compatibility changes
        "qeos_path_specific": output_data_base / DEFAULT_EOS_FILES_SUBDIR_NAME / "QEOS",
        "mseos_path_specific": output_data_base
        / DEFAULT_EOS_FILES_SUBDIR_NAME
        / "MSEOS",
        "qeos_tov_path_specific": output_data_base
        / DEFAULT_TOV_DATA_SUBDIR_NAME
        / "QEOS",
        "mseos_tov_path_specific": output_data_base
        / DEFAULT_TOV_DATA_SUBDIR_NAME
        / "MSEOS",
        "user_tov_data_dir": output_data_base / DEFAULT_TOV_DATA_SUBDIR_NAME,
    }

    if include_slm_paths:
        paths["current_slm_results_dir"] = (
            output_data_base / DEFAULT_RESULTS_SUBDIR_NAME / eos_folder_name / slm_subdir
        )
        paths["current_slm_plots_dir"] = (
            output_data_base / DEFAULT_PLOTS_SUBDIR_NAME / eos_folder_name / slm_subdir
        )
        if is_parametric_run:
            # Parametric tests always live under testData/<EOS>/pSLM
            paths["current_slm_tests_dir"] = (
                output_data_base
                / DEFAULT_TEST_DATA_SUBDIR_NAME
                / eos_folder_name
                / DEFAULT_PSLM_SUBDIR_NAME
            )
        else:
            paths["current_slm_tests_dir"] = (
                output_data_base / DEFAULT_TEST_DATA_SUBDIR_NAME / eos_folder_name
            )

    return paths

create_necessary_dirs(paths: dict[str, Path], additional_dirs: list[Path] | None = None) -> None

Creates necessary directories specified in a dictionary and an optional list.

Iterates through the known output-directory keys of paths (plus any additional_dirs) and creates each directory if it does not already exist (mkdir(parents=True, exist_ok=True)).

Parameters:

Name Type Description Default
paths dict[str, Path]

A dictionary as returned by :func:get_paths, mapping string identifiers to directories.

 required
additional_dirs list[Path]

Additional directories to create (e.g., user-managed data or model directories). Defaults to None.

 None

Returns:

Name Type Description
None None

The function modifies the filesystem but does not return a value.
Source code in src/slmemulator/config.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
def create_necessary_dirs(
    paths: dict[str, Path],
    additional_dirs: list[Path] | None = None,
) -> None:
    """
    Creates necessary directories specified in a dictionary and an optional list.

    Iterates through the known output-directory keys of ``paths`` (plus any
    ``additional_dirs``) and creates each directory if it does not already
    exist (``mkdir(parents=True, exist_ok=True)``).

    Parameters:
        paths (dict[str, pathlib.Path]): A dictionary as returned by
            :func:`get_paths`, mapping string identifiers to directories.
        additional_dirs (list[pathlib.Path], optional): Additional directories
            to create (e.g., user-managed data or model directories).
            Defaults to None.

    Returns:
        None: The function modifies the filesystem but does not return a value.
    """
    output_directory_keys = [
        "plots_dir",
        "results_dir",
        "docs_dir",
        "tests_dir",
        "tutorials_dir",
        "user_eos_data_dir",
        "user_tov_data_dir",
        "test_data_dir",
        "train_data_dir",
        "generated_eos_files_dir",
        # The 'current' paths are the effective targets
        "current_eos_input_dir",
        "current_tov_data_dir",
        "current_slm_results_dir",
        "current_slm_plots_dir",
        "current_slm_tests_dir",
        # Kept for backward compatibility
        "qeos_path_specific",
        "mseos_path_specific",
        "qeos_tov_path_specific",
        "mseos_tov_path_specific",
    ]

    dirs_to_create = {
        paths[key]
        for key in output_directory_keys
        if key in paths and isinstance(paths[key], Path)
    }
    if additional_dirs:
        dirs_to_create.update(p for p in additional_dirs if isinstance(p, Path))

    for path in dirs_to_create:
        try:
            path.mkdir(parents=True, exist_ok=True)
        except OSError as e:
            print(f"Error creating directory {path}: {e}")

pSLM

Parametric SLM built on Banach-GRIM kernel interpolation.

Fits directly on TOV data files (columns: radius, central pressure, mass [, k2]) computed at known EOS parameters, and predicts the log-curves at new parameter values by kernel interpolation of the full curve matrices across min-max-normalized parameter space.

This replaces the earlier k-nearest-neighbour version: k-NN averaging of DMD components is piecewise constant in parameter space and ill-posed for the eigenpairs (mode ordering/sign ambiguity), whereas kernel interpolation of the curves is exact at the training points and smooth in between. The DMD components returned by predict() are computed from the predicted curve, so they are always self-consistent.

For an end-to-end emulator that also generates the training data (EOS generation + TOV solves), see slmemulator.TOVEmulator.

ParametricSLM(fileList, filePath=None, tidal=False, params=None, reg=1e-10, length_scale=None, error_threshold=1e-06, max_r=None)

Parametric SLM over a set of TOV data files.

Example

pslm = ParametricSLM(fileList, filePath, tidal=True, params=[[300, 0.1], [300, 0.3], ...]) pslm.fit() Phi, omega, eigs, b, Xdmd, t = pslm.predict([400.0, 0.2])

Parameters

fileList : sequence of str or Path TOV data files; each contains columns radius (km), central pressure (MeV/fm^3), mass (M_sun) and, if tidal, k2. All files must share the same number of rows. filePath : str or Path, optional Directory prepended to relative file names. tidal : bool Whether the files carry a 4th (k2) column to be used. Default False. params : array-like (n_files, n_params), optional EOS parameters of each file, in the same order as fileList. If not given, parameters are parsed from the file names (every underscore- separated token of the stem that parses as a float) — prefer passing them explicitly, since file-name parsing cannot recover values whose decimal points were replaced by underscores. reg : float Tikhonov regularization of the kernel interpolant. Default 1e-10. length_scale : float, optional RBF length scale in normalized parameter units (default: median pairwise distance of the normalized training parameters). error_threshold, max_r : Passed to SLM() when computing DMD components of predicted curves.

Source code in src/slmemulator/pSLM.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def __init__(
    self,
    fileList,
    filePath=None,
    tidal=False,
    params=None,
    reg=1e-10,
    length_scale=None,
    error_threshold=1e-6,
    max_r=None,
):
    base = Path(filePath) if filePath is not None else None
    self.fileList = [
        (base / f if base is not None and not Path(f).is_absolute() else Path(f))
        for f in fileList
    ]
    self.tidal = bool(tidal)
    self.reg = float(reg)
    self.length_scale = length_scale
    self.error_threshold = error_threshold
    self.max_r = max_r

    self._given_params = None if params is None else np.atleast_2d(
        np.asarray(params, dtype=float)
    )

    # Fitted state
    self.params = None            # (n_files, n_params)
    self.curves_log = None        # (n_files, n_quantities, n_points)
    self.n = None                 # number of quantities
    self.mm1 = None               # number of points per curve
    self._grim = None
    self._p_min = None
    self._p_range = None

params_from_filename(file_path) staticmethod

Extract EOS parameters from a file name: every underscore-separated token of the stem that parses as a float (e.g. 'MR_3.00e+02_1.00e-01.txt' -> [300.0, 0.1]).

Source code in src/slmemulator/pSLM.py 
108
109
110
111
112
113
114
115
116
117
118
119
@staticmethod
def params_from_filename(file_path):
    """Extract EOS parameters from a file name: every underscore-separated
    token of the stem that parses as a float (e.g.
    'MR_3.00e+02_1.00e-01.txt' -> [300.0, 0.1])."""
    values = []
    for token in Path(file_path).stem.split("_"):
        try:
            values.append(float(token))
        except ValueError:
            continue
    return values

fit()

Load all training files and fit the kernel interpolant.

Source code in src/slmemulator/pSLM.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
def fit(self):
    """Load all training files and fit the kernel interpolant."""
    curves, params = [], []
    for i, file_path in enumerate(self.fileList):
        if not file_path.exists():
            print(f"Skipping non-existent file: {file_path}")
            continue
        curve = self._load_curve(file_path)
        if not np.all(np.isfinite(curve)):
            print(f"Skipping {file_path.name}: non-finite values.")
            continue
        if curves and curve.shape != curves[0].shape:
            print(
                f"Skipping {file_path.name}: inconsistent shape "
                f"{curve.shape} (expected {curves[0].shape})."
            )
            continue
        curves.append(curve)
        if self._given_params is not None:
            params.append(self._given_params[i])
        else:
            parsed = self.params_from_filename(file_path)
            if not parsed:
                raise ValueError(
                    f"Could not parse parameters from '{file_path.name}'; "
                    "pass params= explicitly."
                )
            params.append(parsed)

    if not curves:
        raise RuntimeError("No valid training files processed.")

    self.curves_log = np.asarray(curves, dtype=float)
    self.params = np.asarray(params, dtype=float)
    self.n, self.mm1 = self.curves_log.shape[1], self.curves_log.shape[2]

    self._p_min = self.params.min(axis=0)
    self._p_range = self.params.max(axis=0) - self._p_min
    self._p_range[self._p_range == 0.0] = 1.0

    self._grim = BanachGRIMInterpolator(
        reg=self.reg, length_scale=self.length_scale, center=True
    )
    self._grim.fit(list(self.curves_log), self._scale(self.params), tol=1e-12)
    return self

predict_log(theta)

Interpolated log-curve matrix (n_quantities, n_points) at theta.

Source code in src/slmemulator/pSLM.py 
182
183
184
185
186
def predict_log(self, theta):
    """Interpolated log-curve matrix (n_quantities, n_points) at theta."""
    if self._grim is None:
        raise RuntimeError("Model not fitted. Call .fit() first.")
    return self._grim.evaluate(self._scale(np.atleast_1d(theta)))

predict_curves(theta)

Interpolated curves in linear units (n_quantities, n_points).

Source code in src/slmemulator/pSLM.py 
188
189
190
def predict_curves(self, theta):
    """Interpolated curves in linear units (n_quantities, n_points)."""
    return np.exp(self.predict_log(theta))

predict(theta, dt=1.0)

Predict at parameter vector theta.

Returns (Phi, omega, eigs, b, Xdmd, t) — the same shape of result as the earlier k-NN implementation, but the DMD components are computed from the kernel-interpolated curve, so they are self-consistent. Xdmd is the log-space SLM reconstruction; np.exp(Xdmd.real) gives the physical curves.

Source code in src/slmemulator/pSLM.py 
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
def predict(self, theta, dt=1.0):
    """
    Predict at parameter vector theta.

    Returns (Phi, omega, eigs, b, Xdmd, t) — the same shape of result as
    the earlier k-NN implementation, but the DMD components are computed
    from the kernel-interpolated curve, so they are self-consistent.
    Xdmd is the log-space SLM reconstruction; np.exp(Xdmd.real) gives the
    physical curves.
    """
    X_log = self.predict_log(theta)
    Phi, omega, eigs, b, Xdmd, _S, _r = SLM(
        X_log, dt=dt, error_threshold=self.error_threshold, max_r=self.max_r
    )
    t = np.arange(X_log.shape[1]) * dt
    return Phi, omega, eigs, b, Xdmd, t

gaussian_kernel(x1, x2, sigma=1.0)

Gaussian kernel between two parameter vectors (kept for backward compatibility; the interpolation itself uses banach_grim.rbf_kernel).

Source code in src/slmemulator/pSLM.py 
29
30
31
32
33
34
35
def gaussian_kernel(x1, x2, sigma=1.0):
    """Gaussian kernel between two parameter vectors (kept for backward
    compatibility; the interpolation itself uses banach_grim.rbf_kernel)."""
    x1_arr = np.asarray(x1)
    x2_arr = np.asarray(x2)
    dist_sq = np.sum((x1_arr - x2_arr) ** 2)
    return np.exp(-dist_sq / (2 * sigma**2 + 1e-9))

is_on_boundary(param, param_min, param_max, tolerance=1e-05)

Checks if a parameter set is on the boundary of the parameter space.

Source code in src/slmemulator/pSLM.py 
210
211
212
213
214
215
216
217
218
def is_on_boundary(param, param_min, param_max, tolerance=1e-5):
    """Checks if a parameter set is on the boundary of the parameter space."""
    for i in range(len(param)):
        if (
            abs(param[i] - param_min[i]) < tolerance
            or abs(param[i] - param_max[i]) < tolerance
        ):
            return True
    return False

recombination

Recombination thinning for the Banach GRIM algorithm (arXiv:2205.07495).

The implementation lives in banach_grim.py (so that module stays self-contained); this module re-exports it alongside the Gaussian RBF helper.

recombination_thinning(M, weights, tol=1e-12, max_iter=None)

Reduce the support of nonnegative weights while preserving linear moments.

Given M of shape (r, N) and weights w >= 0 of shape (N,), returns new weights b >= 0 with M @ b == M @ w (to numerical precision) and support of at most r atoms (fewer if M's active columns are rank deficient). This is the recombination step of the Banach GRIM algorithm (Lemma 3.1 of arXiv:2205.07495): with an all-ones row included in M, the total mass sum(b) == sum(w) is preserved as well.

Parameters:

Name Type Description Default
M ndarray

Constraint matrix, shape (r, N). Row k holds the values of the k-th linear functional on the N atoms.

 required
weights ndarray

Nonnegative weights, shape (N,).

 required
tol float

Singular values below tol * largest count as null directions.

 1e-12
max_iter int

Safety cap on elimination steps (default N).

 None

Returns:

Name Type Description
b ndarray

Nonnegative weights, shape (N,), with np.count_nonzero(b) <= r and M @ b == M @ weights.
Source code in src/slmemulator/banach_grim.py 
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
def recombination_thinning(M, weights, tol=1e-12, max_iter=None):
    """
    Reduce the support of nonnegative weights while preserving linear moments.

    Given ``M`` of shape (r, N) and ``weights`` w >= 0 of shape (N,), returns
    new weights b >= 0 with ``M @ b == M @ w`` (to numerical precision) and
    support of at most r atoms (fewer if M's active columns are rank
    deficient). This is the recombination step of the Banach GRIM algorithm
    (Lemma 3.1 of arXiv:2205.07495): with an all-ones row included in M, the
    total mass ``sum(b) == sum(w)`` is preserved as well.

    Parameters:
        M (np.ndarray): Constraint matrix, shape (r, N). Row k holds the
            values of the k-th linear functional on the N atoms.
        weights (np.ndarray): Nonnegative weights, shape (N,).
        tol (float): Singular values below ``tol * largest`` count as null
            directions.
        max_iter (int, optional): Safety cap on elimination steps
            (default N).

    Returns:
        b (np.ndarray): Nonnegative weights, shape (N,), with
            ``np.count_nonzero(b) <= r`` and ``M @ b == M @ weights``.
    """
    M = np.atleast_2d(np.asarray(M, dtype=float))
    b = np.asarray(weights, dtype=float).copy()
    if b.ndim != 1 or M.shape[1] != b.shape[0]:
        raise ValueError(f"Shape mismatch: M {M.shape}, weights {b.shape}.")
    if np.any(b < 0):
        raise ValueError("recombination requires nonnegative weights.")

    support = np.flatnonzero(b > 0.0)
    if support.size <= 1:
        return b

    # Null basis of the active columns, computed once: rows of Vt whose
    # singular value is (numerically) zero span the right null space.
    M_s = M[:, support]
    _, s, Vt = np.linalg.svd(M_s, full_matrices=True)
    rank = int(np.sum(s > tol * max(s[0] if s.size else 1.0, 1.0)))
    V = Vt[rank:].T.copy()  # (s_count, q): columns are null vectors

    if max_iter is None:
        max_iter = b.shape[0]

    for _ in range(min(max_iter, V.shape[1] if V.size else 0)):
        if V.size == 0 or support.size <= 1:
            break

        # Step along the first null direction until a weight hits zero,
        # keeping every weight nonnegative. Positive entries of v limit the
        # step; flip v if it has none.
        v = V[:, 0]
        if not np.any(v > tol):
            v = -v
            V[:, 0] = v
        pos = v > tol
        if not np.any(pos):
            # direction is (numerically) zero: discard it
            V = V[:, 1:]
            continue
        ratios = b[support[pos]] / v[pos]
        j_local_pos = int(np.argmin(ratios))
        gamma = ratios[j_local_pos]
        j_local = int(np.flatnonzero(pos)[j_local_pos])

        b[support] = b[support] - gamma * v
        b[support[j_local]] = 0.0  # exactly zero the eliminated atom
        np.clip(b, 0.0, None, out=b)

        # Update the null basis for the reduced support: use v (which has
        # v[j_local] > 0) to cancel the j_local-component of the remaining
        # directions, then drop v and the eliminated row.
        if V.shape[1] > 1:
            V[:, 1:] -= np.outer(v, V[j_local, 1:] / v[j_local])
        V = np.delete(V[:, 1:], j_local, axis=0)
        keep = np.ones(support.size, dtype=bool)
        keep[j_local] = False
        support = support[keep]
        # guard against numerical blow-up in the eliminated basis
        if V.size:
            norms = np.linalg.norm(V, axis=0)
            good = norms > tol
            V = V[:, good] / np.where(norms[good] > 0, norms[good], 1.0)

    return b

gaussian_rbf(r, epsilon)

Gaussian Radial Basis Function.

Source code in src/slmemulator/recombination.py 
14
15
16
def gaussian_rbf(r, epsilon):
    """Gaussian Radial Basis Function."""
    return np.exp(-((epsilon * r) ** 2))

scaledTOV

This code solves TOV equations for mass radius relations. This can also plot the mass-radius curve.

USE: To use the code, here are the steps: 1) Include the file in your main code e.g. import tov_class as tc 2) Load the EoS using the ToV loader, tc.ToV(filename, arraysize) 3) call the solver as tc.ToV.mass_radius(min_pressure, max_pressure) 4) To plot, follow the code in main() on creating the dictionary of inputs

Updates: Solves ToV, can only take inputs of pressure (MeV/fm^3), energy density in MeV, baryon density in fm^-3 in ascending order.

TOV(filename, imax)

Solves TOV equations and gives data-table, mass-radius plot and max. mass, central pressure and central density by loading an EoS datafile.

Source code in src/slmemulator/scaledTOV.py 
47
48
49
50
51
52
53
54
def __init__(self, filename, imax):
    self.file = np.loadtxt(filename, skiprows=3)
    self.e_in = self.file[:, 1] / eps0  # Scaled Energy density
    self.p_in = self.file[:, 2] / pres0  # Scaled pressure
    self.nb_in = self.file[:, 0]
    self.imax = imax
    self.radius = np.empty(self.imax)
    self.mass = np.empty(self.imax)

pressure_from_nb(nb: Union[float, np.ndarray]) -> Union[float, np.ndarray]

Evaluates scaled pressure (\(P/P_0\)) given the baryon number density (\(n_B\)) using linear interpolation of the loaded Equation of State (EOS) data.

This function uses :func:scipy.interpolate.interp1d to create an interpolating function based on the input baryon number density (:attr:self.nb_in) and scaled pressure (:attr:self.p_in) from the loaded EOS table.

Parameters:

Name Type Description Default
nb float or ndarray

The baryon number density (or an array of densities) at which to

required

Returns:

Type Description
Union[float, ndarray]

float or numpy.ndarray:
Union[float, ndarray]

The interpolated scaled pressure (\(P/P_0\)) value(s) corresponding to

Union[float, ndarray]

the input number density nb.
Source code in src/slmemulator/scaledTOV.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def pressure_from_nb(self, nb: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
    """
    Evaluates scaled pressure ($P/P_0$) given the baryon number density ($n_B$) 
    using linear interpolation of the loaded Equation of State (EOS) data.

    This function uses :func:`scipy.interpolate.interp1d` to create an 
    interpolating function based on the input baryon number density 
    (:attr:`self.nb_in`) and scaled pressure (:attr:`self.p_in`) from the 
    loaded EOS table.

    Parameters:
        nb (float or numpy.ndarray): The baryon number density (or an array of densities) at which to 
        evaluate the corresponding scaled pressure.

    Returns:
        float or numpy.ndarray:
        The interpolated scaled pressure ($P/P_0$) value(s) corresponding to 
        the input number density ``nb``.
    """
    p1 = interp1d(
        self.nb_in, self.p_in, axis=0, kind="linear", fill_value="extrapolate"
    )
    return p1(nb)

energy_from_pressure(pressure: Union[float, np.ndarray]) -> Union[float, np.ndarray]

Evaluates scaled energy density (\(\epsilon/\epsilon_0\)) given the scaled pressure (\(P/P_0\)) using linear interpolation of the loaded Equation of State (EOS) data.

This method handles pressures near zero with a special case for numerical stability.

pressure (float or numpy.ndarray): The scaled pressure (\(P/P_0\)) value(s) at which to evaluate the corresponding scaled energy density.

Returns:

Type Description
Union[float, ndarray]

float or numpy.ndarray: The interpolated scaled energy density (\(\epsilon/\epsilon_0\)) value(s).
Source code in src/slmemulator/scaledTOV.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def energy_from_pressure(self, pressure: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
    r"""
    Evaluates scaled energy density ($\epsilon/\epsilon_0$) given the scaled 
    pressure ($P/P_0$) using linear interpolation of the loaded Equation of 
    State (EOS) data.

    This method handles pressures near zero with a special case for numerical stability.

    Parameters:
    pressure (float or numpy.ndarray): The scaled pressure ($P/P_0$) value(s) at which to evaluate the 
        corresponding scaled energy density.

    Returns:
        float or numpy.ndarray: The interpolated scaled energy density ($\epsilon/\epsilon_0$) value(s).

    """
    plow = 1e-10 / pres0
    if pressure < plow:
        return 2.6e-310
    else:
        e1 = interp1d(
            self.p_in, self.e_in, axis=0, kind="linear", fill_value="extrapolate"
        )
        return e1(pressure)

pressure_from_energy(energy: Union[float, np.ndarray]) -> Union[float, np.ndarray]

Evaluates scaled pressure (\(P/P_0\)) given the scaled energy density (\(\epsilon/\epsilon_0\)) using linear interpolation of the loaded Equation of State (EOS) data.

This function defines the inverse of the \(\epsilon(P)\) relation.

energy (float or numpy.ndarray): The scaled energy density (\(\epsilon/\epsilon_0\)) value(s) at which to evaluate the corresponding scaled pressure.

Returns:

Type Description
Union[float, ndarray]

float or numpy.ndarray: The interpolated scaled pressure (\(P/P_0\)) value(s) corresponding to

Union[float, ndarray]

the input scaled energy density.
Source code in src/slmemulator/scaledTOV.py 
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def pressure_from_energy(self, energy: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
    r"""
    Evaluates scaled pressure ($P/P_0$) given the scaled energy density 
    ($\epsilon/\epsilon_0$) using linear interpolation of the loaded 
    Equation of State (EOS) data.

    This function defines the inverse of the $\epsilon(P)$ relation.

    Parameters:
    energy (float or numpy.ndarray): The scaled energy density ($\epsilon/\epsilon_0$) value(s) at which 
        to evaluate the corresponding scaled pressure.

    Returns:
        float or numpy.ndarray: The interpolated scaled pressure ($P/P_0$) value(s) corresponding to 
        the input scaled energy density.
    """
    p1 = interp1d(
        self.e_in, self.p_in, axis=0, kind="linear", fill_value="extrapolate"
    )
    return p1(energy)

baryon_from_energy(energy)

Evaluate number density from energy using interpolation

Source code in src/slmemulator/scaledTOV.py 
126
127
128
129
130
131
def baryon_from_energy(self, energy):
    """Evaluate number density from energy using interpolation"""
    n1 = interp1d(
        self.e_in, self.nb_in, axis=0, kind="linear", fill_value="extrapolate"
    )
    return n1(energy)

RK4(f, x0, t0, te, N)

A simple RK4 solver to avoid overhead of calculating with solve_ivp or any other adaptive step-size function.

Example

tov.RK4(f=func, x0=1., t0=1., te=10., N=100)

Parameters:

Name Type Description Default
f func

A Python function for the ODE(s) to be solved. Able to solve N coupled ODEs.

 required
x0 float

Guess for the function(s) to be solved.

 required
t0 float

Initial point of the grid.

 required
te float

End point of the grid.

 required
N int

The number of steps to take in the range (te-t0).

 required

Returns:

Name Type Description
times array

The grid of solution steps.
solution array

The solutions of each function at each point in the grid.
Source code in src/slmemulator/scaledTOV.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def RK4(self, f, x0, t0, te, N):
    r"""
    A simple RK4 solver to avoid overhead of
    calculating with solve_ivp or any other
    adaptive step-size function.

    Example:
        tov.RK4(f=func, x0=1., t0=1., te=10., N=100)

    Parameters:
        f (func): A Python function for the ODE(s) to be solved.
            Able to solve N coupled ODEs.

        x0 (float): Guess for the function(s) to be solved.

        t0 (float): Initial point of the grid.

        te (float): End point of the grid.

        N (int): The number of steps to take in the range (te-t0).

    Returns:
        times (array): The grid of solution steps.

        solution (array): The solutions of each function
            at each point in the grid.
    """

    h = (te - t0) / N
    times = np.arange(t0, te + h, h)
    solution = []
    x = x0

    for t in times:
        solution.append(np.array(x).T)
        k1 = h * f(t, x)
        k2 = h * f(t + 0.5 * h, x + 0.5 * k1)
        k3 = h * f(t + 0.5 * h, x + 0.5 * k2)
        k4 = h * f(t + h, x + k3)
        x += (k1 + 2 * (k2 + k3) + k4) / 6

    solution = np.asarray(solution, dtype=np.float64).T

    return times, solution

TOV_class

TOV(eos_filepath=None, tidal=False, solver='RK4', solve_ivp_kwargs=None, sol_pts=4000)

inertia using RK4. Also includes uncertainty quantification techniques through the highest posterior density interval (HPD or HDI) calculation. Able to accept one EOS from a single curve or draws from an EOS, such as from a Gaussian Process.

Example

tov = TOV(eos_filepath='path/to/eos', tidal=True)

Parameters:

Name Type Description Default
eos_filepath str or Path

The path to the EOS data table to be used. Supported formats: .table (compOSE), .dat/.txt (columns nb, E, P[, cs2]), and .npz (keys density, edens, pres, cs2).

 None
tidal bool

Whether to calculate tidal deformability or not. Default is False.

 False
solver str

One of "RK4", "RK2", "euler", or "solve_ivp". Default is "RK4".

 'RK4'
solve_ivp_kwargs dict

Keyword arguments passed to scipy.integrate.solve_ivp (only used when solver="solve_ivp").

 None
sol_pts int

Number of integration steps for the fixed-step solvers. Default is 4000.

 4000
Source code in src/slmemulator/TOV_class.py 
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def __init__(
    self,
    eos_filepath=None,
    tidal=False,
    solver="RK4",
    solve_ivp_kwargs=None,
    sol_pts=4000,  # this is a lot
):
    r"""
    Class to calculate the Tolman-Oppenheimer-Volkoff equations,
    including options for the tidal deformability and moment of
    inertia using RK4. Also includes uncertainty quantification
    techniques through the highest posterior density interval (HPD
    or HDI) calculation. Able to accept one EOS from a single curve
    or draws from an EOS, such as from a Gaussian Process.

    Example:
        tov = TOV(eos_filepath='path/to/eos', tidal=True)

    Parameters:
        eos_filepath (str or Path): The path to the EOS data table to be
            used. Supported formats: .table (compOSE), .dat/.txt
            (columns nb, E, P[, cs2]), and .npz (keys density, edens,
            pres, cs2).

        tidal (bool): Whether to calculate tidal deformability or not.
            Default is False.

        solver (str): One of "RK4", "RK2", "euler", or "solve_ivp".
            Default is "RK4".

        solve_ivp_kwargs (dict): Keyword arguments passed to
            scipy.integrate.solve_ivp (only used when solver="solve_ivp").

        sol_pts (int): Number of integration steps for the fixed-step
            solvers. Default is 4000.
    """

    # assign class variables
    self.tidal = tidal
    self.solver = solver
    self.sol_pts = sol_pts
    self.solve_ivp_kwargs = solve_ivp_kwargs  # only used in solve_ivp
    self.tol = 1e-9  # only used in solve_ivp

    # assign scaled variables
    self.eps0 = 1.285e3  # MeV fm-3
    self.pres0 = self.eps0
    self.mass0 = 2.837  # solar masses
    self.rad0 = 8.378  # km

    if eos_filepath is None:
        raise ValueError("No file specified.")

    eos_path = Path(eos_filepath)
    self.eos_file = str(eos_filepath)
    self.eos_name = eos_path.stem
    self.eos_file_extension = eos_path.suffix

    # for data from compOSE
    if self.eos_file_extension == ".table":
        eos_data = np.loadtxt(eos_filepath)

        # check which column is which
        if eos_data.T[4][-1] < eos_data.T[3][-1]:
            self.eps_array = eos_data.T[3] / self.eps0
            self.pres_array = eos_data.T[4] / self.pres0
        else:
            self.eps_array = eos_data.T[4] / self.eps0
            self.pres_array = eos_data.T[3] / self.pres0

        self.nB_array = eos_data.T[1]

        if tidal:
            self.cs2_array = np.gradient(
                self.pres_array, self.eps_array, edge_order=2
            )
        else:
            self.cs2_array = np.zeros(len(self.pres_array))

        # keep unscaled for use in density calculation
        self.pres_array_unscaled = eos_data.T[4]

    # for data generated by codes
    elif self.eos_file_extension in (".dat", ".txt"):
        print("Loading EOS data from file: ", eos_filepath)
        eos_data = np.loadtxt(eos_filepath)
        self.eps_array = eos_data.T[1] / self.eps0
        self.pres_array = eos_data.T[2] / self.pres0
        self.nB_array = eos_data.T[0]
        self.cs2_array = eos_data.T[3] if eos_data.shape[1] > 3 else None

        # keep unscaled for use in density calculation
        self.pres_array_unscaled = eos_data.T[2]

    elif self.eos_file_extension == ".npz":
        eos_data = np.load(eos_filepath)

        # assign to arrays based on header names
        self.eps_array = eos_data["edens"] / self.eps0
        self.pres_array = eos_data["pres"] / self.pres0
        self.nB_array = eos_data["density"]

        if "cs2" in eos_data:
            self.cs2_array = eos_data["cs2"]
        elif tidal:
            # compute dP/deps per draw (columns) or for a single curve
            self.cs2_array = self._dpdeps_gradient(
                self.pres_array, self.eps_array
            )
        else:
            self.cs2_array = None

        # keep unscaled for use in density calculation
        self.pres_array_unscaled = eos_data["pres"]

    # HDF5 files with (possibly multiple) EOS draws in columns
    elif self.eos_file_extension == ".h5":
        print("Loading EOS data from file: ", eos_filepath)
        with h5py.File(eos_filepath, "r") as f:
            self.nB_array = f["dens"][:]
            self.eps_array = f["edens"][:] / self.eps0
            self.pres_array = f["pressure"][:] / self.pres0

        # Sound speed from spline derivatives, cs2 = (dP/dnB)/(deps/dnB).
        # A stored 'soundspeed' dataset is deliberately NOT used: for GP
        # ensembles it is sampled independently of the P and eps draws,
        # and that thermodynamic inconsistency destabilizes the tidal
        # (y) equation. The derivative of the actual P, eps draws is the
        # consistent choice (and matches the production solver).
        nB = (
            self.nB_array
            if self.nB_array.ndim == 1
            else self.nB_array[:, 0]
        )
        pres = np.atleast_2d(self.pres_array.T).T  # (n_points, n_draws)
        eps = np.atleast_2d(self.eps_array.T).T
        dpdeps = np.zeros_like(pres)
        for j in range(pres.shape[1]):
            dP_dnB = CubicSpline(nB, pres[:, j]).derivative()(nB)
            dEps_dnB = CubicSpline(nB, eps[:, j]).derivative()(nB)
            dpdeps[:, j] = dP_dnB / dEps_dnB
        self.cs2_array = dpdeps.reshape(self.pres_array.shape)

        # keep unscaled for use in density calculation
        self.pres_array_unscaled = self.pres_array * self.pres0

    else:
        raise ValueError(
            f"Unsupported EOS file extension: '{self.eos_file_extension}'."
        )

RK4(f, x0, t0, te, N)

A simple RK4 solver to avoid overhead of calculating with solve_ivp or any other adaptive step-size function.

Example

tov.RK4(f=func, x0=1., t0=1., te=10., N=100)

Parameters:

Name Type Description Default
f func

A Python function for the ODE(s) to be solved. Able to solve N coupled ODEs.

 required
x0 float or array - like

Initial value(s) of the function(s) to be solved.

 required
t0 float

Initial point of the grid.

 required
te float

End point of the grid.

 required
N int

The number of steps to take in the range (te-t0).

 required

Returns:

Name Type Description
times array

The grid of solution steps.
solution array

The solutions of each function at each point in the grid.
Source code in src/slmemulator/TOV_class.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def RK4(self, f, x0, t0, te, N):
    r"""
    A simple RK4 solver to avoid overhead of
    calculating with solve_ivp or any other
    adaptive step-size function.

    Example:
        tov.RK4(f=func, x0=1., t0=1., te=10., N=100)

    Parameters:
        f (func): A Python function for the ODE(s) to be solved.
            Able to solve N coupled ODEs.

        x0 (float or array-like): Initial value(s) of the function(s)
            to be solved.

        t0 (float): Initial point of the grid.

        te (float): End point of the grid.

        N (int): The number of steps to take in the range (te-t0).

    Returns:
        times (array): The grid of solution steps.

        solution (array): The solutions of each function
            at each point in the grid.
    """

    h = (te - t0) / N
    times = np.arange(t0, te + h, h)
    solution = []
    x = np.asarray(x0, dtype=np.float64).copy()

    for t in times:
        solution.append(x.copy())
        k1 = h * f(t, x)
        k2 = h * f(t + 0.5 * h, x + 0.5 * k1)
        k3 = h * f(t + 0.5 * h, x + 0.5 * k2)
        k4 = h * f(t + h, x + k3)
        x += (k1 + 2 * (k2 + k3) + k4) / 6

    solution = np.asarray(solution, dtype=np.float64).T

    return times, solution

RK2(f, x0, t0, te, N)

A simple RK2 solver using Heun's method. This is a low-fidelity solver.

Example

tov.RK2(f=func, x0=1., t0=1., te=10., N=100)

Parameters:

Name Type Description Default
f func

A Python function for the ODE(s) to be solved. Able to solve N coupled ODEs.

 required
x0 float or array - like

Initial value(s) of the function(s) to be solved.

 required
t0 float

Initial point of the grid.

 required
te float

End point of the grid.

 required
N int

The number of steps to take in the range (te-t0).

 required

Returns:

Name Type Description
times array

The grid of solution steps.
solution array

The solutions of each function at each point in the grid.
Source code in src/slmemulator/TOV_class.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def RK2(self, f, x0, t0, te, N):
    r"""
    A simple RK2 solver using Heun's method.
    This is a low-fidelity solver.

    Example:
        tov.RK2(f=func, x0=1., t0=1., te=10., N=100)

    Parameters:
        f (func): A Python function for the ODE(s) to be solved.
            Able to solve N coupled ODEs.

        x0 (float or array-like): Initial value(s) of the function(s)
            to be solved.

        t0 (float): Initial point of the grid.

        te (float): End point of the grid.

        N (int): The number of steps to take in the range (te-t0).

    Returns:
        times (array): The grid of solution steps.

        solution (array): The solutions of each function
            at each point in the grid.
    """

    h = (te - t0) / N
    times = np.arange(t0, te + h, h)
    solution = []
    x = np.asarray(x0, dtype=np.float64).copy()

    for t in times:
        solution.append(x.copy())
        k1 = f(t, x)
        k2 = f(t + h, x + k1 * h)
        x += h * (k1 + k2) * 0.5

    solution = np.asarray(solution, dtype=np.float64).T

    return times, solution

euler(f, x0, t0, te, N)

A simple forward Euler solver to avoid overhead of calculating with solve_ivp or any other adaptive step-size function. This is a low fidelity solver!

Example

tov.euler(f=func, x0=1., t0=1., te=10., N=100)

Parameters:

Name Type Description Default
f func

A Python function for the ODE(s) to be solved. Able to solve N coupled ODEs.

 required
x0 float or array - like

Initial value(s) of the function(s) to be solved.

 required
t0 float

Initial point of the grid.

 required
te float

End point of the grid.

 required
N int

The number of steps to take in the range (te-t0).

 required

Returns:

Name Type Description
times array

The grid of solution steps.
solution array

The solutions of each function at each point in the grid.
Source code in src/slmemulator/TOV_class.py 
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
def euler(self, f, x0, t0, te, N):
    r"""
    A simple forward Euler solver to avoid overhead of
    calculating with solve_ivp or any other
    adaptive step-size function.
    This is a low fidelity solver!

    Example:
        tov.euler(f=func, x0=1., t0=1., te=10., N=100)

    Parameters:
        f (func): A Python function for the ODE(s) to be solved.
            Able to solve N coupled ODEs.

        x0 (float or array-like): Initial value(s) of the function(s)
            to be solved.

        t0 (float): Initial point of the grid.

        te (float): End point of the grid.

        N (int): The number of steps to take in the range (te-t0).

    Returns:
        times (array): The grid of solution steps.

        solution (array): The solutions of each function
            at each point in the grid.
    """

    h = (te - t0) / N
    times = np.arange(t0, te + h, h)
    solution = []
    x = np.asarray(x0, dtype=np.float64).copy()

    for t in times:
        solution.append(x.copy())
        x += h * f(t, x)

    solution = np.asarray(solution, dtype=np.float64).T

    return times, solution

tov_equations_scaled(x, y0)

The Tolman-Oppenheimer-Volkoff equations in scaled format, to be solved with the RK4 routine. If selected, the tidal deformability and moment of inertia will be included and solved simultaneously.

Example

tov.tov_equations_scaled(x=0.2, y0=[m_init, p_init])

Parameters:

Name Type Description Default
x float

A point in the scaled radius grid.

 required
y0 list

The list of initial guesses for each function solved.

 required

Returns:

Type Description

The solutions, in array format, of each function to be solved.
Source code in src/slmemulator/TOV_class.py 
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
def tov_equations_scaled(self, x, y0):
    r"""
    The Tolman-Oppenheimer-Volkoff equations in scaled format, to be
    solved with the RK4 routine. If selected, the tidal deformability
    and moment of inertia will be included and solved
    simultaneously.

    Example:
        tov.tov_equations_scaled(x=0.2,
            y0=[m_init, p_init])

    Parameters:
        x (float): A point in the scaled radius grid.

        y0 (list): The list of initial guesses for each function
            solved.

    Returns:
        The solutions, in array format, of each function to be
            solved.
    """

    # unpack the initial conditions
    if self.tidal:
        pres, mass, y = y0
        cs2 = self.cs2_interp(pres)
    else:
        pres, mass = y0

    eps = self.eps_interp(pres)

    # must also receive monotonically increasing P(n) results
    if pres > 0.0:

        # pressure equation
        dpdx = (
            -0.5 * (pres + eps) * (mass + 3 * x**3.0 * pres) / (x**2.0 - x * mass)
        )

        # mass equation
        dmdx = 3.0 * x**2.0 * eps

        # tidal deformability equation
        if self.tidal:
            f = self.f_x(x, mass, pres, eps)
            q = self.q_x(x, mass, pres, eps, cs2)

            dydx = -(1.0 / x) * (y*y + f * y + q)

            # check for invalid tidal run
            # if not np.isfinite(dydx):
            #     raise RuntimeError("Invalid tidal run.")

    else:
        dpdx = 0.0
        dmdx = 0.0

        if self.tidal:
            dydx = 0.0

    if self.tidal:
        return np.array([dpdx, dmdx, dydx], dtype=np.float64)

    return np.array([dpdx, dmdx], dtype=np.float64)

f_x(x, mass, pres, eps)

A function in the tidal deformability calculation.

Example

tov.f_x(x=0.2, mass=1.06, pres=2.34, eps=6.0)

Parameters:

Name Type Description Default
x float

The current gridpoint in scaled radius.

 required
mass float

The current mass.

 required
pres float

The current pressure from the EOS.

 required
eps float

The current energy density from the EOS.

 required

Returns:

Type Description

The value of F(x) at the current radius.
Source code in src/slmemulator/TOV_class.py 
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
def f_x(self, x, mass, pres, eps):
    r"""
    A function in the tidal deformability calculation.

    Example:
        tov.f_x(x=0.2, mass=1.06, pres=2.34, eps=6.0)

    Parameters:
        x (float): The current gridpoint in scaled radius.

        mass (float): The current mass.

        pres (float): The current pressure from the EOS.

        eps (float): The current energy density from the EOS.

    Returns:
        The value of F(x) at the current radius.
    """

    one = 1.0 - (3.0 / 2.0) * (eps - pres) * x**2.0
    two = 1.0 - mass / x
    return one / two

q_x(x, mass, pres, eps, cs2)

A function in the calculation of the tidal deformability.

Example

tov.q_x(x=0.1, mass=2.0, pres=1.0, eps=3.0, cs2=0.33)

Parameters:

Name Type Description Default
x float

The current gridpoint in scaled radius.

 required
mass float

The current mass.

 required
pres float

The current pressure from the EOS.

 required
eps float

The current energy density from the EOS.

 required
cs2 float

The current speed of sound from the EOS.

 required

Returns:

Type Description

The value of Q(x) at the current radius.
Source code in src/slmemulator/TOV_class.py 
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
def q_x(self, x, mass, pres, eps, cs2):
    r"""
    A function in the calculation of the tidal deformability.

    Example:
        tov.q_x(x=0.1, mass=2.0, pres=1.0, eps=3.0, cs2=0.33)

    Parameters:
        x (float): The current gridpoint in scaled radius.

        mass (float): The current mass.

        pres (float): The current pressure from the EOS.

        eps (float): The current energy density from the EOS.

        cs2 (float): The current speed of sound from the EOS.

    Returns:
        The value of Q(x) at the current radius.
    """
    pre = (3.0 / 2.0) * x**2.0 / (1.0 - mass / x)
    one = 5.0 * eps + 9.0 * pres + ((eps + pres) / cs2) - (4.0 / x**2.0)
    two = mass + 3.0 * x**3.0 * pres
    three = x - mass
    return pre * one - (two / three) ** 2.0

tidal_def(yR, mass, radius)

The calculation of the tidal deformability, Lambda, and the tidal Love number, k2. This function is calculated after the RK4 routine has been completed.

Example

tov.tidal_def(yR=np.array, mass=np.array, radius=np.array)

Parameters:

Name Type Description Default
yR float

The array of y at the maximum radii points.

 required
mass float

The array of mass at the maximum radii.

 required
radius float

The maximum radii array.

 required

Returns:

Name Type Description
tidal_deform array

The tidal deformability solved at each point in the maximum radius.
k2 array

The value of the Love number calculated at the compactness M/R and the value of y at maximum radius.
Source code in src/slmemulator/TOV_class.py 
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
def tidal_def(self, yR, mass, radius):
    r"""
    The calculation of the tidal deformability, Lambda, and
    the tidal Love number, k2. This function is calculated after
    the RK4 routine has been completed.

    Example:
        tov.tidal_def(yR=np.array, mass=np.array, radius=np.array)

    Parameters:
        yR (float): The array of y at the maximum radii points.

        mass (float): The array of mass at the maximum radii.

        radius (float): The maximum radii array.

    Returns:
        tidal_deform (array): The tidal deformability solved at
            each point in the maximum radius.

        k2 (array): The value of the Love number calculated at the
            compactness M/R and the value of y at maximum radius.
    """

    # love number calculation
    beta = (mass / radius) * (self.rad0 / (2.0 * self.mass0))
    k2 = (
        (8.0 / 5.0)
        * beta**5.0
        * (1.0 - 2.0 * beta) ** 2.0
        * (2.0 - yR + 2.0 * beta * (yR - 1.0))
        * (
            2.0 * beta * (6.0 - 3.0 * yR + 3.0 * beta * (5.0 * yR - 8.0))
            + 4.0
            * beta**3.0
            * (
                13.0
                - 11.0 * yR
                + beta * (3.0 * yR - 2.0)
                + 2.0 * beta**2.0 * (1.0 + yR)
            )
            + 3.0
            * (1.0 - 2.0 * beta) ** 2.0
            * (2.0 - yR + 2.0 * beta * (yR - 1.0))
            * np.log(1.0 - 2.0 * beta)
        )
        ** (-1.0)
    )

    # tidal deformability calculation
    tidal_deform = (
        (2.0 / 3.0) * k2 * (2.0 * self.mass0 * radius / (self.rad0 * mass)) ** 5.0
    )

    return tidal_deform, k2

tovsolve(pcent)

Solves the TOV equations for a single central pressure and returns the stellar profile up to the surface.

Parameters:

Name Type Description Default
pcent float

Central pressure in scaled units (P/pres0).

 required

Returns:

Name Type Description
solns ndarray

Columns of radius (km), pressure (MeV/fm^3) and mass (M_sun); the metric function y is appended as a fourth column when tidal is True.
Source code in src/slmemulator/TOV_class.py 
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
def tovsolve(self, pcent):
    r"""
    Solves the TOV equations for a single central pressure and returns
    the stellar profile up to the surface.

    Parameters:
        pcent (float): Central pressure in scaled units (P/pres0).

    Returns:
        solns (np.ndarray): Columns of radius (km), pressure (MeV/fm^3)
            and mass (M_sun); the metric function y is appended as a
            fourth column when tidal is True.
    """
    self._set_eos_interpolants()

    initial = [pcent, 1e-12]
    if self.tidal:
        initial.append(2.0)

    xval, sol = self.RK4(
        self.tov_equations_scaled, initial, 1e-3, 4.0, self.sol_pts
    )

    # truncate at the stellar surface (pressure drops to ~zero)
    positive = np.flatnonzero(sol[0] > 1e-10)
    surface = positive[-1] if positive.size > 0 else 0

    columns = [
        xval[: surface + 1] * self.rad0,  # radius
        sol[0][: surface + 1] * self.pres0,  # pressure
        sol[1][: surface + 1] * self.mass0,  # mass
    ]
    if self.tidal:
        columns.append(sol[2][: surface + 1])

    return np.column_stack(columns)

tov_routine(verbose=False, write_to_file=False)

The TOV routine to solve each set of coupled ODEs and to output the quantities needed to display the M-R curve, as well as the tidal deformability and moment of inertia if desired.

Example

tov.tov_routine(verbose=True, write_to_file=True)

Parameters:

Name Type Description Default
verbose bool

Whether to plot quantities and display the full maximum mass array. Default is False.

 False
write_to_file bool

Choice to write the TOV results to a file located in a folder of the user's choice. Default is False.

 False

Returns:

Type Description

self.total_radius (array): The array of total maximum radius values.

self.total_pres_central (array): The array of total central pressure values.

self.total_mass (array): The array of total maximum mass values.

When tidal is True, self.k2 and self.tidal_deformability are

returned as well.
Source code in src/slmemulator/TOV_class.py 
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
def tov_routine(self, verbose=False, write_to_file=False):
    r"""
    The TOV routine to solve each set of coupled ODEs and to output
    the quantities needed to display the M-R curve, as well as the
    tidal deformability and moment of inertia if desired.

    Example:
        tov.tov_routine(verbose=True, write_to_file=True)

    Parameters:
        verbose (bool): Whether to plot quantities and display
            the full maximum mass array. Default is False.

        write_to_file (bool): Choice to write the TOV results to
            a file located in a folder of the user's choice.
            Default is False.

    Returns:
        self.total_radius (array): The array of total maximum
            radius values.

        self.total_pres_central (array): The array of total
            central pressure values.

        self.total_mass (array): The array of total
            maximum mass values.

        When tidal is True, self.k2 and self.tidal_deformability are
        returned as well.
    """


    # list for storing results
    self.sols_varying_p0 = []

    # initial pressure
    pres_init = min(2.0, float(np.max(self.pres_array)))
    mass_init = 0.0

    if self.tidal:
        y_init = 2.0

    # central-pressure grid, shared by all EOS draws
    low_pres = max(1e-3, float(np.min(self.pres_array)))
    x = np.geomspace(1e-8, 2.5, 50)
    pres_space = np.geomspace(low_pres, pres_init, len(x))
    self.pres_space = pres_space

    # one column of results per EOS draw (single files give one draw)
    samples = 1 if self.pres_array.ndim == 1 else self.pres_array.shape[1]
    total_mass = np.zeros((len(x), samples))
    total_radius = np.zeros((len(x), samples))
    total_pres_central = np.zeros((len(x), samples))
    yR_all = np.zeros((len(x), samples))
    if self.tidal:
        tidal_all = np.zeros((len(x), samples))
        k2_all = np.zeros((len(x), samples))
    max_mass_arr = np.zeros(samples)
    max_radius_arr = np.zeros(samples)
    max_pres_arr = np.zeros(samples)

    # loop over the EOS draws
    for j in range(samples):
        self._set_eos_interpolants(draw=j)

        # per-draw arrays (impermanent so this will work)
        max_mass = np.zeros(len(x))
        pres_central = np.zeros(len(x))
        max_radius = np.zeros(len(x))
        yR = np.zeros(len(x))

        # loop over the TOV equations
        for i in range(len(x)):
            # initial conditions
            init_guess = [pres_space[i], mass_init]
            if self.tidal:
                init_guess.append(y_init)

            # high fidelity (four function evals per sol_pt)
            if self.solver == "RK4":
                xval, sol = self.RK4(
                    self.tov_equations_scaled, init_guess, 1e-3, 4.0,
                    self.sol_pts
                )

            # low fidelity (two function evals per sol_pt)
            elif self.solver == "RK2":
                xval, sol = self.RK2(
                    self.tov_equations_scaled, init_guess, 1e-3, 4.0,
                    self.sol_pts
                )

            # low fidelity (one function eval per sol_pt)
            elif self.solver == "euler":
                xval, sol = self.euler(
                    self.tov_equations_scaled, init_guess, 1e-3, 4.0,
                    self.sol_pts
                )

            # adaptive (typically high fidelity)
            elif self.solver == "solve_ivp":
                if self.solve_ivp_kwargs is None:
                    self.solve_ivp_kwargs = {
                        "method": "RK45",
                        "atol": 5e-14,
                        "rtol": self.tol,
                        "max_step": 0.01,
                        "dense_output": True,
                    }
                span = [1e-3, 2.5] if self.tidal else [1e-8, 2.5]
                result = solve_ivp(
                    self.tov_equations_scaled,
                    span,
                    init_guess,
                    **self.solve_ivp_kwargs,
                )
                if not result.success:
                    print("Solver failed.")
                    print(result.message)
                xval = result.t
                sol = result.y
            else:
                raise ValueError(
                    f'Solver, {self.solver} unknown. Must be "RK4", "RK2", '
                    f'"euler", or "solve_ivp".'
                )

            # index of the stellar surface: last point with positive pressure
            pressure_positive_indices = np.flatnonzero(sol[0] > 1e-10)
            if pressure_positive_indices.size > 0:
                index_mass = pressure_positive_indices[-1]
                max_mass[i] = sol[1, index_mass]
            else:
                # pressure never stayed positive; treat the center as surface
                index_mass = 0
                max_mass[i] = 0.0

            # central pressure
            pres_central[i] = np.max(sol[0])

            # maximum radius
            max_radius[i] = xval[index_mass]

            if self.tidal:
                yR[i] = sol[2][index_mass]

            # collect the results for each central pressure
            columns = [
                xval[: index_mass + 1] * self.rad0,  # radius
                sol[0][: index_mass + 1] * self.pres0,  # pressure
                sol[1][: index_mass + 1] * self.mass0,  # mass
            ]
            if self.tidal:
                columns.append(sol[2][: index_mass + 1])
            self.sols_varying_p0.append(np.column_stack(columns).T)

        # scale results and store this draw's column
        max_mass = max_mass * self.mass0
        max_radius = max_radius * self.rad0
        pres_central = pres_central * self.pres0
        total_mass[:, j] = max_mass
        total_radius[:, j] = max_radius
        total_pres_central[:, j] = pres_central
        yR_all[:, j] = yR

        # max mass calculation, radius, and central pressure
        corr_radius_index = int(np.argmax(max_mass))
        max_mass_arr[j] = max_mass[corr_radius_index]
        max_radius_arr[j] = max_radius[corr_radius_index]
        max_pres_arr[j] = pres_central[corr_radius_index]

        print(
            "Max mass: ",
            max_mass_arr[j],
            "Radius: ",
            max_radius_arr[j],
            "Central pressure: ",
            max_pres_arr[j],
        )

        # tidal deformability for this draw
        if self.tidal:
            tidal_all[:, j], k2_all[:, j] = self.tidal_def(
                yR, max_mass, max_radius
            )

    # single-draw results stay 1-D / scalar for backward compatibility
    single = samples == 1
    self.total_mass = total_mass[:, 0] if single else total_mass
    self.total_radius = total_radius[:, 0] if single else total_radius
    self.total_pres_central = (
        total_pres_central[:, 0] if single else total_pres_central
    )
    self.yR = yR_all[:, 0] if single else yR_all
    if self.tidal:
        self.tidal_deformability = tidal_all[:, 0] if single else tidal_all
        self.k2 = k2_all[:, 0] if single else k2_all

    self.maximum_mass = float(max_mass_arr[0]) if single else max_mass_arr
    self.corr_radius = float(max_radius_arr[0]) if single else max_radius_arr
    self.corr_pres = float(max_pres_arr[0]) if single else max_pres_arr

    # save these results
    self.max_mass_arr = self.maximum_mass
    self.max_radius_arr = self.corr_radius
    self.max_pres_arr = self.corr_pres

    if verbose:
        print("Max mass array: ", max_mass)

        # plot stuff (last draw for multi-draw input)
        plt.plot(max_radius, max_mass, label=r"TOV")
        plt.xlabel("Radius [km]")
        plt.ylabel("Max Mass [M_solar]")
        plt.legend()
        plt.show()
        plt.plot(max_radius, pres_central)
        plt.xlabel("Radius [km]")
        plt.ylabel("Central pressure [MeV/fm^3]")
        plt.show()

        if self.tidal:
            plt.plot(max_mass / max_radius, yR)
            plt.xlabel(r"$\beta$")
            plt.ylabel(r"y(r)")
            plt.savefig("yr_scaled.png")
            plt.show()
            plt.plot(max_radius, tidal_all[:, -1], color="k")
            plt.xlabel("Radius", fontsize=14)
            plt.ylabel(r"$\Lambda(R)$", fontsize=14)
            plt.xticks(fontsize=12)
            plt.yticks(fontsize=12)
            plt.savefig("tidal.png")
            plt.show()
            plt.plot(max_mass, tidal_all[:, -1], color="k")
            plt.xlabel(r"Mass [$M_{\odot}$]", fontsize=14)
            plt.ylabel(r"$\Lambda(M)$", fontsize=14)
            plt.xticks(fontsize=12)
            plt.yticks(fontsize=12)
            plt.savefig("tidal_mass.png")
            plt.show()
            plt.plot(max_mass / max_radius, k2_all[:, -1], color="k")
            plt.xticks(fontsize=12)
            plt.yticks(fontsize=12)
            plt.xlabel(r"$\beta$", fontsize=14)
            plt.ylabel(r"$k_{2}(\beta)$", fontsize=14)
            plt.savefig("k2.png")
            plt.show()
            plt.plot(max_radius, k2_all[:, -1], color="k")
            plt.xticks(fontsize=12)
            plt.yticks(fontsize=12)
            plt.xlabel(r"$R$ [km]", fontsize=14)
            plt.ylabel(r"$k_{2}(R)$", fontsize=14)
            plt.xlim(9.0, 15.0)
            plt.savefig("k2R.png")
            plt.show()

        # check the solution
        print(
            "Radius: ",
            self.corr_radius,
            "Maximum mass: ",
            self.maximum_mass,
            "Central pressure: ",
            self.corr_pres,
        )

    # if desired, write to a file
    if write_to_file:
        tov_data = np.column_stack(
            [self.total_radius, self.total_mass, self.total_pres_central]
        )
        out_dir = Path("TOV_data")
        out_dir.mkdir(parents=True, exist_ok=True)
        file_name = out_dir / f"rpm_results_{self.eos_name}.txt"
        header = "Radius[km] Mass[Msol] Central_Pressure[MeV/fm3]"
        np.savetxt(file_name, tov_data, header=header, delimiter=" ")

    if self.tidal:
        return (
            self.total_radius,
            self.total_pres_central,
            self.total_mass,
            self.k2,
            self.tidal_deformability,
        )

    return self.total_radius, self.total_pres_central, self.total_mass

max_arrays()

Returns the max arrays needed for the interval calculation.

Returns:

Type Description

self.max_radius_arr (array): Maximum radius array.

self.max_pres_arr (array): Maximum central pressure array.

self.max_mass_arr (array): Maximum mass array.
Source code in src/slmemulator/TOV_class.py 
847
848
849
850
851
852
853
854
855
856
857
858
859
def max_arrays(self):
    r"""
    Returns the max arrays needed for the interval calculation.

    Parameters:
        None.

    Returns:
        self.max_radius_arr (array): Maximum radius array.
        self.max_pres_arr (array): Maximum central pressure array.
        self.max_mass_arr (array): Maximum mass array.
    """
    return self.max_radius_arr, self.max_pres_arr, self.max_mass_arr

central_dens(pres_arr=None)

Calculation to determine the central density of the star at the maximum mass and radius determined from the tov_routine().

Example

tov.central_dens()

Parameters:

Name Type Description Default
pres_arr array

An optional pressure array to use for calculating central densities at places other than the absolute TOV maximum mass of each curve. Default is None, and code will use absolute TOV maximum mass central pressure class array.

 None

Returns:

Name Type Description
c_dens float or array

The central density, one value per EOS draw.
Source code in src/slmemulator/TOV_class.py 
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
def central_dens(self, pres_arr=None):
    r"""
    Calculation to determine the central density of the star
    at the maximum mass and radius determined from the tov_routine().

    Example:
        tov.central_dens()

    Parameters:
        pres_arr (array): An optional pressure array to use for
            calculating central densities at places other
            than the absolute TOV maximum mass of each curve.
            Default is None, and code will use absolute TOV
            maximum mass central pressure class array.

    Returns:
        c_dens (float or array): The central density, one value per
            EOS draw.
    """

    # one column per EOS draw [MeV/fm^3]
    press = self.pres_array_unscaled
    single = press.ndim == 1
    if single:
        press = press[:, None]
    nB = self.nB_array if self.nB_array.ndim > 1 else self.nB_array[:, None]

    _, samples = press.shape
    target = self.corr_pres if pres_arr is None else pres_arr
    target = np.broadcast_to(np.atleast_1d(target), (samples,))

    c_dens = np.zeros(samples)
    for j in range(samples):
        # interpolate the EOS to find the proper central densities
        p_n_interp = interp1d(
            press[:, j],
            nB[:, min(j, nB.shape[1] - 1)],
            kind="cubic",
            fill_value="extrapolate",
        )
        # solve at the proper central pressure for nB_central
        c_dens[j] = p_n_interp(target[j])

    return c_dens[0] if single else c_dens

canonical_NS_radius()

Calculation of the radius of a 1.4 M_sol neutron star.

Example

tov.canonical_NS_radius()

Returns:

Name Type Description
rad_14 float or array

The 1.4 M_sol radius, one value per EOS draw.
Source code in src/slmemulator/TOV_class.py 
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
def canonical_NS_radius(self):
    r"""
    Calculation of the radius of a 1.4 M_sol neutron star.

    Example:
        tov.canonical_NS_radius()

    Parameters:
        None.

    Returns:
        rad_14 (float or array): The 1.4 M_sol radius, one value per
            EOS draw.
    """

    mass, radius = self.total_mass, self.total_radius
    if mass.ndim == 1:
        m_r_interp = interp1d(
            mass, radius, kind="linear", fill_value="extrapolate"
        )
        return m_r_interp(1.4)

    # one column per EOS draw
    samples = mass.shape[1]
    rad_14 = np.zeros(samples)
    for j in range(samples):
        m_r_interp = interp1d(
            mass[:, j],
            radius[:, j],
            kind="linear",
            fill_value="extrapolate",
        )
        rad_14[j] = m_r_interp(1.4)

    return rad_14

tovScaledRev

Information about the code: This code solves TOV equations for mass radius relations. This can also plot the mass-radius curve.

The code solves dr/dp and dm/dp instead of the regular way.

USE: To use the code, here are the steps: 1) Include the file in your main code e.g. import tov_class as tc 2) Load the EoS using the ToV loader, tc.ToV(filename, arraysize) 3) call the solver as tc.ToV.mass_radius(min_pressure, max_pressure) 4) To plot, follow the code in main() on creating the dictionary of inputs

Updates: Version 0.0.1-1 Solves ToV, can only take inputs of pressure (MeV/fm^3), energy density in MeV, baryon density in fm^-3 in ascending order.

TOV(filename, imax, tidal=False)

Solves TOV equations and gives data-table, mass-radius plot and max. mass, central pressure and central density by loading an EoS datafile.

Source code in src/slmemulator/tovScaledRev.py 
42
43
44
45
46
47
48
49
50
51
52
def __init__(self, filename, imax, tidal=False):
    self.file = np.loadtxt(filename)
    self.tidal = tidal
    self.e_in = self.file[:, 1] / eps0  # Scaled Energy density
    self.p_in = self.file[:, 2] / pres0  # Scaled pressure
    self.nb_in = self.file[:, 0]  # Scaled baryon density
    if self.tidal:
        self.cs2_in = self.file[:, 3]  # sound speed
    self.imax = imax
    self.radius = np.empty(self.imax)
    self.mass = np.empty(self.imax)

pressure_from_nb(nb)

Evaluate pressure from number density using interpolation

Source code in src/slmemulator/tovScaledRev.py 
54
55
56
57
58
59
def pressure_from_nb(self, nb):
    """Evaluate pressure from number density using interpolation"""
    p1 = interp1d(
        self.nb_in, self.p_in, axis=0, kind="linear", fill_value="extrapolate"
    )
    return p1(nb)

energy_from_pressure(pressure)

Evaluate energy density from pressure using interpolation

Source code in src/slmemulator/tovScaledRev.py 
61
62
63
64
65
66
67
68
69
70
def energy_from_pressure(self, pressure):
    """Evaluate energy density from pressure using interpolation"""
    plow = 1e-10 / pres0
    if pressure < plow:
        return 2.6e-310
    else:
        e1 = interp1d(
            self.p_in, self.e_in, axis=0, kind="linear", fill_value="extrapolate"
        )
        return e1(pressure)

pressure_from_energy(energy)

Evaluate pressure from energy density using interpolation

Source code in src/slmemulator/tovScaledRev.py 
72
73
74
75
76
77
def pressure_from_energy(self, energy):
    """Evaluate pressure from energy density using interpolation"""
    p1 = interp1d(
        self.e_in, self.p_in, axis=0, kind="linear", fill_value="extrapolate"
    )
    return p1(energy)

baryon_from_energy(energy)

Evaluate number density from energy using interpolation

Source code in src/slmemulator/tovScaledRev.py 
79
80
81
82
83
84
def baryon_from_energy(self, energy):
    """Evaluate number density from energy using interpolation"""
    n1 = interp1d(
        self.e_in, self.nb_in, axis=0, kind="linear", fill_value="extrapolate"
    )
    return n1(energy)