utils.py

import numpy as np
import pandas as pd
from numba import njit, prange, float64, int64
from shapely.geometry import LineString, Point
from sklearn.neighbors import KDTree
from sklearn.decomposition import PCA
from astropy.coordinates import SkyCoord
import h5py
from scipy.odr import *
from numpy import linalg as LA


@njit(parallel=True, nogil=True, fastmath=True)
def hough(vec, thetas):
    """
    This function transforms points in the image space to the Hough line equation parameter space.
    :param vec: input vector, 2d array
    :param thetas: range of angles in the Hough parameter space
    :return: Hough curves of (theta, rho) parameter pairs
    """
    # Cache some reusable values
    cos_t = np.cos(thetas)
    sin_t = np.sin(thetas)
    num_thetas = len(thetas)

    curves = np.empty((len(vec), num_thetas, 2), dtype=float64)
    for i in range(len(vec)):
        curves[i, :, 0] = thetas  # theta parameter
        curves[i, :, 1] = vec[i][0] * cos_t + vec[i][1] * sin_t  # rho parameter

    return curves


def norm_curves(curves):
    """
    Z-score normalize the hough curves.
    :param curves: hough function output. Hough space curves for each image space point.
    :return curves_conc: concatenated hough curves.
    :return curves_conc_normed: normalized concatenated curves
    :return curves_norm: normalized curves split according to image space points
    :return mean_all: mean of all curves
    :return std_all: standard deviation of all curves
    """
    # get shape parameters
    num_classes, num_points_per_class, num_dimensions = curves.shape
    # concatenate the curves
    curves = np.ascontiguousarray(curves)
    curves_conc = curves.reshape(-1, num_dimensions)
    # calculate mean and standard deviation
    mean_all = np.mean(curves_conc, axis=0)
    std_all = np.std(curves_conc, axis=0)
    # z-score normalize all curves
    curves_conc_norm = (curves_conc - mean_all) / std_all
    # split curves back according to image space points
    curves_norm = curves_conc_norm.reshape(
        num_classes, num_points_per_class, num_dimensions
    )

    return curves_conc, curves_conc_norm, curves_norm, mean_all, std_all


def find_intersections(curves):
    """
    Finds the intersections of any two curves in the Hough space
    :param curves: hough function output. Hough space curves for each image space point.
    :return intersection points: intersection points of the hough space curves
    """
    # Convert the curves to Shapely LineString objects
    lines = [LineString(curve) for curve in curves]
    # Find the intersection points for all pairs of curves
    intersection_points = []
    for i in range(len(lines)):
        for j in range(i + 1, len(lines)):
            inter = lines[i].intersection(lines[j])
            if isinstance(inter, Point):
                intersection_points.append((inter.x, inter.y))
            elif inter.geom_type == "MultiPoint":
                for point in inter:
                    intersection_points.append((point.x, point.y))
    return np.array(intersection_points)


def func_PCA(x, y):
    """
    This function performs Principal Component Analysis.
    :param x: x-coordinate, 1d-array
    :param y: y-coordinate, 1d-array
    :return: axis ratio, short over long axis
    """
    # Create a 2D array of x and y coordinates
    xy = np.column_stack((x, y))
    # Perform principal component analysis
    pca = PCA()
    pca.fit(xy)
    # Get the eigenvalues and eigenvectors
    eigenvalues = pca.explained_variance_
    eigenvectors = pca.components_
    # Get the axis lengths and axis ratio
    major_axis = 2 * np.sqrt(eigenvalues[0])
    minor_axis = 2 * np.sqrt(eigenvalues[1])
    axis_ratio = minor_axis / major_axis
    return axis_ratio, minor_axis / 2, major_axis / 2, eigenvectors


def find_max(
    tree,
    class_info,
    curves_conc_norm,
    num_points_per_class,
    radius,
    rho_min=None,
    rho_max=None,
):
    """
    Find the most crowded point in the Hough space
    :param tree: KDtree
    :param class_info: which point in the Hough space belongs to which point in the image space
    :param curves_conc_norm: concatenated and z-score normalized Hough curves
    :param num_points_per_class: number of Hough points in a curve = number of thetas
    :param radius: search radius around a given point in the Hough space
    :param rho_min: minimum rho to consider
    :param rho_max: maximum rho to consider
    :return: total number of points around the most crowded point, number of different image space points responsible,
    indices of the neighboring points corresponding to individual image space points, Hough parameter pair of the most
    crowded point.
    """
    neighbors_indices = tree.query_radius(
        curves_conc_norm, r=radius
    )  # query the KDTree for points in a given radius
    max_total_count = 0
    max_class_count = 0
    max_neighbor_indices = np.array([])
    most_crowded_point = None
    # loop over the neighbors of all points
    for i, neighbors in enumerate(neighbors_indices):
        class_idx, point_idx = divmod(
            i, num_points_per_class
        )  # image space point and index within the Hough curve
        rho = curves_conc_norm[i, 1]
        if (rho_min is None or (rho_min <= rho)) and (
            rho_max is None or (rho <= rho_max)
        ):
            neighbor_classes = np.unique(
                class_info[neighbors]
            )  # indices of different image space points
            class_count = len(
                neighbor_classes
            )  # number of different image space points within the search radius
            total_count = len(
                neighbors
            )  # total number of Hough points within the search radius
            # Update the total number of points and total number of different classes if larger than for the previous point
            if total_count > max_total_count or (
                total_count == max_total_count and class_count > max_class_count
            ):
                max_total_count = total_count
                max_class_count = class_count
                max_neighbor_indices = neighbor_classes.copy()
                most_crowded_point = (class_idx, point_idx)
    return max_total_count, max_class_count, max_neighbor_indices, most_crowded_point


def tree_skl(vec, curves_shape, curves_conc_norm, radius, rho_min=None, rho_max=None):
    """
    Construct a KDTree and find nearest neighbors of all points in the Hough space
    :param vec: satellite coordinates in the restframe of the host, 2d-array
    :param curves_shape: shape of the hough curves, (number of satellites, len(thetas), number of dimensions), tuple
    :param curves_conc_norm: concatenated z-score normalized Hough curves
    :param radius: search radius around a given point
    :param rho_min: minimum rho to consider
    :param rho_max: maximum rho to consider
    :return: total number of points around the most crowded point, number of different image space points responsible,
    indices of the neighboring points corresponding to individual image space points, Hough parameter pair of the most
    crowded point.
    """
    # Get the shape of the normalized data
    num_classes, num_points_per_class, num_dimensions = curves_shape
    # Create a KD-Tree
    tree = KDTree(curves_conc_norm)
    # Class information corresponding to the flattened data
    class_info = np.repeat(np.arange(num_classes), num_points_per_class)
    # Find the most crowded point
    max_total_count, max_class_count, max_neighbor_indices, most_crowded_point = (
        find_max(
            tree,
            class_info,
            curves_conc_norm,
            num_points_per_class,
            radius,
            rho_min,
            rho_max,
        )
    )
    # Calculate the axis ratio of the resulting structure via PCA
    ar, _, _, _ = func_PCA(
        vec[:, 0][max_neighbor_indices], vec[:, 1][max_neighbor_indices]
    )
    return (
        max_total_count,
        max_class_count,
        max_neighbor_indices,
        most_crowded_point,
        ar,
    )


def tree_skl_loop(
    vec, curves_shape, curves_conc_norm, radii, rho_min=None, rho_max=None
):
    """
    Construct a KDTree and query around every Hough point for a range of different search radii.
    :param vec: Image space coordinates, 2d-array
    :param curves_shape: shape of the Hough curves
    :param curves_conc_norm: concatenated z-score Hough curves
    :param radii: range of search radii to query
    :param rho_min: minimum rho to consider
    :param rho_max: maximum rho to consider
    :return: Results from the search radius which maximizes both flatness and member population
    """
    # Get the shape of the normalized data
    num_classes, num_points_per_class, num_dimensions = curves_shape
    # Create a KD-Tree
    tree = KDTree(curves_conc_norm)
    # Class information corresponding to the flattened data
    class_info = np.repeat(np.arange(num_classes), num_points_per_class)

    total_count_master = np.empty(1)
    class_count_master = np.empty(1)
    neighbor_indices_master = np.array([])
    most_crowded_point_master = ()
    ratio_master = 1
    radius_master = np.empty(1)
    ar_master = np.empty(1)
    sem_minor_master = np.empty(1)
    sem_major_master = np.empty(1)

    # Query the tree to get points within a range of radii r for each point
    for idx, radius in enumerate(radii):
        max_total_count, max_class_count, max_neighbor_indices, most_crowded_point = (
            find_max(
                tree,
                class_info,
                curves_conc_norm,
                num_points_per_class,
                radius,
                rho_min,
                rho_max,
            )
        )
        # Calculate the axis ratio of the resulting structure via PCA
        ar, sem_minor, sem_major, _ = func_PCA(
            vec[:, 0][max_neighbor_indices], vec[:, 1][max_neighbor_indices]
        )
        # flatness over number of members squared
        ratio = ar / (len(max_neighbor_indices) ** 2)
        # update the values of the new ratio is smaller than the previous one
        if ratio < ratio_master:
            total_count_master = max_total_count
            class_count_master = max_class_count
            neighbor_indices_master = max_neighbor_indices
            most_crowded_point_master = most_crowded_point
            ratio_master = ratio
            radius_master = radius
            ar_master = ar
            sem_minor_master = sem_minor
            sem_major_master = sem_major

    return (
        total_count_master,
        class_count_master,
        neighbor_indices_master,
        most_crowded_point_master,
        radius_master,
        ar_master,
        sem_minor_master,
        sem_major_master,
    )


def get_fit_params(curves, max_point):
    """
    Transform Hough parameters to image space best fit parameters.
    :param curves: Original un-normalized Hough curves
    :param max_point: Index pair of the most crowded point, tuple
    :return: best fit slope and intercept
    """
    # Get best fit parameters
    theta_best = curves[max_point][0]
    rho_best = curves[max_point][1]
    # Back-transform to image space
    slope_best = -np.cos(theta_best) / np.sin(theta_best)
    inter_best = rho_best / np.sin(theta_best)
    return slope_best, inter_best


def res_var(vec, slope, inter):
    """
    Calculate the perpendicular residual variance between the Hough members and the best fit line
    :param vec: image space coordinates, 2d-array
    :param slope: best fit slope, float
    :param inter: best fit intercept, float
    :return: residual variance
    """
    x_data, y_data = vec[:, 0], vec[:, 1]
    orthogonal_residuals = np.abs((slope * x_data - y_data + inter)) / np.sqrt(
        slope**2 + 1
    )
    variance = np.sum(orthogonal_residuals**2) / (len(vec) - 2)
    return variance


def arcsec_to_kpc(arcsec, d):
    """
    Transform arcseconds to kpc
    :param arcsec: Angular distance
    :param d: distance to object
    :return: physical distance in kpc
    """
    radians = np.radians(arcsec / 3600)
    return radians * d * 1000


@njit(float64[:, :](float64[:], float64[:], int64))
def monte_carlo(central_gal, dist, n):
    """
    Return root-mean-square distance of n artificial satellites to the best fit line
    :param central_gal: coordinates of the host galaxy, 2d-array
    :param dist: 2d separation between satellite and host
    :param n: number of satellites
    :return: artificial satellite system, radial distribution preserved
    """
    x_art = np.empty(n)  # create array for new x values
    y_art = np.empty(n)  # create array for new y values
    for i in prange(n):
        theta = (
            2 * np.pi * np.random.random()
        )  # generate random angles in the interval 0 to 2pi
        x_art[i] = dist[i] * np.sin(theta) + central_gal[0]
        y_art[i] = dist[i] * np.cos(theta) + central_gal[1]
    sat_art = np.array(
        list(zip(x_art, y_art))
    )  # combine new values to array of (x,y) pairs

    return sat_art


def host_restframe(ras, decs, ra_host, dec_host, d_host):
    """
    Transform satellite world coordinates to the host restframe.
    :param ras: right ascension in degrees, 1d-array
    :param decs: declination in degrees, 1d-array
    :param ra_host: host right ascension in degrees, float
    :param dec_host: host declination in degrees, float
    :param d_host: distance to host in Mpc, float
    :return: satellite coordinates in kpc relative to the host, 2d-array
    """
    # projects the satellite coordinates to the same RA as the host
    coords_with_host_ra = SkyCoord(
        np.ones(len(ras)) * ra_host, decs, unit="deg", frame="icrs"
    )
    # projects the satellite coordinates to the same Dec as the host
    coords_with_host_dec = SkyCoord(
        ras, np.ones(len(ras)) * dec_host, unit="deg", frame="icrs"
    )
    coord_host = SkyCoord(ra_host, dec_host, unit="deg", frame="icrs")

    x_rel = (
        d_host
        * np.tan(np.deg2rad(coord_host.separation(coords_with_host_dec).deg))
        * 1000
    )
    y_rel = (
        d_host
        * np.tan(np.deg2rad(coord_host.separation(coords_with_host_ra).deg))
        * 1000
    )

    # give sign according to position. By convention those with lower RA/DEC than host have a minus sign
    for i in np.arange(0, len(x_rel)):
        if ras[i] < ra_host:
            x_rel[i] = -x_rel[i]
    for i in np.arange(0, len(y_rel)):
        if decs[i] < dec_host:
            y_rel[i] = -y_rel[i]

    return np.array(list(zip(x_rel, y_rel)))


def import_data(
    host,
    companion,
    table_path,
    spec="",
    restframe_ancor="host",
    dSph_only=False,
    restframe_shift=False,
):
    """
    Import data from table.
    :param host: host name
    :param companion: companion name, if any
    :param table_path: path to tables
    :param spec: specification on the dataframe
    :param restframe_ancor: defines what should be used as the restframe ancor (host or mean of host and companion)
    :param dSph_only: only consider dwarf spheroidals
    :param restframe_shift: perform shift to host rest frame if not already done
    :return: satellite coordinates in the rest frame of the host in kpc
    """
    df = pd.read_csv(table_path + f"{host}_system{spec}.csv")
    if dSph_only:
        df = df[df["type"] == "dSph"].reset_index(drop=True)
    if companion is not None:
        dwarf_df = df.loc[~df["ID"].isin([host, companion])].reset_index(drop=True)
        host_df = df.loc[df["ID"].eq(host)].reset_index(drop=True)
        companion_df = df.loc[df["ID"].eq(companion)].reset_index(drop=True)
        companion_host_df = df.loc[df["ID"].isin([host, companion])].reset_index(
            drop=True
        )
    else:
        dwarf_df = df.loc[df["ID"] != host].reset_index(drop=True)
        host_df = df.loc[df["ID"].eq(host)].reset_index(drop=True)
    if restframe_shift:
        if restframe_ancor == "host":
            ras, decs, ra_host, dec_host, d_host = (
                np.array(df.ra),
                np.array(df.dec),
                host_df.ra.values,
                host_df.dec.values,
                host_df.D.values,
            )
        elif restframe_ancor == "centroid_host_companion":
            ras, decs, d_host = np.array(df.ra), np.array(df.dec), host_df.D.values
            ra_host, dec_host = np.array(
                [np.mean(companion_host_df.ra), np.mean(companion_host_df.dec)]
            )
            print(f"Companion host centroid is at: {ra_host}, {dec_host}.")
        else:
            print("Not implemented. Choose between host and centroid_host_companion")
        ras_sat, decs_sat = np.array(dwarf_df.ra), np.array(dwarf_df.dec)
        if host == "M31":
            d_host = d_host / 1000.0
        vec_all = host_restframe(ras, decs, ra_host, dec_host, d_host)
        x_rel_all, y_rel_all = vec_all[:, 0], vec_all[:, 1]
        df["x_rel"], df["y_rel"] = x_rel_all, y_rel_all
        df.to_csv(table_path + f"{host}_system{spec}.csv", index=False)
        vec_sat = host_restframe(ras_sat, decs_sat, ra_host, dec_host, d_host)
        x_rel, y_rel = vec_sat[:, 0], vec_sat[:, 1]
        dwarf_df["x_rel"], dwarf_df["y_rel"] = x_rel, y_rel
    else:
        x_rel, y_rel = dwarf_df["x_rel"], dwarf_df["y_rel"]
    vec = np.array(list(zip(x_rel, y_rel)))  # 2d-array
    return vec, dwarf_df


def import_sim_data(host, table_path, viewing_angles):
    """
    Imports simulation data from hdf5 file.
    :param host: host name
    :param table_path: path to folder
    :param viewing_angles: number of random viewing angles in mock observations
    :return: dictionaries containing simulation data, one dictionary for each extension
    """
    with h5py.File(table_path + "{}_analogs.hdf5".format(host), "r") as f:
        # Create empty dictionaries to store data for each group
        tng50_data = {}
        tng50_ext_data = {}

        # Loop through datasets in the first group
        for dataset_name, dataset in f["TNG50"].items():
            data = dataset[()]
            tng50_data[dataset_name] = data

        # Loop through datasets in the second group
        for dataset_name, dataset in f["TNG50 Extended"].items():
            data = dataset[()]
            tng50_ext_data[dataset_name] = data
    if viewing_angles == 10:
        return tng50_data
    if viewing_angles == 100:
        return tng50_ext_data


def kin_corr(
    host,
    companion,
    vec,
    voter_indices,
    vel,
    slope,
    inter,
    voters_data,
    n_vel_data,
    table_path,
    anchor,
    run_on,
    data_incomplete=False,
    pick_vels=None,
    sim_mass_arr=None,
):
    """
    Find the number of corotating satellites.
    :param host: system name
    :param companion: companion name, if any
    :param vec: satellite vector in the host restframe, 2d-array
    :param voter_indices: list of indices pointing to the structure members
    :param vel: satellite velocity
    :param slope: best fit slope
    :param inter: best fit intercept
    :param voters_data: number of on-plane satllites in the data
    :param n_vel_data: number of satellites with velocities in the data
    :param table_path: table path
    :param anchor: reference point for the phase-space of the system (satellite system/host/weighted group mean)
    :param run_on: run on data or simulations; data/sim
    :param data_incomplete: define if there are velocities missing in the data; True/False
    :param pick_vels: pick number of analog velocities based on absolute number or percentual number of on-plane
    velocities in the data. This is an option because the Hough transform does not always pick the exact same number
    of on-plane dwarfs as in the data.
    :param sim_mass_arr: satellite mass from simulation, to pick the n highest mass satellites
    positions and velocities
    :return: number of corotating satellites, modified data frame if passed as input
    """
    vec_voters = vec[voter_indices]
    # set velocities of non-voters to nan
    vel[np.setdiff1d(np.arange(len(vel)), voter_indices)] = np.nan
    number_of_vels = len(voter_indices)
    if data_incomplete:
        if pick_vels == "absolute":
            number_of_vels = n_vel_data
        if pick_vels == "percent":
            number_of_vels = int(
                np.round(len(voter_indices) * n_vel_data / voters_data)
            )  # round to get an absolute number
        # get dark masses of hough voters
        sim_mass_arr_voters = sim_mass_arr[voter_indices]
        # for the voter subset pick the brightest x satellites which have velocities available in the data
        idx_brightest = idx_n_largest(sim_mass_arr_voters, number_of_vels)
        # print(f'N_on_plane {len(voter_indices)}, N_w_vels {number_of_vels}')
        # set all velocities which are not the brightest x hough voters to nan
        vel[np.setdiff1d(np.arange(len(vel)), voter_indices[idx_brightest])] = np.nan
    if run_on == "data":
        df_w_host = pd.read_csv(table_path + "{}_system.csv".format(host))
        vel_all = np.array(df_w_host["v"])
        host_vel = np.array(df_w_host["v"].loc[df_w_host["ID"] == host])
    if run_on == "sim":
        host_vel = 0.0
    if anchor == "sat_system":
        r0 = np.array([np.mean(vec_voters[:, 0]), np.mean(vec_voters[:, 1])])
        v_anchor = np.nanmean(vel)
    elif anchor == "host":
        # r0 = np.array([0., 0.])
        r0 = np.array([np.mean(vec_voters[:, 0]), np.mean(vec_voters[:, 1])])
        v_anchor = host_vel
    else:
        print(
            "Choose satellite system or host as a reference, alternative options not yet implemented."
        )

    # calculate satellite distances from the minor axis
    sl_ort = 1 / -slope
    inter_ort = (slope + 1 / slope) * r0[0] + inter
    d_maj = (vec[:, 1] - (sl_ort * vec[:, 0] + inter_ort)) / np.sqrt(sl_ort**2 + 1)
    v_rel = vel - v_anchor
    coord_rel_host = np.array([0.0, 0.0])
    d_maj_host = dist_along_maj(r0, coord_rel_host, slope, inter)
    if run_on == "data":
        # add v_rel and d_maj to original df for pv plot
        df_w_host.loc[:, "v_rel"] = vel_all - v_anchor
        is_host = df_w_host["ID"].isin([host])  # host mask
        is_companion = df_w_host["ID"].isin([companion])  # companion mask
        is_sat = ~df_w_host["ID"].isin([host, companion])  # satellite mask
        coord_rel_companion = np.array(
            [df_w_host["x_rel"].loc[is_companion], df_w_host["y_rel"].loc[is_companion]]
        )
        # calculate distance from minor axis along the major axis for host and companion
        d_maj_companion = dist_along_maj(r0, coord_rel_companion, slope, inter)
        # for host
        df_w_host.loc[is_host, "d_maj"] = d_maj_host
        # for companion
        df_w_host.loc[is_companion, "d_maj"] = d_maj_companion
        # for satellites
        df_w_host.loc[is_sat, "d_maj"] = d_maj
        df_w_host.to_csv(table_path + "{}_system_pv.csv".format(host), index=False)

    if anchor == "sat_system":
        anchor_point = 0.0
    if anchor == "host":
        anchor_point = d_maj_host

    neg_neg_count = np.sum((d_maj < anchor_point) & (v_rel < 0))
    neg_pos_count = np.sum((d_maj < anchor_point) & (v_rel > 0))
    pos_pos_count = np.sum((d_maj > anchor_point) & (v_rel > 0))
    pos_neg_count = np.sum((d_maj > anchor_point) & (v_rel < 0))

    diag_pos_count = neg_neg_count + pos_pos_count
    diag_neg_count = neg_pos_count + pos_neg_count

    n_corot = max(diag_pos_count, diag_neg_count)

    if run_on == "sim":
        return n_corot, number_of_vels
    else:
        return n_corot


def idx_n_largest(arr, n):
    """
    Get the indices of the n largest entries in an array.
    :param arr: imput array, 1d
    :param n: number of largest entries
    :return: indices of n largest entries, 1d-array
    """
    indices = np.argpartition(arr, -n)[-n:]
    return np.sort(indices)


def dist_along_maj(r0, vec, slope, inter):
    """
    Calculate perpendicular distance from the minor axis along the major axis of the structure
    :param r0: reference point
    :param vec: coordinates, 2d-array
    :param slope: best-fit slope
    :param inter: best-fit intercept
    :return: distance along the major axis
    """
    vec = vec.reshape(-1, 2)
    sl_ort = 1 / -slope
    inter_ort = (slope + 1 / slope) * r0[0] + inter
    d_maj = (vec[:, 1] - (sl_ort * vec[:, 0] + inter_ort)) / np.sqrt(sl_ort**2 + 1)
    return d_maj


def linear_func(p, x):
    """Defines a linear function"""
    m, c = p
    return m * x + c


def tls(x, y):
    """
    Perform total least square fitting with input data (x,y)
    :param x: x coordinates
    :param y: y coordinates
    :return: slope, intercept and residual variance
    """
    linear_model = Model(linear_func)
    data = RealData(x, y)
    odr_job = ODR(data, linear_model, beta0=[0.0, 1.0])
    out = odr_job.run()
    return out.beta[0], out.beta[1], out.res_var


def toi(vec, r0):
    I = np.array([[1, 0], [0, 1]])
    T0 = np.array([[0.0, 0.0], [0.0, 0.0]])
    for i in range(len(vec)):
        T0 += LA.norm(np.subtract(vec[i], r0)) ** 2 * I - np.dot(
            np.transpose(np.subtract(vec[i], r0)), np.subtract(vec[i], r0)
        )
    eig_val, evec = LA.eig(T0)
    return eig_val, evec


def round_to_nearest_10_or_100(number):
    """
    This function rounds to the nearest multiple of 10 if the number is below 100 and the nearest multiple of 100 if
    the number is above 100.
    :param number: Input number, float or int
    :return: Nearest multiple of 10 or 100, int
    """
    if number <= 100:
        return int((number + 5) // 10 * 10)  # Round to nearest 10
    else:
        return int((number + 50) // 100 * 100)  # Round to nearest 100


def generate_markers(num_markers):
    markers = [
        ".",
        ",",
        "o",
        "v",
        "^",
        "<",
        ">",
        "1",
        "2",
        "3",
        "4",
        "8",
        "s",
        "p",
        "P",
        "*",
        "h",
        "H",
        "+",
        "x",
        "X",
        "D",
        "d",
        "|",
        "_",
        0,
        1,
        2,
        3,
        4,
        5,
        6,
        7,
        8,
        9,
        10,
        11,
    ]
    return markers[:num_markers]