scmkl.estimate_sigma

 1import numpy as np
 2import scipy
 3
 4from scmkl.calculate_z import _process_data
 5
 6
 7def estimate_sigma(adata, n_features = 5000):
 8    '''
 9    Calculate kernel widths to inform distribution for projection of 
10    Fourier Features. Calculates one sigma per group of features.
11
12    Parameters
13    ----------
14    **adata** : *AnnData* 
15        > Created by `create_adata`.
16    
17    **n_features** : *int*  
18        > Number of random features to include when estimating sigma. 
19        Will be scaled for the whole pathway set according to a 
20        heuristic. Used for scalability.
21    
22    Returns
23    -------
24    **adata** : *AnnData*
25        > Key added `adata.uns['sigma']`.
26
27    Examples
28    --------
29    >>> adata = scmkl.estimate_sigma(adata)
30    >>> adata.uns['sigma']
31    array([10.4640895 , 10.82011454,  6.16769438,  9.86156855, ...])
32    '''
33 
34    sigma_list = []
35
36    # Loop over every group in group_dict
37    for group_features in adata.uns['group_dict'].values():
38
39        # Select only features in that group and downsample for scalability
40        num_group_features = len(group_features)
41        group_array = np.array(list(group_features))
42        n_feats = min([n_features, num_group_features])
43        group_features = adata.uns['seed_obj'].choice(group_array, n_feats, 
44                                                      replace = False) 
45
46        # Use on the train data to estimate sigma
47        X_train = adata[adata.uns['train_indices'], group_features].X
48        X_train = _process_data(X_train = X_train, 
49                                scale_data = adata.uns['scale_data'], 
50                                return_dense = True)
51        
52        # Sample cells for scalability
53        sample_idx = np.arange(X_train.shape[0])
54        n_samples = np.min((2000, X_train.shape[0]))
55        distance_indices = adata.uns['seed_obj'].choice(sample_idx, n_samples)
56
57        # Calculate Distance Matrix with specified metric
58        sigma = scipy.spatial.distance.cdist(X_train[distance_indices,:], 
59                                             X_train[distance_indices,:], 
60                                             adata.uns['distance_metric'])
61        sigma = np.mean(sigma)
62
63        # sigma = 0 is numerically unusable in later steps
64        # Using such a small sigma will result in wide distribution, and 
65        # typically a non-predictive Z
66        if sigma == 0:
67            sigma += 1e-5
68
69        if n_features < num_group_features:
70            # Heuristic we calculated to account for fewer features used in 
71            # distance calculation
72            sigma = sigma * num_group_features / n_features 
73
74        sigma_list.append(sigma)
75    
76    adata.uns['sigma'] = np.array(sigma_list)
77        
78    return adata