78def calculate_z(adata, n_features=5000, batches=10, 
 79                batch_size=100) -> ad.AnnData:
 80    """
 81    Function to calculate Z matrices for all groups in both training 
 82    and testing data.
 83
 84    Parameters
 85    ----------
 86    adata : ad.AnnData
 87        created by `scmkl.create_adata()` with `adata.uns.keys()`: 
 88        `'train_indices'`, and `'test_indices'`. 
 89
 90    n_features : int
 91        Number of random feature to use when calculating Z; used for 
 92        scalability.
 93
 94    batches : int
 95        The number of batches to use for the distance calculation.
 96        This will average the result of `batches` distance calculations
 97        of `batch_size` randomly sampled cells. More batches will converge
 98        to population distance values at the cost of scalability.
 99
100    batch_size : int
101        The number of cells to include per batch for distance
102        calculations. Higher batch size will converge to population
103        distance values at the cost of scalability.
104        If `batches*batch_size > num_training_cells`,
105        `batch_size` will be reduced to 
106        `int(num_training_cells / batches)`.
107
108    Returns
109    -------
110    adata : ad.AnnData
111        `adata` with Z matrices accessible with `adata.uns['Z_train']` 
112        and `adata.uns['Z_test']`.
113
114    Examples
115    --------
116    >>> adata = scmkl.estimate_sigma(adata)
117    >>> adata = scmkl.calculate_z(adata)
118    >>> adata.uns.keys()
119    dict_keys(['Z_train', 'Z_test', 'sigmas', 'train_indices', 
120    'test_indices'])
121    """
122    # Number of groupings taking from group_dict
123    n_pathway = len(adata.uns['group_dict'].keys())
124    D = adata.uns['D']
125
126    sq_i_d = np.sqrt(1/D)
127
128    # Capturing training and testing sizes
129    train_len = len(adata.uns['train_indices'])
130    test_len = len(adata.uns['test_indices'])
131
132    if batch_size * batches > len(adata.uns['train_indices']):
133        old_batch_size = batch_size
134        batch_size = int(len(adata.uns['train_indices'])/batches)
135        print("Specified batch size required too many cells for "
136                "independent batches. Reduced batch size from "
137                f"{old_batch_size} to {batch_size}")
138
139    if 'sigma' not in adata.uns.keys():
140        n_samples = np.min((2000, adata.uns['train_indices'].shape[0]))
141        sample_range = np.arange(n_samples)
142        batch_idx = get_batches(sample_range, adata.uns['seed_obj'], 
143                                batches=batches, batch_size=batch_size)
144        sigma_indices = sample_cells(adata.uns['train_indices'], n_samples, adata.uns['seed_obj'])
145
146    # Create Arrays to store concatenated group Zs
147    # Each group of features will have a corresponding entry in each array
148    n_cols = 2*adata.uns['D']*n_pathway
149    Z_train = np.zeros((train_len, n_cols))
150    Z_test = np.zeros((test_len, n_cols))
151
152
153    # Setting kernel function 
154    match adata.uns['kernel_type'].lower():
155        case 'gaussian':
156            proj_func = gaussian_trans
157        case 'laplacian':
158            proj_func = laplacian_trans
159        case 'cauchy':
160            proj_func = cauchy_trans
161
162
163    # Loop over each of the groups and creating Z for each
164    sigma_list = list()
165    for m, group_features in enumerate(adata.uns['group_dict'].values()):
166
167        n_group_features = len(group_features)
168
169        X_train, X_test = get_group_mat(adata, n_features, group_features, 
170                                        n_group_features, process_test=True)
171        
172        if adata.uns['tfidf']:
173            X_train, X_test = tfidf_train_test(X_train, X_test)
174
175        # Data filtering, and transformation according to given data_type
176        # Will remove low variance (< 1e5) features regardless of data_type
177        # If scale_data will log scale and z-score the data
178        X_train, X_test = process_data(X_train=X_train, X_test=X_test, 
179                                       scale_data=adata.uns['scale_data'], 
180                                       return_dense=True)    
181
182        # Getting sigma
183        if 'sigma' in adata.uns.keys():
184            sigma = adata.uns['sigma'][m]
185        else:
186            sigma = est_group_sigma(adata, X_train, n_group_features, 
187                                    n_features, batch_idx=batch_idx)
188            sigma_list.append(sigma)
189            
190        assert sigma > 0, "Sigma must be more than 0"
191        train_projection, test_projection = calc_groupz(X_train, X_test, 
192                                                        adata, D, sigma, 
193                                                        proj_func)
194
195        # Store group Z in whole-Z object
196        # Preserves order to be able to extract meaningful groups
197        cos_idx, sin_idx = get_z_indices(m, D)
198
199        Z_train[0:, cos_idx] = np.cos(train_projection)
200        Z_train[0:, sin_idx] = np.sin(train_projection)
201
202        Z_test[0:, cos_idx] = np.cos(test_projection)
203        Z_test[0:, sin_idx] = np.sin(test_projection)
204
205    adata.uns['Z_train'] = Z_train*sq_i_d
206    adata.uns['Z_test'] = Z_test*sq_i_d
207
208    if 'sigma' not in adata.uns.keys():
209        adata.uns['sigma'] = np.array(sigma_list)
210
211    return adata

198def calculate_d(num_samples : int):
199    """
200    This function calculates the optimal number of dimensions for 
201    performance. See https://doi.org/10.48550/arXiv.1806.09178 for more
202    information.
203
204    Parameters
205    ----------
206    num_samples : int
207        The number of samples in the data set including both training
208        and testing sets.
209
210    Returns
211    -------
212    d : int
213        The optimal number of dimensions to run scMKL with the given 
214        data set.
215
216    Examples
217    --------
218    >>> raw_counts = scipy.sparse.load_npz('MCF7_counts.npz')
219    >>>
220    >>> num_cells = raw_counts.shape[0]
221    >>> d = scmkl.calculate_d(num_cells)
222    >>> d
223    161
224    """
225    d = int(np.sqrt(num_samples)*np.log(np.log(num_samples)))
226
227    return int(np.max([d, 100]))

264def create_adata(X: scipy.sparse._csc.csc_matrix | np.ndarray | pd.DataFrame, 
265                 feature_names: np.ndarray, cell_labels: np.ndarray, 
266                 group_dict: dict, obs_names: None | np.ndarray=None, 
267                 scale_data: bool=True, split_data: np.ndarray | None=None, 
268                 D: int | None=None, remove_features: bool=True, 
269                 train_ratio: float=0.8, distance_metric: str='euclidean', 
270                 kernel_type: str='Gaussian', random_state: int=1, 
271                 allow_multiclass: bool = False, 
272                 class_threshold: str | int = 'median',
273                 reduction: str | None = None, tfidf: bool = False):
274    """
275    Function to create an AnnData object to carry all relevant 
276    information going forward.
277
278    Parameters
279    ----------
280    X : scipy.sparse.csc_matrix | np.ndarray | pd.DataFrame
281        A data matrix of cells by features (sparse array 
282        recommended for large datasets).
283
284    feature_names : np.ndarray
285        Array of feature names corresponding with the features 
286        in `X`.
287
288    cell_labels : np.ndarray
289        A numpy array of cell phenotypes corresponding with 
290        the cells in `X`.
291
292    group_dict : dict 
293        Dictionary containing feature grouping information (i.e. 
294        `{geneset1: np.array([gene_1, gene_2, ..., gene_n]), geneset2: 
295        np.array([...]), ...}`.
296
297    obs_names : None | np.ndarray
298        The cell names corresponding to `X` to be assigned to output 
299        object `.obs_names` attribute.
300
301    scale_data : bool  
302        If `True`, data matrix is log transformed and standard 
303        scaled. 
304        
305    split_data : None | np.ndarray
306        If `None`, data will be split stratified by cell labels. 
307        Else, is an array of precalculated train/test split 
308        corresponding to samples. Can include labels for entire
309        dataset to benchmark performance or for only training
310        data to classify unknown cell types (i.e. `np.array(['train', 
311        'test', ..., 'train'])`.
312
313    D : int 
314        Number of Random Fourier Features used to calculate Z. 
315        Should be a positive integer. Higher values of D will 
316        increase classification accuracy at the cost of computation 
317        time. If set to `None`, will be calculated given number of 
318        samples. 
319    
320    remove_features : bool
321        If `True`, will remove features from `X` and `feature_names` 
322        not in `group_dict` and remove features from groupings not in 
323        `feature_names`.
324
325    train_ratio : float
326        Ratio of number of training samples to entire data set. Note:
327        if a threshold is applied, the ratio training samples may 
328        decrease depending on class balance and `class_threshold`
329        parameter if `allow_multiclass = True`.
330
331    distance_metric : str
332        The pairwise distance metric used to estimate sigma. Must
333        be one of the options used in `scipy.spatial.distance.cdist`.
334
335    kernel_type : str
336        The approximated kernel function used to calculate Zs.
337        Must be one of `'Gaussian'`, `'Laplacian'`, or `'Cauchy'`.
338
339    random_state : int
340        Integer random_state used to set the seed for 
341        reproducibilty.
342
343    allow_multiclass : bool
344        If `False`, will ensure that cell labels are binary.
345
346    class_threshold : str | int
347        Number of samples allowed in the training data for each cell
348        class in the training data. If `'median'`, the median number 
349        of cells per cell class will be the threshold for number of 
350        samples per class.
351
352    reduction: str | None
353        Choose which dimension reduction technique to perform on 
354        features within a group. 'svd' will run 
355        `sklearn.decomposition.TruncatedSVD`, 'linear' will multiply 
356        by an array of 1s down to 50 dimensions.
357        
358    tfidf: bool
359        Whether to calculate TFIDF transformation on peaks within 
360        groupings.
361        
362    Returns
363    -------
364    adata : ad.AnnData
365        AnnData with the following attributes and keys:
366
367        `adata.X` (array_like):
368            Data matrix.
369    
370        `adata.var_names` (array_like): 
371            Feature names corresponding to `adata.X`.
372
373        `adata.obs['labels']` (array_like):
374            cell classes/phenotypes from `cell_labels`.
375
376        `adata.uns['train_indices']` (array_like):
377            Indices for training data. 
378
379        `adata.uns['test_indices']` (array_like)
380            Indices for testing data.
381
382        `adata.uns['group_dict']` (dict):
383            Grouping information.
384
385        `adata.uns['seed_obj']` (np.random._generator.Generator): 
386            Seed object with seed equal to 100 * `random_state`.
387
388        `adata.uns['D']` (int):
389            Number of dimensions to scMKL with.
390
391        `adata.uns['scale_data']` (bool):
392            Whether or not data is log and z-score transformed.
393
394        `adata.uns['distance_metric']` (str): 
395            Distance metric as given.
396    
397        `adata.uns['kernel_type']` (str): 
398            Kernel function as given.
399
400        `adata.uns['svd']` (bool): 
401            Whether to calculate SVD reduction.
402
403        `adata.uns['tfidf']` (bool): 
404            Whether to calculate TF-IDF per grouping.
405
406    Examples
407    --------
408    >>> data_mat = scipy.sparse.load_npz('MCF7_RNA_matrix.npz')
409    >>> gene_names = np.load('MCF7_gene_names.pkl', allow_pickle = True)
410    >>> group_dict = np.load('hallmark_genesets.pkl', 
411    >>>                      allow_pickle = True)
412    >>> 
413    >>> adata = scmkl.create_adata(X = data_mat, 
414    ...                            feature_names = gene_names, 
415    ...                            group_dict = group_dict)
416    >>> adata
417    AnnData object with n_obs × n_vars = 1000 × 4341
418    obs: 'labels'
419    uns: 'group_dict', 'seed_obj', 'scale_data', 'D', 'kernel_type', 
420    'distance_metric', 'train_indices', 'test_indices'
421    """
422
423    assert X.shape[1] == len(feature_names), ("Different number of features "
424                                              "in X than feature names")
425    
426    if not allow_multiclass:
427        assert len(np.unique(cell_labels)) == 2, ("cell_labels must contain "
428                                                  "2 classes")
429    if D is not None:    
430        assert isinstance(D, int) and D > 0, 'D must be a positive integer'
431
432    kernel_options = ['gaussian', 'laplacian', 'cauchy']
433    assert kernel_type.lower() in kernel_options, ("Given kernel type not "
434                                                   "implemented. Gaussian, "
435                                                   "Laplacian, and Cauchy "
436                                                   "are the acceptable "
437                                                   "types.")
438
439    # Create adata object and add column names
440    adata = ad.AnnData(X)
441    adata.var_names = feature_names
442
443    if isinstance(obs_names, (np.ndarray)):
444        adata.obs_names = obs_names
445
446    filtered_feature_names, group_dict = _filter_features(feature_names, 
447                                                          group_dict)
448
449    if remove_features:
450        warnings.filterwarnings('ignore', category = ad.ImplicitModificationWarning)
451        adata = adata[:, filtered_feature_names]
452    
453    gc.collect()
454
455    # Add metadata to adata object
456    adata.uns['group_dict'] = group_dict
457    adata.uns['seed_obj'] = np.random.default_rng(100*random_state)
458    adata.uns['scale_data'] = scale_data
459    adata.uns['D'] = D if D is not None else calculate_d(adata.shape[0])
460    adata.uns['kernel_type'] = kernel_type
461    adata.uns['distance_metric'] = distance_metric
462    adata.uns['reduction'] = reduction if isinstance(reduction, str) else 'None'
463    adata.uns['tfidf'] = tfidf
464
465    if (split_data is None):
466        assert X.shape[0] == len(cell_labels), ("Different number of cells "
467                                                "than labels")
468        adata.obs['labels'] = cell_labels
469
470        if (allow_multiclass == False):
471            split = _binary_split(cell_labels, 
472                                  seed_obj = adata.uns['seed_obj'],
473                                  train_ratio = train_ratio)
474            train_indices, test_indices = split
475
476        elif (allow_multiclass == True):
477            split = _multi_class_split(cell_labels, 
478                                       seed_obj = adata.uns['seed_obj'], 
479                                       class_threshold = class_threshold,
480                                       train_ratio = train_ratio)
481            train_indices, test_indices = split
482
483        adata.uns['labeled_test'] = True
484
485    else:
486        x_eq_labs = X.shape[0] == len(cell_labels)
487        train_eq_labs = X.shape[0] == len(cell_labels)
488        assert x_eq_labs or train_eq_labs, ("Must give labels for all cells "
489                                            "or only for training cells")
490        
491        train_indices = np.where(split_data == 'train')[0]
492        test_indices = np.where(split_data == 'test')[0]
493
494        if len(cell_labels) == len(train_indices):
495
496            padded_cell_labels = np.zeros((X.shape[0])).astype('object')
497            padded_cell_labels[train_indices] = cell_labels
498            padded_cell_labels[test_indices] = 'padded_test_label'
499
500            adata.obs['labels'] = padded_cell_labels
501            adata.uns['labeled_test'] = False
502
503        elif len(cell_labels) == len(split_data):
504            adata.obs['labels'] = cell_labels
505            adata.uns['labeled_test'] = True
506
507    # Ensuring all train samples are first in adata object followed by test
508    sort_idx, train_indices, test_indices = sort_samples(train_indices, 
509                                                         test_indices)
510
511    adata = adata[sort_idx]
512
513    if not isinstance(obs_names, (np.ndarray)):
514        adata.obs = adata.obs.reset_index(drop=True)
515        adata.obs.index = adata.obs.index.astype('O')
516
517    adata.uns['train_indices'] = train_indices
518    adata.uns['test_indices'] = test_indices
519
520    if not scale_data:
521        print("WARNING: Data will not be log transformed and scaled. "
522              "To change this behavior, set scale_data to True")
523
524    return adata

164def estimate_sigma(adata: ad.AnnData,
165                   n_features: int = 5000,
166                   batches: int = 10, 
167                   batch_size: int = 100) -> ad.AnnData:
168    """
169    Calculate kernel widths to inform distribution for projection of 
170    Fourier Features. Calculates one sigma per group of features.
171
172    Parameters
173    ----------
174    adata : ad.AnnData
175        Created by `create_adata`.
176    
177    n_features : int
178        Number of random features to include when estimating sigma. 
179        Will be scaled for the whole pathway set according to a 
180        heuristic. Used for scalability.
181
182    batches : int
183        The number of batches to use for the distance calculation.
184        This will average the result of `batches` distance calculations
185        of `batch_size` randomly sampled cells. More batches will converge
186        to population distance values at the cost of scalability.
187
188    batch_size : int
189        The number of cells to include per batch for distance
190        calculations. Higher batch size will converge to population
191        distance values at the cost of scalability.
192        If `batches` * `batch_size` > # training cells,
193        `batch_size` will be reduced to `int(num training cells / 
194        batches)`.
195        
196    Returns
197    -------
198    adata : ad.AnnData
199        Key added `adata.uns['sigma']` with grouping kernel widths.
200
201    Examples
202    --------
203    >>> adata = scmkl.estimate_sigma(adata)
204    >>> adata.uns['sigma']
205    array([10.4640895 , 10.82011454,  6.16769438,  9.86156855, ...])
206    """
207    assert batch_size <= len(adata.uns['train_indices']), ("Batch size must be "
208                                                          "smaller than the "
209                                                          "training set.")
210
211    if batch_size * batches > len(adata.uns['train_indices']):
212        old_batch_size = batch_size
213        batch_size = int(len(adata.uns['train_indices'])/batches)
214        print("Specified batch size required too many cells for "
215                "independent batches. Reduced batch size from "
216                f"{old_batch_size} to {batch_size}")
217
218    if batch_size > 2000:
219        print("Warning: Batch sizes over 2000 may "
220               "result in long run-time.")
221        
222    # Getting subsample indices
223    sample_size = np.min((2000, adata.uns['train_indices'].shape[0]))
224    indices = sample_cells(adata.uns['train_indices'], sample_size=sample_size, 
225                           seed_obj=adata.uns['seed_obj'])
226
227    # Getting batch indices
228    sample_range = np.arange(sample_size)
229    batch_idx = get_batches(sample_range, adata.uns['seed_obj'], 
230                            batches, batch_size)
231
232    # Loop over every group in group_dict
233    sigma_array = np.zeros((len(adata.uns['group_dict'])))
234    for m, group_features in enumerate(adata.uns['group_dict'].values()):
235
236        n_group_features = len(group_features)
237
238        # Filtering to only features in grouping using filtered view of adata
239        X_train = get_group_mat(adata[indices], n_features=n_features, 
240                            group_features=group_features, 
241                            n_group_features=n_group_features)
242        
243        if adata.uns['tfidf']:
244            X_train = tfidf(X_train, mode='normalize')
245
246        # Data filtering, and transformation according to given data_type
247        # Will remove low variance (< 1e5) features regardless of data_type
248        # If scale_data will log scale and z-score the data
249        X_train = process_data(X_train=X_train,
250                               scale_data=adata.uns['scale_data'], 
251                               return_dense=True)    
252
253        # Estimating sigma
254        sigma = est_group_sigma(adata, X_train, n_group_features, 
255                                n_features, batch_idx=batch_idx)
256
257        sigma_array[m] = sigma
258    
259    adata.uns['sigma'] = sigma_array
260        
261    return adata

264def extract_results(results: dict, metric: str):
265    """
266    
267    """
268    summary = {'Alpha' : list(),
269               metric : list(),
270               'Number of Selected Groups' : list(),
271               'Top Group' : list()}
272    
273    alpha_list = list(results['Metrics'].keys())
274
275    # Creating summary DataFrame for each model
276    for alpha in alpha_list:
277        cur_alpha_rows = results['Norms'][alpha]
278        top_weight_rows = np.max(results['Norms'][alpha])
279        top_group_index = np.where(cur_alpha_rows == top_weight_rows)
280        num_selected = len(results['Selected_groups'][alpha])
281        top_group_name = np.array(results['Group_names'])[top_group_index]
282        
283        if 0 == num_selected:
284            top_group_name = ["No groups selected"]
285
286        summary['Alpha'].append(alpha)
287        summary[metric].append(results['Metrics'][alpha][metric])
288        summary['Number of Selected Groups'].append(num_selected)
289        summary['Top Group'].append(*top_group_name)
290    
291    return pd.DataFrame(summary)

203def find_candidates(organism: str='human', key_terms: str | list='', blacklist: str | list | bool=False):
204    """
205    Given `organism` and `key_terms`, will search for gene 
206    groupings that could fit the datasets/classification task. 
207    `blacklist` terms undesired in group names.
208
209    Parameters
210    ----------
211    organism : str
212        The species the gene grouping is for. Options are 
213        `{'Human', 'Mouse', 'Yeast', 'Fly', 'Fish', 'Worm'}`
214
215    key_terms : str | list
216        The types of cells or other specifiers the gene set is for 
217        (example: 'CD4 T').
218
219    blacklist : str | list | bool
220        Term(s) undesired in group names. Ignored unless provided.
221
222    Returns
223    -------
224    libraries : list
225        A list of gene set library names that could serve for the 
226        dataset/classification task.        
227
228    Examples
229    --------
230    >>> scmkl.find_candidates('human', key_terms=' b ')
231                                Library  No. Gene Sets
232    0                    Azimuth_2023           1241
233    1         Azimuth_Cell_Types_2021            341
234    2   Cancer_Cell_Line_Encyclopedia            967
235    3                 CellMarker_2024           1134
236    No. Key Type Matching
237    9
238    9
239    0
240    21
241    """
242    check_organism(organism)
243
244    if organism.lower() in global_lib_orgs:
245        glo = global_lib_orgs.copy()
246        glo.remove(organism)
247        other_org = glo[0]
248        libs = human_genesets
249        libs = [name for name in libs if not other_org in name.lower()]
250    else:
251        libs = gp.get_library_name(organism)
252        other_org = ''
253
254    libs = {name : gp.get_library(name, organism)
255            for name in libs}
256    
257    libs_df, _ = check_libs(libs, key_terms, blacklist, other_org)
258
259    return libs_df

527def format_adata(adata: ad.AnnData | str, cell_labels: np.ndarray | str, 
528                 group_dict: dict | str, use_raw: bool=False, 
529                 scale_data: bool=True, split_data: np.ndarray | None=None, 
530                 D: int | None=None, remove_features: bool=True, 
531                 train_ratio: float=0.8, distance_metric: str='euclidean', 
532                 kernel_type: str='Gaussian', random_state: int=1, 
533                 allow_multiclass: bool = False, 
534                 class_threshold: str | int = 'median', 
535                 reduction: str | None = None, tfidf: bool = False):
536    """
537    Function to format an `ad.AnnData` object to carry all relevant 
538    information going forward. `adata.obs_names` will be retained.
539
540    **NOTE: Information not needed for running `scmkl` will be 
541    removed.**
542
543    Parameters
544    ----------
545    adata : ad.AnnData
546        Object with data for `scmkl` to be applied to. Only requirment 
547        is that `.var_names` is correct and data matrix is in `adata.X` 
548        or `adata.raw.X`. A h5ad file can be provided as a `str` and it 
549        will be read in.
550
551    cell_labels : np.ndarray | str
552        If type `str`, the labels for `scmkl` to learn are captured 
553        from `adata.obs['cell_labels']`. Else, a `np.ndarray` of cell 
554        phenotypes corresponding with the cells in `adata.X`.
555
556    group_dict : dict | str
557        Dictionary containing feature grouping information (i.e. 
558        `{geneset1: np.array([gene_1, gene_2, ..., gene_n]), geneset2: 
559        np.array([...]), ...}`. A pickle file can be provided as a `str` 
560        and it will be read in.
561
562    obs_names : None | np.ndarray
563        The cell names corresponding to `X` to be assigned to output 
564        object `.obs_names` attribute.
565
566    use_raw : bool
567        If `False`, will use `adata.X` to create new `adata`. Else, 
568        will use `adata.raw.X`.
569
570    scale_data : bool  
571        If `True`, data matrix is log transformed and standard 
572        scaled. 
573        
574    split_data : None | np.ndarray
575        If `None`, data will be split stratified by cell labels. 
576        Else, is an array of precalculated train/test split 
577        corresponding to samples. Can include labels for entire
578        dataset to benchmark performance or for only training
579        data to classify unknown cell types (i.e. `np.array(['train', 
580        'test', ..., 'train'])`.
581
582    D : int 
583        Number of Random Fourier Features used to calculate Z. 
584        Should be a positive integer. Higher values of D will 
585        increase classification accuracy at the cost of computation 
586        time. If set to `None`, will be calculated given number of 
587        samples. 
588    
589    remove_features : bool
590        If `True`, will remove features from `X` and `feature_names` 
591        not in `group_dict` and remove features from groupings not in 
592        `feature_names`.
593
594    train_ratio : float
595        Ratio of number of training samples to entire data set. Note:
596        if a threshold is applied, the ratio training samples may 
597        decrease depending on class balance and `class_threshold`
598        parameter if `allow_multiclass = True`.
599
600    distance_metric : str
601        The pairwise distance metric used to estimate sigma. Must
602        be one of the options used in `scipy.spatial.distance.cdist`.
603
604    kernel_type : str
605        The approximated kernel function used to calculate Zs.
606        Must be one of `'Gaussian'`, `'Laplacian'`, or `'Cauchy'`.
607
608    random_state : int
609        Integer random_state used to set the seed for 
610        reproducibilty.
611
612    allow_multiclass : bool
613        If `False`, will ensure that cell labels are binary.
614
615    class_threshold : str | int
616        Number of samples allowed in the training data for each cell
617        class in the training data. If `'median'`, the median number 
618        of cells per cell class will be the threshold for number of 
619        samples per class.
620
621    reduction: str | None
622        Choose which dimension reduction technique to perform on 
623        features within a group. 'svd' will run 
624        `sklearn.decomposition.TruncatedSVD`, 'linear' will multiply 
625        by an array of 1s down to 50 dimensions.
626        
627    tfidf: bool
628        Whether to calculate TFIDF transformation on peaks within 
629        groupings.
630        
631    Returns
632    -------
633    adata : ad.AnnData
634        AnnData with the following attributes and keys:
635
636        `adata.X` (array_like):
637            Data matrix.
638    
639        `adata.var_names` (array_like): 
640            Feature names corresponding to `adata.X`.
641
642        `adata.obs['labels']` (array_like):
643            cell classes/phenotypes from `cell_labels`.
644
645        `adata.uns['train_indices']` (array_like):
646            Indices for training data. 
647
648        `adata.uns['test_indices']` (array_like)
649            Indices for testing data.
650
651        `adata.uns['group_dict']` (dict):
652            Grouping information.
653
654        `adata.uns['seed_obj']` (np.random._generator.Generator): 
655            Seed object with seed equal to 100 * `random_state`.
656
657        `adata.uns['D']` (int):
658            Number of dimensions to scMKL with.
659
660        `adata.uns['scale_data']` (bool):
661            Whether or not data is log and z-score transformed.
662
663        `adata.uns['distance_metric']` (str): 
664            Distance metric as given.
665    
666        `adata.uns['kernel_type']` (str): 
667            Kernel function as given.
668
669        `adata.uns['svd']` (bool): 
670            Whether to calculate SVD reduction.
671
672        `adata.uns['tfidf']` (bool): 
673            Whether to calculate TF-IDF per grouping.
674
675    Examples
676    --------
677    >>> adata = ad.read_h5ad('MCF7_rna.h5ad')
678    >>> group_dict = np.load('hallmark_genesets.pkl', 
679    >>>                      allow_pickle = True)
680    >>> 
681    >>> 
682    >>> # The labels in adata.obs we want to learn are 'celltypes'
683    >>> adata = scmkl.format_adata(adata, 'celltypes', 
684    ...                            group_dict)
685    >>> adata
686    AnnData object with n_obs × n_vars = 1000 × 4341
687    obs: 'labels'
688    uns: 'group_dict', 'seed_obj', 'scale_data', 'D', 'kernel_type', 
689    'distance_metric', 'train_indices', 'test_indices'
690    """
691    if str == type(adata):
692        adata = ad.read_h5ad(adata)
693
694    if str == type(group_dict):
695        group_dict = np.load(group_dict, allow_pickle=True)
696        
697    if str == type(cell_labels):
698        err_msg = f"{cell_labels} is not in `adata.obs`"
699        assert cell_labels in adata.obs.keys(), err_msg
700        cell_labels = adata.obs[cell_labels].to_numpy()
701    
702    if use_raw:
703        assert adata.raw, "`adata.raw` is empty, set `use_raw` to `False`"
704        X = adata.raw.X
705    else:
706        X = adata.X
707
708    adata = create_adata(X, adata.var_names.to_numpy().copy(), cell_labels, 
709                         group_dict, adata.obs_names.to_numpy().copy(), 
710                         scale_data, split_data, D, remove_features, 
711                         train_ratio, distance_metric, kernel_type, 
712                         random_state, allow_multiclass, class_threshold, 
713                         reduction, tfidf)
714
715    return adata

 90def format_group_names(group_names: list | pd.Series | np.ndarray, 
 91                       rm_words: list=list()):
 92    """
 93    Takes an ArrayLike object of group names and formats them.
 94
 95    Parameters
 96    ----------
 97    group_names : array_like
 98        An array of group names to format.
 99
100    rm_words : list
101        Words to remove from all group names.
102
103    Returns
104    -------
105    new_group_names : list
106        Formatted version of the input group names.
107
108    Examples
109    --------
110    >>> groups = ['HALLMARK_E2F_TARGETS', 'HALLMARK_HYPOXIA']
111    >>> new_groups = scmkl.format_group_names(groups)
112    >>> new_groups
113    ['Hallmark E2F Targets', 'Hallmark Hypoxia']
114    """
115    new_group_names = list()
116    rm_words = [word.lower() for word in rm_words]
117
118    for name in group_names:
119        new_name = list()
120        for word in re.split(r'_|\s', name):
121            if word.isalpha() and (len(word) > 3):
122                word = word.capitalize()
123            if word.lower() not in rm_words:
124                new_name.append(word)
125        new_name = ' '.join(new_name)
126        new_group_names.append(new_name)
127
128    return new_group_names

262def get_gene_groupings(lib_name: str, organism: str='human', key_terms: str | list='', 
263                       blacklist: str | list | bool=False, min_overlap: int=2,
264                      genes: list | tuple | pd.Series | np.ndarray | set=[]):
265    """
266    Takes a gene set library name and filters to groups containing 
267    element(s) in `key_terms`. If genes is provided, will 
268    ensure that there are at least `min_overlap` number of genes in 
269    each group. Resulting groups will meet all of the before-mentioned 
270    criteria if `isin_logic` is `'and'` | `'or'`.
271
272    Parameters
273    ----------
274    lib_name : str
275        The desired library name.
276
277    organism : str
278        The species the gene grouping is for. Options are 
279        `{'Human', 'Mouse', 'Yeast', 'Fly', 'Fish', 'Worm'}`.
280
281    key_terms : str | list
282        The types of cells or other specifiers the gene set is for 
283        (example: 'CD4 T').
284
285    genes : array_like
286        A vector of genes from the reference/query datasets. If not 
287        assigned, function will not filter groups based on feature 
288        overlap.
289
290    min_overlap : int
291        The minimum number of genes that must be present in a group 
292        for it to be kept. If `genes` is not given, ignored.
293
294    Returns
295    -------
296    lib : dict
297        The filtered library as a `dict` where keys are group names 
298        and keys are features.
299
300    Examples
301    --------
302    >>> dataset_feats = [
303    ...    'FUCA1', 'CLIC4', 'STMN1', 'SYF2', 'TAS1R1', 
304    ...    'NOL9', 'TAS1R3', 'SLC2A5', 'THAP3', 'IGHM', 
305    ...    'MARCKS', 'BANK1', 'TNFRSF13B', 'IGKC', 'IGHD', 
306    ...    'LINC01857', 'CD24', 'CD37', 'IGHD', 'RALGPS2'
307    ...    ]
308    >>> rna_grouping = scmkl.get_gene_groupings(
309    ...   'Azimuth_2023', key_terms=[' b ', 'b cell', 'b '], 
310    ...   genes=dataset_feats)
311    >>>
312    >>> rna_groupings.keys()
313    dict_keys(['PBMC-L1-B Cell', 'PBMC-L2-Intermediate B Cell', ...])
314    """
315    check_organism(organism)
316    
317    lib = gp.get_library(lib_name, organism)
318
319    if organism.lower() in global_lib_orgs:
320        glo = global_lib_orgs.copy()
321        glo.remove(organism)
322        other_org = glo[0]
323    else:
324        other_org = ''
325
326    group_names = list(lib.keys())
327    res = check_groups(group_names, key_terms, blacklist, other_org)
328    del res['num_groups']
329
330    # Finding groups where group name matches key_terms
331    g_summary = pd.DataFrame(res)
332
333    if key_terms:
334        kept = g_summary['key_terms_in']
335        kept_groups = g_summary['name'][kept].to_numpy()
336        g_summary = g_summary[kept]
337    else:
338        print("Not filtering with `key_terms` parameter.")
339        kept_groups = g_summary['name'].to_numpy()
340
341    if blacklist:
342        kept = ~g_summary['blacklist_in']
343        kept_groups = g_summary['name'][kept].to_numpy()
344    else:
345        print("Not filtering with `blacklist` parameter.")
346    
347    # Filtering library
348    lib = {group : lib[group] for group in kept_groups}
349
350    if 0 < len(genes):
351        del_groups = list()
352        genes = list(set(genes.copy()))
353        for group, features in lib.items():
354            overlap = np.isin(features, genes)
355            overlap = np.sum(overlap)
356            if overlap < min_overlap:
357                print(overlap, flush=True)
358                del_groups.append(group)
359
360        # Removing genes without enough overlap
361        for group in del_groups:
362            print(f'Removing {group} from grouping.')
363            del lib[group]
364
365    else:
366        print("Not checking overlap between group and dataset features.")
367
368    return lib

398def get_metrics(results: dict, include_as: bool=False) -> pd.DataFrame:
399    """
400    Takes either a single scMKL result or a dictionary where each 
401    entry cooresponds to one result. Returns a dataframe with cols 
402    ['Alpha', 'Metric', 'Value']. If `include_as == True`, another 
403    col of booleans will be added to indicate whether or not the run 
404    respective to that alpha was chosen as optimal via CV. If 
405    `include_key == True`, another column will be added with the name 
406    of the key to the respective file (only applicable with multiple 
407    results).
408
409    Parameters
410    ----------
411    results : dict | None
412        A dictionary with the results of a single run from 
413        `scmkl.run()`. Must be `None` if `rfiles is not None`.
414
415    rfiles : dict | None
416        A dictionary of results dictionaries containing multiple 
417        results from `scmkl.run()`. 
418
419    include_as : bool
420        When `True`, will add a bool col to output pd.DataFrame 
421        where rows with alphas cooresponding to alpha_star will be 
422        `True`.
423
424    Returns
425    -------
426    df : pd.DataFrame
427        A pd.DataFrame containing all of the metrics present from 
428        the runs input.
429
430    Examples
431    --------
432    >>> # For a single file
433    >>> results = scmkl.run(adata)
434    >>> metrics = scmkl.get_metrics(results = results)
435
436    >>> # For multiple runs saved in a dict
437    >>> output_dir = 'scMKL_outputs/'
438    >>> rfiles = scmkl.read_files(output_dir)
439    >>> metrics = scmkl.get_metrics(rfiles=rfiles)
440    """
441    # Checking which data is being worked with 
442    is_mult, is_many = _parse_result_type(results)
443
444    # Initiating col list with minimal columns
445    cols = ['Alpha', 'Metric', 'Value']
446
447    if include_as:
448        cols.append('Alpha Star')
449    if is_mult:
450        cols.append('Class')
451
452    if is_many:
453        cols.append('Key')
454        df = pd.DataFrame(columns = cols)
455        for key, result in results.items():
456                cur_df = parse_metrics(results = result, key = key, 
457                                        include_as = include_as)
458                df = pd.concat([df, cur_df.copy()])
459            
460    else:
461        df = parse_metrics(results = results, include_as = include_as)
462
463    return df