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
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
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.

Source code in src/slmemulator/SLM.py
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.
    """
    # 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: str = None)

Cleans up specified files and directories within a given directory.

directory (str): The path to the directory to clean. Defaults to the current working directory if None.

Source code in src/slmemulator/cleanData.py
def clean_directory(directory: str = None):
    """
    Cleans up specified files and directories within a given directory.

    Parameters:
    directory (str): The path to the directory to clean. Defaults to the
                     current working directory if None.
    """
    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_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: Path = None, use_mseos: bool = True, is_parametric_run: bool = True) -> dict

Returns a dictionary of resolved paths for data input/output and other project directories. Prioritizes explicit arguments for output_base_dir, then sensible defaults. Conditionally sets 'current' EOS/TOV/results paths based on the use_mseos flag. Adjusts base for general data directories (EOS_Data, TOV_data, testData, trainData) to output_base_dir if provided (for local working directories), otherwise defaults them to be under src_dir (for package internal use).

Parameters:

Name Type Description Default
output_base_dir Path

The base directory for all generated outputs. Defaults to the project root.

None
use_mseos bool

If True, 'current' paths will point to MSEOS-related directories. If False, 'current' paths will point to QEOS-related directories. Defaults to True.

True
is_parametric_run bool

If True, indicates a parametric run. This parameter is available for future path logic, but currently the primary factor for data directory location is output_base_dir. Defaults to True.

True

Returns:

Name Type Description
dict dict

A dictionary containing all relevant path configurations.

Source code in src/slmemulator/config.py
def get_paths(
    output_base_dir: Path = None,
    use_mseos: bool = True,
    is_parametric_run: bool = True,
) -> dict:
    """
    Returns a dictionary of resolved paths for data input/output and other project directories.
    Prioritizes explicit arguments for output_base_dir, then sensible defaults.
    Conditionally sets 'current' EOS/TOV/results paths based on the use_mseos flag.
    Adjusts base for general data directories (EOS_Data, TOV_data, testData, trainData)
    to output_base_dir if provided (for local working directories),
    otherwise defaults them to be under src_dir (for package internal use).

    Args:
        output_base_dir (Path, optional): The base directory for all generated outputs.
                                         Defaults to the project root.
        use_mseos (bool): If True, 'current' paths will point to MSEOS-related directories.
                          If False, 'current' paths will point to QEOS-related directories.
                          Defaults to True.
        is_parametric_run (bool): If True, indicates a parametric run.
                                  This parameter is available for future path logic,
                                  but currently the primary factor for
                                  data directory location is `output_base_dir`.
                                  Defaults to True.

    Returns:
        dict: A dictionary containing all relevant path configurations.
    """
    # Use explicit output_base_dir or default to the calculated PROJECT_ROOT
    output_base = output_base_dir or PROJECT_ROOT

    # src_dir will always be under the primary base (output_base or PROJECT_ROOT)
    src_dir = output_base / DEFAULT_SRC_SUBDIR_NAME

    # Determine the base for general data directories (EOS_Data, TOV_data, testData, trainData)
    # If output_base_dir is explicitly provided, these data directories should be
    # relative to it. Otherwise, they remain within the src_dir of the package.
    if output_base_dir is not None:
        data_root_dir = output_base
    else:
        data_root_dir = src_dir

    paths = {
        "project_root": output_base,
        "src_dir": src_dir,
        # EOS_Codes_dir should point to the location within the *installed package*
        # This ensures that whether installed or in editable mode, it finds the right place.
        "eos_codes_dir": pkg_resources.files("slmemulator").joinpath(
            DEFAULT_EOS_CODES_SUBDIR_NAME
        ),
        "plots_dir": data_root_dir / DEFAULT_PLOTS_SUBDIR_NAME,
        "results_dir": data_root_dir / DEFAULT_RESULTS_SUBDIR_NAME,
        "docs_dir": output_base / DEFAULT_DOCS_SUBDIR_NAME,
        "tests_dir": output_base / DEFAULT_TESTS_SUBDIR_NAME,
        "tutorials_dir": output_base / DEFAULT_TUTORIALS_SUBDIR_NAME,
        # Data directories whose location depends on `output_base_dir`
        "eos_data_dir": data_root_dir
        / DEFAULT_EOS_DATA_SUBDIR_NAME,  # For user-provided/existing EOS data
        "tov_data_dir": data_root_dir / DEFAULT_TOV_DATA_SUBDIR_NAME,
        "test_data_dir": data_root_dir / DEFAULT_TEST_DATA_SUBDIR_NAME,
        "train_path": data_root_dir / DEFAULT_TRAIN_DATA_SUBDIR_NAME,
        # EOS_files_dir contains generated EOS files (MSEOS/QEOS subdirs)
        # This is where output of MSEOS.py/Quarkyonia.py goes, typically under src/
        "eos_files_dir": data_root_dir / DEFAULT_EOS_FILES_SUBDIR_NAME,
        "eos_data_dir": pkg_resources.files("slmemulator").joinpath(
            DEFAULT_EOS_DATA_SUBDIR_NAME
        ),
        # Specific paths for QEOS and MSEOS (always defined relative to their respective bases)
        "qeos_path_specific": (data_root_dir / DEFAULT_EOS_FILES_SUBDIR_NAME) / "QEOS",
        "mseos_path_specific": (data_root_dir / DEFAULT_EOS_FILES_SUBDIR_NAME)
        / "MSEOS",
        # Correcting specific TOV data paths to use data_root_dir (for user data)
        "qeos_tov_path_specific": (data_root_dir / DEFAULT_TOV_DATA_SUBDIR_NAME)
        / "QEOS",
        "mseos_tov_path_specific": (data_root_dir / DEFAULT_TOV_DATA_SUBDIR_NAME)
        / "MSEOS",
        # General SLM specific result/plot subdirectories (non-parametric)
        "slm_res_mseos_specific": (
            data_root_dir
            / DEFAULT_RESULTS_SUBDIR_NAME
            / "MSEOS"
            / DEFAULT_SLM_SUBDIR_NAME
        ),
        "slm_res_qeos_specific": (
            data_root_dir
            / DEFAULT_RESULTS_SUBDIR_NAME
            / "QEOS"
            / DEFAULT_SLM_SUBDIR_NAME
        ),
        "slm_plots_mseos_specific": (
            data_root_dir
            / DEFAULT_PLOTS_SUBDIR_NAME
            / "MSEOS"
            / DEFAULT_SLM_SUBDIR_NAME
        ),
        "slm_plots_qeos_specific": (
            data_root_dir / DEFAULT_PLOTS_SUBDIR_NAME / "QEOS" / DEFAULT_SLM_SUBDIR_NAME
        ),
        # Parametric SLM specific result/plot subdirectories
        "slm_res_mseos_parametric_specific": (
            data_root_dir
            / DEFAULT_RESULTS_SUBDIR_NAME
            / "MSEOS"
            / DEFAULT_PSLM_SUBDIR_NAME
        ),
        "slm_res_qeos_parametric_specific": (
            data_root_dir
            / DEFAULT_RESULTS_SUBDIR_NAME
            / "QEOS"
            / DEFAULT_PSLM_SUBDIR_NAME
        ),
        "slm_plots_mseos_parametric_specific": (
            data_root_dir
            / DEFAULT_PLOTS_SUBDIR_NAME
            / "MSEOS"
            / DEFAULT_PSLM_SUBDIR_NAME
        ),
        "slm_plots_qeos_parametric_specific": (
            data_root_dir
            / DEFAULT_PLOTS_SUBDIR_NAME
            / "QEOS"
            / DEFAULT_PSLM_SUBDIR_NAME
        ),
    }

    # Conditionally set the "current" paths based on the use_mseos flag AND `is_parametric_run`
    if use_mseos:
        paths["current_eos_input_dir"] = paths["mseos_path_specific"]
        paths["current_tov_data_dir"] = paths["mseos_tov_path_specific"]
        if is_parametric_run:
            paths["current_slm_results_dir"] = paths[
                "slm_res_mseos_parametric_specific"
            ]
            paths["current_slm_plots_dir"] = paths[
                "slm_plots_mseos_parametric_specific"
            ]
        else:  # Non-parametric SLM run for MSEOS
            paths["current_slm_results_dir"] = paths["slm_res_mseos_specific"]
            paths["current_slm_plots_dir"] = paths["slm_plots_mseos_specific"]
    else:  # QEOS
        paths["current_eos_input_dir"] = paths["qeos_path_specific"]
        paths["current_tov_data_dir"] = paths["qeos_tov_path_specific"]
        if is_parametric_run:
            paths["current_slm_results_dir"] = paths["slm_res_qeos_parametric_specific"]
            paths["current_slm_plots_dir"] = paths["slm_plots_qeos_parametric_specific"]
        else:  # Non-parametric SLM run for QEOS
            paths["current_slm_results_dir"] = paths["slm_res_qeos_specific"]
            paths["current_slm_plots_dir"] = paths["slm_plots_qeos_specific"]

    # Ensure all paths are Path objects
    return {k: Path(v) for k, v in paths.items()}

pSLM

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

Source code in src/slmemulator/pSLM.py
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
@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
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
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
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
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
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

Information about the code: 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: 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)

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

Evaluate pressure from number density using interpolation

Source code in src/slmemulator/scaledTOV.py
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/scaledTOV.py
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/scaledTOV.py
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/scaledTOV.py
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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)