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 nonlinear terms.
Source code in src/slmemulator/SLM.py 
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
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 nonlinear terms.
    """
    n_variables, n_snapshots = X.shape

    # Calculate the number of augmented terms
    num_augmented_terms = sum(
        1 for _ in combinations_with_replacement(range(n_variables), 2)
    )

    # Pre-allocate the augmented matrix
    augmented_X = np.empty(
        (n_variables + num_augmented_terms, n_snapshots), dtype=X.dtype
    )
    augmented_X[:n_variables, :] = X

    # Add quadratic terms and cross-products
    row_idx = n_variables
    for i, j in combinations_with_replacement(range(n_variables), 2):
        augmented_X[row_idx, :] = X[i, :] * X[j, :]
        row_idx += 1

    return augmented_X

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

Dynamic Mode decomposition for 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, it will go up to the maximum possible rank (min(X.shape)).

 None
Source code in src/slmemulator/SLM.py 
 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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def SLM(X, dt, error_threshold=1e-4, max_r=None):
    r"""
    Dynamic Mode decomposition for 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 (np.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, it will go up to
            the maximum possible rank (min(X.shape)).
    """
    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)

    # Determine the maximum possible rank
    max_possible_r = min(X1.shape)
    if max_r is None:
        max_r_to_check = max_possible_r
    else:
        max_r_to_check = min(max_r, max_possible_r)

    # Initialize r and best_error
    r_optimal = 1
    min_error = float("inf")
    best_Xdmd = None
    best_Phi = None
    best_omega = None
    best_lambda_vals = None
    best_b = None

    # Iterate through possible ranks to find the optimal 'r'
    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, :]

        # Compute Atilde
        Atilde = U_r.T @ X2 @ V_r.T @ S_r_inv

        # Compute eigenvectors and eigenvalues
        D, W_r = np.linalg.eig(Atilde)

        Phi_current = X2 @ V_r.T @ S_r_inv @ W_r  # DMD modes
        lambda_vals_current = D  # discrete-time eigenvalues
        omega_current = np.log(lambda_vals_current) / dt  # continuous-time eigenvalue

        # Compute DMD mode amplitudes b
        x1 = X1[:, 0]
        b_current = np.linalg.lstsq(Phi_current, x1, rcond=None)[0]

        # DMD reconstruction for the current 'r'
        mm1 = X1.shape[1] + 1
        t = np.arange(mm1) * dt

        time_dynamics_current = b_current[:, np.newaxis] * np.exp(
            omega_current[:, np.newaxis] * t
        )
        Xdmd2_current = Phi_current @ time_dynamics_current

        # Truncate to original number of variables (log-transformed)
        Xdmd_current_original_vars = Xdmd2_current[:n, :]

        # Calculate error (max absolute difference)
        # Using original X (log-transformed) for comparison
        current_error = np.max(np.abs(X - Xdmd_current_original_vars))

        # print(f"Testing r={r_current}: Max absolute error = {current_error:.6f}")

        if current_error <= error_threshold:
            r_optimal = r_current
            min_error = current_error
            # Store the results for the optimal r
            best_Xdmd = Xdmd_current_original_vars
            best_Phi = Phi_current[:n, :]  # Truncate Phi to original variables
            best_omega = omega_current
            best_lambda_vals = lambda_vals_current
            best_b = b_current
            break  # Found the smallest r that satisfies the threshold

        # If we didn't meet the threshold, but this r gives the best error so far, keep its results
        if current_error < min_error:
            min_error = current_error
            r_optimal = r_current
            best_Xdmd = Xdmd_current_original_vars
            best_Phi = Phi_current[:n, :]
            best_omega = omega_current
            best_lambda_vals = lambda_vals_current
            best_b = b_current

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

    # If no 'r' met the threshold, use the one that gave the minimum error
    if best_Xdmd is None:  # This should not happen if max_r_to_check >= 1
        # Fallback to a default if no optimal r is found (e.g., r=1)
        r_optimal = 1
        U_r = U_full[:, :r_optimal]
        S_r_inv = np.diag(1.0 / S_full[:r_optimal])
        V_r = Vt_full[:r_optimal, :]
        Atilde = U_r.T @ X2 @ V_r.T @ S_r_inv
        D, W_r = np.linalg.eig(Atilde)
        best_Phi = X2 @ V_r.T @ S_r_inv @ W_r
        best_lambda_vals = D
        best_omega = np.log(best_lambda_vals) / dt
        x1 = X1[:, 0]
        best_b = np.linalg.lstsq(best_Phi, x1, rcond=None)[0]
        mm1 = X1.shape[1] + 1
        t = np.arange(mm1) * dt
        time_dynamics = best_b[:, np.newaxis] * np.exp(best_omega[:, np.newaxis] * t)
        best_Xdmd = (best_Phi @ time_dynamics)[:n, :]
        print(
            f"Warning: No r met the threshold. Using r={r_optimal} with max error {np.max(np.abs(X - best_Xdmd)):.6f}"
        )

    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, mass and central pressure

Parameters:

Name Type Description Default
fileName str

Filename containing the EOS in the format nb (fm^-3), E (MeV), P (MeV/fm^3)

 required

Returns:

Name Type Description
dataArray array

Data array containing radii, central pressure and mass (includes tidal deformability k_2 if set to true).
Source code in src/slmemulator/SLM.py 
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
229
230
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
def solve_tov(fileName, tidal=False, parametric=False, mseos=True):
    r"""
    Solves the TOV equation and returns radius, mass and central pressure

    Parameters:
        fileName (str): Filename containing the EOS in the format nb (fm^-3),
            E (MeV), P (MeV/fm^3)

    Returns:
        dataArray (array): Data array containing radii, central pressure
            and mass (includes tidal deformability k_2 if set to true).
    """
    # Get current paths from the config module
    paths = get_paths()

    eos_file_path = None
    tov_path_target = None

    if parametric is False:
        # Access EOS_Data as package resource
        try:
            slm_package_base = pkg_resources.files(_SLM_PACKAGE_NAME)
            # Get a Path-like object to the EOS_Data file within the installed package
            eos_file_path = pkg_resources.files(_INTERNAL_EOS_DATA_PATH).joinpath(
                fileName
            )
            eos_data_resource_dir = slm_package_base.joinpath("EOS_Data")
            eos_file_path = eos_data_resource_dir.joinpath(fileName)
            # When passing to TOV, ensure it's a string path if the solver expects it
            file = TOV(str(eos_file_path), tidal=tidal)
            tov_path_target = paths["tov_data_dir"]
        except FileNotFoundError:
            raise FileNotFoundError(
                f"Internal EOS file '{fileName}' not found in package data."
            )
    else:
        # Access EOS_files from external configured path
        if mseos is True:
            eos_file_path = paths["mseos_path"] / fileName
            tov_path_target = paths["mseos_tov_path"]
        else:
            eos_file_path = paths["qeos_path"] / fileName
            tov_path_target = paths["qeos_tov_path"]
            print("Path:", eos_file_path)

        file = TOV(str(eos_file_path), tidal=tidal)  # Pass as string path

    # Now tov_path_target is guaranteed to be assigned
    if not tov_path_target.exists():
        tov_path_target.mkdir(parents=True, exist_ok=True)

    file.tov_routine(verbose=False, write_to_file=False)
    print("R of 1.4 solar mass star: ", file.canonical_NS_radius())
    dataArray = [
        file.total_radius.flatten(),
        file.total_pres_central.flatten(),
        file.total_mass.flatten(),
    ]
    if tidal is True:
        dataArray.append(file.k2.flatten())
        # dataArray.append(file.tidal_deformability.flatten()[::-1])

    dataArray = np.asarray(dataArray, dtype=np.float64)

    # Construct the output filename more robustly
    name_parts = os.path.basename(fileName).strip(".txt").split("_")
    print("Name parts:", name_parts)
    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"

    # Save directly to the target path without changing directory
    output_full_path = tov_path_target / output_file_name
    np.savetxt(output_full_path, dataArray.T, fmt="%1.8e")
    return dataArray

cleanData

clean_directory(directory: Optional[str] = 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 
 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
def clean_directory(directory: Optional[str] = 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

get_paths(output_base_dir: Optional[Path] = 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 (Sparse Linear Modeling) 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 current working directory is implicitly used as the root for relative paths. 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 typically include 'pSLM'; otherwise, they'll use 'SLM'. Defaults to True.

 True
include_slm_paths bool

If True, the dictionary will include the specific directories (current_results_dir, current_plots_dir, current_test_dir) related to the SLM/pSLM outputs. Defaults to True.

 True

Returns:

Type Description
Dict[str, Path]

dict[str, pathlib.Path]: A dictionary containing all relevant path configurations. Keys include:
Source code in src/slmemulator/config.py 
 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
170
171
172
173
174
175
176
177
178
179
180
181
def get_paths(
    output_base_dir: Optional[Path] = 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 (Sparse Linear Modeling) 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 current working directory 
            is implicitly used as the root for relative paths. 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 typically include 
            'pSLM'; otherwise, they'll use 'SLM'. Defaults to ``True``.
        include_slm_paths (bool, optional):
            If ``True``, the dictionary will include the specific directories 
            (``current_results_dir``, ``current_plots_dir``, ``current_test_dir``) 
            related to the SLM/pSLM outputs. Defaults to ``True``.

    Returns:
        dict[str, pathlib.Path]: A dictionary containing all relevant path configurations. Keys include:
    """

    project_root = PROJECT_ROOT
    output_data_base = output_base_dir or project_root
    src_dir = project_root / DEFAULT_SRC_SUBDIR_NAME

    # Define extensions to strip
    EXTENSIONS_TO_STRIP = [".txt", ".dat", ".table"]

    # Strip extensions and convert to uppercase 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 # Stop after stripping the first matching extension

    # --- Dynamic EOS-specific folder name, ensuring it's uppercase for directories ---
    EOS_FOLDER_NAME = clean_eos_name.upper()

    # --- Define SLM subdirectory based on run type (SLM or pSLM) ---
    SLM_SUBDIR = DEFAULT_PSLM_SUBDIR_NAME if is_parametric_run else DEFAULT_SLM_SUBDIR_NAME

    paths = {
        "project_root": project_root,
        "src_dir": src_dir,
        # Internal package resources (read-only)
        "package_eos_codes_dir": pkg_resources.files("slmemulator").joinpath(
            DEFAULT_EOS_CODES_SUBDIR_NAME
        ),
        "package_eos_data_dir": pkg_resources.files("slmemulator").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_path": output_data_base / DEFAULT_TRAIN_DATA_SUBDIR_NAME,
        "val_path": 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,

        # 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,
    }

    # Conditionally add SLM-specific result/plot/test subdirectories
    if include_slm_paths:

        # Define the dynamic current SLM results and plots 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
        )

        # --- NEW PATH: Dynamic SLM Test Data Path ---
        if is_parametric_run:
            # Only create the /pSLM folder for tests if it's a parametric run
            paths["current_slm_tests_dir"] = (
                output_data_base
                / DEFAULT_TEST_DATA_SUBDIR_NAME
                / EOS_FOLDER_NAME
                / DEFAULT_PSLM_SUBDIR_NAME  # Always use 'pSLM' for parametric tests
            )
        else:
            # If not parametric, use the base testData/EOS_NAME folder
            paths["current_slm_tests_dir"] = (
                output_data_base
                / DEFAULT_TEST_DATA_SUBDIR_NAME
                / EOS_FOLDER_NAME
            )

    # Ensure all paths are Path objects and filter out any None/non-Path values
    final_paths = {}
    for k, v in paths.items():
        if isinstance(v, Path):
            final_paths[k] = v
        elif include_slm_paths and k.startswith("current_slm_") and k in paths:
             final_paths[k] = paths[k]

    return final_paths

create_necessary_dirs(paths: Dict[str, Path], additional_dirs: Optional[List[Path]] = None) -> None

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

This function iterates through all Path objects provided in the input dictionary's values and the optional list, ensuring that each directory is created if it does not already exist. It uses pathlib.Path.mkdir with parents=True and exist_ok=True, meaning it will create parent directories if necessary and will not raise an error if the directory already exists.

Parameters:

Name Type Description Default
paths dict[str, Path]

A dictionary where keys are string identifiers (e.g., 'output_path') and values are :class:pathlib.Path objects representing the directories to be created.

 required
additional_dirs list[Path]

An optional list of additional :class:pathlib.Path objects representing directories to be created (e.g., user-managed data or model directories). The default is 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 
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
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
def create_necessary_dirs(
    paths: Dict[str, Path], 
    additional_dirs: Optional[List[Path]] = None
) -> None:
    """
    Creates necessary directories specified in a dictionary and an optional list.

    This function iterates through all Path objects provided in the input dictionary's
    values and the optional list, ensuring that each directory is created if it
    does not already exist. It uses pathlib.Path.mkdir with parents=True and
    exist_ok=True, meaning it will create parent directories if necessary and
    will not raise an error if the directory already exists.

    Parameters:
        paths (dict[str, pathlib.Path]): A dictionary where keys are string identifiers (e.g., 'output_path')
            and values are :class:`pathlib.Path` objects representing the directories
            to be created.
        additional_dirs (list[pathlib.Path], optional): An optional list of additional :class:`pathlib.Path` objects representing
            directories to be created (e.g., user-managed data or model directories). 
            The default is None.

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

    """

    # Define a list of keys corresponding to directories that should be created.
    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_path",
        "generated_eos_files_dir",

        # The 'current' paths are the effective targets, so ensure they are created
        "current_eos_input_dir",
        "current_tov_data_dir",
        "current_slm_results_dir",
        "current_slm_plots_dir",
        "current_slm_tests_dir", # <-- ADDED NEW PATH

        # Kept for backward compatibility
        "qeos_path_specific",
        "mseos_path_specific",
        "qeos_tov_path_specific",
        "mseos_tov_path_specific",
    ]

    unique_dirs_to_create = set()
    for key in output_directory_keys:
        if key in paths and isinstance(paths[key], Path):
            unique_dirs_to_create.add(paths[key])

    # Add any explicitly provided additional directories
    if additional_dirs:
        for path_obj in additional_dirs:
            if isinstance(path_obj, Path):
                unique_dirs_to_create.add(path_obj)

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

pSLM

ParametricSLM(fileList, filePath, tidal=False, error_threshold=0.0001, max_r=None)

Source code in src/slmemulator/pSLM.py 
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
def __init__(
    self, fileList, filePath, tidal=False, error_threshold=1e-4, max_r=None
):
    # fileList is expected to be a list of Path objects from testpSLM.py
    self.fileList = [Path(f) for f in fileList]  # Ensure all are Path objects
    self.filePath = Path(filePath)  # Ensure filePath is a Path object
    self.tidal = tidal
    self.error_threshold = error_threshold
    self.max_r = max_r
    self.n = None  # Number of quantities/rows in the data X
    self.mm1 = (
        None  # Number of time points/columns in the data X (m-1 snapshots for DMD)
    )

    # Lists to store DMD components for each training file
    self.Phi_list = []
    self.omega_list = []
    self.D_list = []
    self.b_list = []
    self.train_data_shapes = []  # Store shapes (n, mm1) of training data

    # Store extracted parameters for each file, in the same order as fileList
    self.file_params = []
    self.param_scaler = None  # Scaler for parameters, fitted in .fit()
    self.scaled_params = None  # Scaled parameters of training data

augment_data(X) staticmethod

Augments the input data X by adding quadratic terms (X_i * X_j).

Parameters:

Name Type Description Default
X ndarray

Original data matrix (n, m).

 required

Returns:

Type Description

np.ndarray (np.ndarray): Augmented data matrix.
Source code in src/slmemulator/pSLM.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
@staticmethod
def augment_data(X):
    """
    Augments the input data X by adding quadratic terms (X_i * X_j).

    Parameters:
        X (np.ndarray): Original data matrix (n, m).

    Returns:
        np.ndarray (np.ndarray): Augmented data matrix.
    """
    n, m = X.shape
    out = [X]
    # Add quadratic terms (combinations with replacement of rows of X)
    for i, j in combinations_with_replacement(range(n), 2):
        out.append(X[i] * X[j])
    # Add cubic terms if needed (example commented out)
    # for i, j, k in combinations_with_replacement(range(n), 3):
    #     out.append(X[i] * X[j] * X[k])
    return np.vstack(out)

fit()

Fits the Parametric SLM model by processing each file in fileList, performing DMD, and storing the DMD components along with extracted parameters.

Source code in src/slmemulator/pSLM.py 
268
269
270
271
272
273
274
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
317
318
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
383
384
385
386
387
def fit(self):
    """
    Fits the Parametric SLM model by processing each file in fileList,
    performing DMD, and storing the DMD components along with extracted parameters.
    """
    self.Phi_list = []
    self.omega_list = []
    self.D_list = []
    self.b_list = []
    self.file_params = []  # Reset for clarity, though init should handle it
    self.train_data_shapes = []  # Store shapes for consistency check in predict

    processed_files_count = 0
    for file_path_obj in self.fileList:  # Iterate through Path objects
        # Validate file_path_obj is a Path object and exists
        if not isinstance(file_path_obj, Path) or not file_path_obj.exists():
            print(f"Skipping invalid or non-existent file: {file_path_obj}")
            continue

        # Extract parameters from the filename (e.g., TOV_MS_Ls_Lv_zeta_xi.txt)
        # We assume the file name format is consistent with _extract_params_from_tov_filename
        # We determine mseos based on 'MS' or 'QEOS' in the filename
        is_mseos = "MS" in file_path_obj.stem  # Check if filename contains 'MS'
        # Call the static method using self.
        params = self._extract_params_from_tov_filename(
            file_path_obj.stem, is_mseos
        )

        if params is None:  # _extract_params_from_tov_filename prints warning
            continue  # Skip if parameters could not be extracted

        try:
            # Load data from the file
            # The file is a Path object, np.loadtxt can directly handle it
            data = np.loadtxt(file_path_obj)

            if not np.all(np.isfinite(data)):
                print(
                    f"Skipping {file_path_obj.name}: contains NaN or Inf values. Please check data quality."
                )
                continue

            # Process data (assuming data has columns like Radius, Mass, etc.)
            # X should be (num_quantities, num_points)
            epsilon = 1e-9  # Small constant to prevent log(0)
            if self.tidal:
                # Assuming 4 columns: Radius, Mass, Central Pressure, Tidal Deformability
                # Ensure data has at least 4 columns
                if data.shape[1] < 4:
                    print(
                        f"Skipping {file_path_obj.name}: Expected at least 4 columns for tidal data, got {data.shape[1]}."
                    )
                    continue
                X = np.log(data[:, [0, 1, 2, 3]].T + epsilon)
            else:
                # Assuming 3 columns: Radius, Mass, Central Pressure (or other relevant)
                # Ensure data has at least 3 columns
                if data.shape[1] < 3:
                    print(
                        f"Skipping {file_path_obj.name}: Expected at least 3 columns for non-tidal data, got {data.shape[1]}."
                    )
                    continue
                X = np.log(data[:, [0, 1, 2]].T + epsilon)

            X = np.asarray(
                X, dtype=np.float64
            )  # Ensure float64 for numerical stability

            # Store shape for consistency checks in predict
            if (
                self.n is None
            ):  # Only set for the very first successfully processed file
                self.n, self.mm1 = X.shape[0], X.shape[1]
                self.augmented_n = self.augment_data(
                    np.zeros((self.n, self.mm1))
                ).shape[0]
            elif X.shape[0] != self.n or X.shape[1] != self.mm1:
                print(
                    f"Warning: {file_path_obj.name} has inconsistent shape {X.shape}. Expected ({self.n}, {self.mm1}). Skipping."
                )
                continue  # Skip files with inconsistent shapes for DMD

            dt = 1.0  # Assuming uniform time step of 1.0 between data points
            result = self._SLM_auto_r(X, dt)  # Perform SLM and get optimal rank

            # Check if SLM_auto_r returned a valid result
            if result["err"] == float("inf"):
                print(
                    f"Failed to fit SLM model for {file_path_obj.name}. Skipping."
                )
                continue

            # Store DMD components and parameters for successful fits
            self.Phi_list.append(result["Phi"])
            self.omega_list.append(result["omega"])
            self.D_list.append(result["D"])
            self.b_list.append(result["b"])
            self.file_params.append(
                params
            )  # Store parameters for the processed file
            self.train_data_shapes.append(
                X.shape
            )  # Store shape for this specific training data

            processed_files_count += 1
            # print(f"Successfully processed {file_path_obj.name}. Max error: {result['err']:.4e}")

        except Exception as e:
            print(f"Error processing file {file_path_obj.name}: {e}. Skipping.")
            continue

    if processed_files_count == 0:
        raise RuntimeError(
            "No valid training files processed! Please check your input data and `error_threshold`."
        )

    # Fit StandardScaler on the collected parameters
    self.params = np.array(self.file_params)
    self.param_scaler = StandardScaler().fit(self.params)
    self.scaled_params = self.param_scaler.transform(self.params)

predict(param_values, k=1, output_interp=False, distance_threshold=None)

Predicts the dynamics (Xdmd) for a given set of parameters by finding nearest neighbors in the parameter space and averaging their DMD components or using the closest one.

Parameters:

Name Type Description Default
param_values list or ndarray

A list or array of parameter values to predict for (e.g., [Ls, Lv, zeta, xi] for MSEOS). This argument is mandatory.

 required
k int

Number of nearest neighbors to use for prediction. Defaults to 1 (pure nearest neighbor).

 1
output_interp bool

If True, interpolate the final Xdmd outputs by averaging curves. If False (default) and k>1, averages the DMD components (Phi, omega, D, b).

 False
distance_threshold float

If the distance to the closest neighbor exceeds this threshold, a warning is printed.

 None

Returns:

Name Type Description
tuple

(Phi_avg/Phi_nn, omega_avg/omega_nn, D_avg/D_nn, b_avg/b_nn, Xdmd_predicted, t) - Reconstructed data (Xdmd_predicted) and related averaged DMD components.

Raises: ValueError: If the model has not been fitted or no training data is available.

Source code in src/slmemulator/pSLM.py 
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
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
436
437
438
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
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
def predict(self, param_values, k=1, output_interp=False, distance_threshold=None):
    """
    Predicts the dynamics (Xdmd) for a given set of parameters by finding nearest neighbors
    in the parameter space and averaging their DMD components or using the closest one.

    Parameters:
        param_values (list or np.ndarray): A list or array of parameter values to predict for
                                                     (e.g., [Ls, Lv, zeta, xi] for MSEOS). This argument is mandatory.
        k (int): Number of nearest neighbors to use for prediction.
                 Defaults to 1 (pure nearest neighbor).
        output_interp (bool): If True, interpolate the final Xdmd outputs by averaging curves.
                              If False (default) and k>1, averages the DMD components (Phi, omega, D, b).
        distance_threshold (float, optional): If the distance to the closest neighbor
                                              exceeds this threshold, a warning is printed.

    Returns:
        tuple: (Phi_avg/Phi_nn, omega_avg/omega_nn, D_avg/D_nn, b_avg/b_nn, Xdmd_predicted, t)
               - Reconstructed data (Xdmd_predicted) and related averaged DMD components.
    Raises:
        ValueError: If the model has not been fitted or no training data is available.
    """
    if not self.Phi_list:
        raise ValueError("Model not fitted. Call .fit() first.")
    if not self.file_params:
        raise ValueError("No training data parameters available from fit().")

    # If param_values are provided (the main use case for prediction)
    theta = np.asarray(param_values, dtype=np.float64)

    # Scale the input parameters using the fitted scaler
    try:
        scaled_theta = self.param_scaler.transform([theta])[0]
    except ValueError as e:
        raise ValueError(
            f"Parameter scaling failed for {theta}. Ensure input matches training dimensions. Error: {e}"
        )

    # Calculate Euclidean distances to all scaled training parameters
    dists = np.linalg.norm(self.scaled_params - scaled_theta, axis=1)

    # Get indices of k smallest distances
    # Ensure k does not exceed the number of available training samples
    k_actual = min(k, len(self.file_params))
    if k_actual == 0:
        raise ValueError("No training data available to make a prediction.")

    sorted_indices = np.argsort(dists)
    idxs = sorted_indices[:k_actual]

    # Check if the closest neighbor is too far (optional warning/handling)
    if distance_threshold is not None and dists[idxs[0]] > distance_threshold:
        print(
            f"Warning: Closest neighbor distance ({dists[idxs[0]]:.4f}) exceeds threshold ({distance_threshold}). "
            f"Prediction might be unreliable for parameters {param_values}."
        )

    # Use the shape of the training data corresponding to the first nearest neighbor
    # as the reference shape for reconstruction (assuming consistency)
    if not self.train_data_shapes:
        raise ValueError(
            "Training data shapes not recorded during fit(). Cannot predict."
        )

    # Use the original (non-augmented) shape of training data for reconstruction
    ref_n_original, ref_mm1_original = self.train_data_shapes[
        idxs[0]
    ]  # Use shape from first nearest neighbor

    dt = 1.0  # Assuming original dt=1.0

    if output_interp and k_actual > 1:
        # Option 1: Interpolate outputs by averaging reconstructed curves
        curves = []
        for idx in idxs:
            Phi = self.Phi_list[idx]
            omega = self.omega_list[idx]
            b = self.b_list[idx]

            t_reconstruction = np.arange(ref_mm1_original) * dt
            time_dynamics = b[:, np.newaxis] * np.exp(
                omega[:, np.newaxis] * t_reconstruction[np.newaxis, :]
            )
            Xdmd_full_aug = Phi @ time_dynamics
            Xdmd = Xdmd_full_aug[
                :ref_n_original, :
            ]  # Trim back to original dimensions
            curves.append(Xdmd)

        Xdmd_predicted = np.mean(curves, axis=0)  # Average the reconstructed curves

        # For consistency, return averaged DMD components even if output is interpolated
        Phi_avg = np.mean([self.Phi_list[i] for i in idxs], axis=0)
        omega_avg = np.mean([self.omega_list[i] for i in idxs], axis=0)
        D_avg = np.mean([self.D_list[i] for i in idxs], axis=0)
        b_avg = np.mean([self.b_list[i] for i in idxs], axis=0)

        return (Phi_avg, omega_avg, D_avg, b_avg, Xdmd_predicted, t_reconstruction)

    else:
        # Option 2: Average DMD components (Phi, omega, D, b) from k nearest neighbors
        # (or just use the single nearest neighbor if k=1 or output_interp is False)
        Phi_avg = np.mean([self.Phi_list[i] for i in idxs], axis=0)
        omega_avg = np.mean([self.omega_list[i] for i in idxs], axis=0)
        D_avg = np.mean([self.D_list[i] for i in idxs], axis=0)
        b_avg = np.mean([self.b_list[i] for i in idxs], axis=0)

        t_reconstruction = (
            np.arange(ref_mm1_original) * dt
        )  # Time points for prediction
        time_dynamics = b_avg[:, np.newaxis] * np.exp(
            omega_avg[:, np.newaxis] * t_reconstruction[np.newaxis, :]
        )
        Xdmd_full_aug = Phi_avg @ time_dynamics
        Xdmd_predicted = Xdmd_full_aug[
            :ref_n_original, :
        ]  # Trim back to original dimensions

        return (Phi_avg, omega_avg, D_avg, b_avg, Xdmd_predicted, t_reconstruction)

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 
519
520
521
522
523
524
525
526
527
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

gaussian_rbf(r, epsilon)

Gaussian Radial Basis Function.

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

recombination_thinning(A_matrix, y_vector, initial_solution_x=None, tolerance=1e-09)

Conceptual implementation of the recombination thinning process based on the paper. This function aims to find a sparse solution x' to the linear system Ax = y.

Parameters:

Name Type Description Default
A_matrix array

The matrix A in the linear system Ax = y. In the context of GRIM, this matrix relates features and linear functionals.

 required
y_vector array

The vector y in the linear system Ax = y. In the context of GRIM, this relates to the target function evaluated by linear functionals.

 required
initial_solution_x array

An initial solution vector x. If None, a least squares solution is computed.

 None
tolerance float

Tolerance for numerical stability, especially for SVD and checking non-zero components.

 1e-09

Returns:

Type Description

np.array: A sparser solution vector x' for the system Ax = y. Returns None if the system is inconsistent or other issues.
Source code in src/slmemulator/recombination.py 
14
15
16
17
18
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
def recombination_thinning(A_matrix, y_vector, initial_solution_x=None, tolerance=1e-9):
    """
    Conceptual implementation of the recombination thinning process based on the paper.
    This function aims to find a sparse solution x' to the linear system Ax = y.

    Parameters:
        A_matrix (np.array): The matrix A in the linear system Ax = y.
                             In the context of GRIM, this matrix relates features and
                             linear functionals.
        y_vector (np.array): The vector y in the linear system Ax = y.
                             In the context of GRIM, this relates to the target function
                             evaluated by linear functionals.
        initial_solution_x (np.array, optional): An initial solution vector x. If None,
                                                  a least squares solution is computed.
        tolerance (float): Tolerance for numerical stability, especially for SVD and
                           checking non-zero components.

    Returns:
        np.array: A sparser solution vector x' for the system Ax = y.
                  Returns None if the system is inconsistent or other issues.
    """

    if A_matrix.shape[0] != y_vector.shape[0]:
        raise ValueError("Dimensions of A_matrix and y_vector are not compatible.")

    # 2. Find an initial solution x if not provided
    if initial_solution_x is None:
        try:
            initial_solution_x = np.linalg.lstsq(A_matrix, y_vector, rcond=None)[0]
        except np.linalg.LinAlgError:
            print(
                "Warning: Least squares solution failed for initial_solution_x. Returning None."
            )
            return None

    # The `recombination_thinning` as provided aims to return an initial least squares solution.
    # For more advanced thinning/sparsity, additional iterative logic or L1 regularization
    # would typically be integrated here. As per your provided `bgrim.py` and prior discussions,
    # we are using this function's current behavior.
    return initial_solution_x

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]:
    """
    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]:
    """
    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 = TOVSolver(eos_filepath='path/to/eos', tidal=True, moment=True)

Parameters:

Name Type Description Default
eos_filepath str

The path to the EOS data table to be used.

 None
tidal bool

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

 False
moment bool

Whether to calculate moment of inertia or not. Default is False.

 required

Returns:

Type Description

None.
Source code in src/slmemulator/TOV_class.py 
 18
 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
def __init__(
    self,
    eos_filepath=None,
    tidal=False,
    solver="RK4",
    solve_ivp_kwargs=None,
    sol_pts=4000,
):
    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 = TOVSolver(eos_filepath='path/to/eos', tidal=True, moment=True)

    Parameters:
        eos_filepath (str): The path to the EOS data table to be used.

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

        moment (bool): Whether to calculate moment of inertia or not.
            Default is False.

    Returns:
        None.
    """

    # assign class variables
    self.tidal = tidal
    self.solver = solver
    self.sol_pts = sol_pts
    self.solve_ivp_kwargs = solve_ivp_kwargs  # this is only used in solve_ivp
    self.tol = 1e-9  # this is 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

    # load in the data (expects file with data headers for now)
    if eos_filepath is not None:

        file_name, self.eos_file_extension = os.path.splitext(eos_filepath)
        self.eos_name = (file_name.split("/"))[-1]

        # for data from compOSE
        if self.eos_file_extension == ".table":
            # extract from .table file
            self.eos_file = eos_filepath
            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 is True:
                dpdeps = np.gradient(self.pres_array, self.eps_array, edge_order=2)
                self.cs2_array = dpdeps
            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 == ".dat" or self.eos_file_extension == ".txt":

            # extract data and assign to arrays (assumes one sample/draw of EOS)
            self.eos_file = eos_filepath
            print("Loading EOS data from file: ", eos_filepath)
            eos_data = np.loadtxt(eos_filepath)  # , skiprows=2)
            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]
            if len(eos_data[1]) > 3:
                self.cs2_array = eos_data.T[3]
            else:
                self.cs2_array = None

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

        elif self.eos_file_extension == ".npz":

            self.eos_file = eos_filepath
            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"]
            self.cs2_array = eos_data["cs2"]

        # print congrats
        print("Woo it worked!")

    else:
        raise ValueError("No file specified.")

    return None

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/TOV_class.py 
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
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

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

A simple RK2 solver using the 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

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/TOV_class.py 
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
def RK2(self, f, x0, t0, te, N):
    r"""
    A simple RK2 solver using the 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): 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 = 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

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/TOV_class.py 
219
220
221
222
223
224
225
226
227
228
229
230
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
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): 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)
        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 
261
262
263
264
265
266
267
268
269
270
271
272
273
274
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
317
318
319
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 is True:
        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 is True:
            f = self.f_x(x, mass, pres, eps)
            q = self.q_x(x, mass, pres, eps, cs2)
            dydx = -(1.0 / x) * (y**2.0 + f * y + q)

    else:
        dpdx = 0.0
        dmdx = 0.0

        if self.tidal is True:
            dydx = 0.0

    if self.tidal is True:
        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 
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
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 
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
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 
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
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

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_max_mass (array): The array of total maximum mass values.
Source code in src/slmemulator/TOV_class.py 
429
430
431
432
433
434
435
436
437
438
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
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
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
556
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
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_max_mass (array): The array of total
            maximum mass values.
    """
    # list for storing results
    self.sols_varying_p0 = []
    # initial pressure
    pres_init = min(2.0, np.max(self.pres_array))  # Use np.max for clarity
    mass_init = 0.0

    if self.tidal is True:
        y_init = 2.0

    # set initial radius array
    low_pres = max(1e-3, np.min(self.pres_array))  # Use np.min for clarity
    x = np.geomspace(1e-8, 2.5, 50)
    pres_space = np.geomspace(low_pres, pres_init, len(x))
    self.pres_space = pres_space

    if self.tidal is True:
        self.yR = np.zeros(len(x))

    # set up arrays for the final results
    self.total_mass = np.zeros(len(x))
    self.total_radius = np.zeros(len(x))
    self.total_pres_central = np.zeros(len(x))

    # set up arrays (impermanent so this will work)
    max_mass = np.zeros(len(x))
    pres_central = np.zeros(len(x))
    max_radius = np.zeros(len(x))

    # interpolate the energy density
    self.eps_interp = interp1d(
        self.pres_array,
        self.eps_array,
        axis=0,
        kind="linear",
        fill_value="extrapolate",
    )

    if self.tidal is True:
        self.cs2_interp = interp1d(
            self.pres_array,
            self.cs2_array,
            axis=0,
            kind="linear",
            fill_value="extrapolate",
        )

    # loop over the TOV equations
    for i in range(len(x)):
        init_guess = []
        # initial conditions
        mass_arg = mass_init
        pres_arg = pres_space[i]
        init_guess.append(pres_arg)
        init_guess.append(mass_arg)

        if self.tidal is True:
            y_arg = y_init
            init_guess.append(y_arg)

        # Improve init_guess as one entry for all solutions
        # 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,
                }
            if self.tidal is True:
                sol = solve_ivp(
                    self.tov_equations_scaled,
                    [1e-3, 2.5],
                    [pres_arg, mass_arg, y_arg],
                    **self.solve_ivp_kwargs,
                )
            else:
                sol = solve_ivp(
                    self.tov_equations_scaled,
                    [1e-8, 2.5],
                    [pres_arg, mass_arg],
                    **self.solve_ivp_kwargs,
                )
            if not sol.success:
                print("Solver failed.")
                print(sol.message)
            xval = sol.t
            sol = sol.y
        else:
            raise ValueError(
                f'Solver, {self.solver} unknown. Must be "RK4", "RK2", "euler", or "solve_ivp".'
            )

        # maximum mass
        # Find the index where pressure becomes very small (approaching zero)
        # Use np.flatnonzero for more robust indexing
        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:
            # Handle cases where pressure never goes above threshold or immediately drops
            max_mass[i] = 0.0  # Or some other appropriate default
            # print(f"Warning: Pressure did not stay positive for central pressure {pres_arg}")

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

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

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

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

    # scale results (send back totals at the end using these)
    max_mass = max_mass * self.mass0
    max_radius = max_radius * self.rad0
    pres_central = pres_central * self.pres0
    self.total_mass = max_mass
    self.total_radius = max_radius
    self.total_pres_central = pres_central
    # max mass calculation, radius, and central pressure
    self.maximum_mass = np.max(max_mass)
    corr_radius_index = np.where(max_mass == self.maximum_mass)[0][0]
    self.corr_radius = max_radius[corr_radius_index]
    self.corr_pres = pres_central[corr_radius_index]

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

    print(
        "Max mass: ",
        self.maximum_mass,
        "Radius: ",
        self.corr_radius,
        "Central pressure: ",
        self.corr_pres,
    )

    # tidal deformability
    if self.tidal is True:
        self.tidal_deformability, self.k2 = self.tidal_def(
            self.yR, max_mass, max_radius
        )
        # print(
        # "Tidal deformability at all points: {}".format(self.tidal_deformability)
        # )
        # print("k2 at all points: {}".format(self.k2))

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

        # plot stuff
        plt.plot(max_radius, max_mass, label=r"TOV")
        # plt.plot(dense_rad, mrad_arr, label=r'Interpolant')
        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 is True:
            plt.plot(max_mass / max_radius, self.yR)
            plt.xlabel(r"$\beta$")
            plt.ylabel(r"y(r)")
            plt.savefig("yr_scaled.png")
            plt.show()
            # fig = plt.figure(figsize=(8,6), dpi=200)
            plt.plot(max_radius, self.tidal_deformability, 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()
            # fig = plt.figure(figsize=(8,6), dpi=200)
            plt.plot(max_mass, self.tidal_deformability, 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()
            # fig = plt.figure(figsize=(8,6), dpi=200)
            plt.plot(max_mass / max_radius, self.k2, 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()
            # fig = plt.figure(figsize=(8,6), dpi=200)
            plt.plot(max_radius, self.k2, 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 is True:
        tov_data = np.column_stack(
            [self.total_radius, self.total_mass, self.total_pres_central]
        )
        file_name = "TOV_data/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 is True:
        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 
721
722
723
724
725
726
727
728
729
730
731
732
733
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 array

The array of central densities for each EOS used.
Source code in src/slmemulator/TOV_class.py 
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
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 (array): The array of central densities for
            each EOS used.
    """

    # interpolate the EOS to find the proper central densities
    p_n_interp = interp1d(
        self.pres_array_unscaled,
        self.nB_array,
        kind="cubic",
        fill_value="extrapolate",
    )

    # solve at the proper central pressure for nB_central
    c_dens = p_n_interp(self.corr_pres)

    return 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 array

The array of values of the radius for each EOS used.
Source code in src/slmemulator/TOV_class.py 
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
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 (array): The array of values of the radius
            for each EOS used.
    """

    # interpolate the mass radius result?
    m_r_interp = interp1d(
        self.total_mass,
        self.total_radius,
        kind="linear",
        fill_value="extrapolate",
    )
    return m_r_interp(1.4)

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)