|
| 1 | +import logging |
| 2 | + |
| 3 | +import numpy as np |
| 4 | +from scipy.sparse import csr_matrix |
| 5 | + |
| 6 | +from tsne import _tsne |
| 7 | +from tsne.nearest_neighbors import KDTree, NNDescent, KNNIndex |
| 8 | + |
| 9 | +log = logging.getLogger(__name__) |
| 10 | + |
| 11 | + |
| 12 | +class Affinities: |
| 13 | + """Compute the affinities among some initial data and new data. |
| 14 | +
|
| 15 | + tSNE takes as input an affinity matrix P, and does not really care about |
| 16 | + the space in which the original data points lie. This means we are not |
| 17 | + limited to problems with numeric matrices (although that is the most common |
| 18 | + use-case) but can also optimize graph layouts. |
| 19 | +
|
| 20 | + We use perplexity, as defined by Van der Maaten in the original paper as a |
| 21 | + continuous analogue to the number of neighbor affinities we want to |
| 22 | + preserve during optimization. |
| 23 | +
|
| 24 | + """ |
| 25 | + def __init__(self, perplexity=30): |
| 26 | + self.perplexity = perplexity |
| 27 | + self.P = None |
| 28 | + |
| 29 | + def to_new(self, data, perplexity=None, return_distances=False): |
| 30 | + """Compute the affinities of new data points to the existing ones. |
| 31 | +
|
| 32 | + This is especially useful for `transform` where we need the conditional |
| 33 | + probabilities from the existing to the new data. |
| 34 | +
|
| 35 | + """ |
| 36 | + |
| 37 | + |
| 38 | +class NearestNeighborAffinities(Affinities): |
| 39 | + """Compute affinities using the nearest neighbors defined by perplexity.""" |
| 40 | + def __init__(self, data, perplexity=30, method='approx', metric='euclidean', |
| 41 | + symmetrize=True, n_jobs=1): |
| 42 | + self.n_samples = data.shape[0] |
| 43 | + |
| 44 | + perplexity = self.check_perplexity(perplexity) |
| 45 | + k_neighbors = min(self.n_samples - 1, int(3 * perplexity)) |
| 46 | + |
| 47 | + # Support shortcuts for built-in nearest neighbor methods |
| 48 | + methods = {'exact': KDTree, 'approx': NNDescent} |
| 49 | + if isinstance(method, KNNIndex): |
| 50 | + knn_index = method |
| 51 | + |
| 52 | + elif method not in methods: |
| 53 | + raise ValueError('Unrecognized nearest neighbor algorithm `%s`. ' |
| 54 | + 'Please choose one of the supported methods or ' |
| 55 | + 'provide a valid `KNNIndex` instance.') |
| 56 | + else: |
| 57 | + knn_index = methods[method](metric=metric, n_jobs=n_jobs) |
| 58 | + |
| 59 | + knn_index.build(data) |
| 60 | + neighbors, distances = knn_index.query_train(data, k=k_neighbors) |
| 61 | + |
| 62 | + # Store the results on the object |
| 63 | + self.perplexity = perplexity |
| 64 | + self.knn_index = knn_index |
| 65 | + self.P = joint_probabilities_nn( |
| 66 | + neighbors, distances, perplexity, symmetrize=symmetrize, n_jobs=n_jobs) |
| 67 | + |
| 68 | + self.n_jobs = n_jobs |
| 69 | + |
| 70 | + def to_new(self, data, perplexity=None, return_distances=False): |
| 71 | + perplexity = perplexity or self.perplexity |
| 72 | + perplexity = self.check_perplexity(perplexity) |
| 73 | + k_neighbors = min(self.n_samples - 1, int(3 * perplexity)) |
| 74 | + |
| 75 | + neighbors, distances = self.knn_index.query(data, k_neighbors) |
| 76 | + |
| 77 | + P = joint_probabilities_nn( |
| 78 | + neighbors, distances, perplexity, symmetrize=False, |
| 79 | + n_reference_samples=self.n_samples, n_jobs=self.n_jobs, |
| 80 | + ) |
| 81 | + |
| 82 | + if return_distances: |
| 83 | + return P, neighbors, distances |
| 84 | + |
| 85 | + return P |
| 86 | + |
| 87 | + def check_perplexity(self, perplexity): |
| 88 | + """Check for valid perplexity value.""" |
| 89 | + if self.n_samples - 1 < 3 * perplexity: |
| 90 | + old_perplexity, perplexity = perplexity, (self.n_samples - 1) / 3 |
| 91 | + log.warning('Perplexity value %d is too high. Using perplexity %.2f' % |
| 92 | + (old_perplexity, perplexity)) |
| 93 | + |
| 94 | + return perplexity |
| 95 | + |
| 96 | + |
| 97 | +class GraphAffinities(Affinities): |
| 98 | + def __init__(self, data, use_directed=True, use_weights=True): |
| 99 | + super().__init__() |
| 100 | + |
| 101 | + def to_new(self, data): |
| 102 | + pass |
| 103 | + |
| 104 | + |
| 105 | +def joint_probabilities_nn(neighbors, distances, perplexity, symmetrize=True, |
| 106 | + n_reference_samples=None, n_jobs=1): |
| 107 | + """Compute the conditional probability matrix P_{j|i}. |
| 108 | +
|
| 109 | + This method computes an approximation to P using the nearest neighbors. |
| 110 | +
|
| 111 | + Parameters |
| 112 | + ---------- |
| 113 | + neighbors : np.ndarray |
| 114 | + A `n_samples * k_neighbors` matrix containing the indices to each |
| 115 | + points' nearest neighbors in descending order. |
| 116 | + distances : np.ndarray |
| 117 | + A `n_samples * k_neighbors` matrix containing the distances to the |
| 118 | + neighbors at indices defined in the neighbors parameter. |
| 119 | + perplexity : double |
| 120 | + The desired perplexity of the probability distribution. |
| 121 | + symmetrize : bool |
| 122 | + Whether to symmetrize the probability matrix or not. Symmetrizing is |
| 123 | + used for typical t-SNE, but does not make sense when embedding new data |
| 124 | + into an existing embedding. |
| 125 | + n_reference_samples : int |
| 126 | + The number of samples in the existing (reference) embedding. Needed to |
| 127 | + properly construct the sparse P matrix. |
| 128 | + n_jobs : int |
| 129 | + Number of threads. |
| 130 | +
|
| 131 | + Returns |
| 132 | + ------- |
| 133 | + csr_matrix |
| 134 | + A `n_samples * n_reference_samples` matrix containing the probabilities |
| 135 | + that a new sample would appear as a neighbor of a reference point. |
| 136 | +
|
| 137 | + """ |
| 138 | + n_samples, k_neighbors = distances.shape |
| 139 | + |
| 140 | + if n_reference_samples is None: |
| 141 | + n_reference_samples = n_samples |
| 142 | + |
| 143 | + # Compute asymmetric pairwise input similarities |
| 144 | + conditional_P = _tsne.compute_gaussian_perplexity( |
| 145 | + distances, perplexity, num_threads=n_jobs) |
| 146 | + conditional_P = np.asarray(conditional_P) |
| 147 | + |
| 148 | + P = csr_matrix((conditional_P.ravel(), neighbors.ravel(), |
| 149 | + range(0, n_samples * k_neighbors + 1, k_neighbors)), |
| 150 | + shape=(n_samples, n_reference_samples)) |
| 151 | + |
| 152 | + # Symmetrize the probability matrix |
| 153 | + if symmetrize: |
| 154 | + P = (P + P.T) / 2 |
| 155 | + |
| 156 | + # Convert weights to probabilities using pair-wise normalization scheme |
| 157 | + P /= np.sum(P) |
| 158 | + |
| 159 | + return P |
0 commit comments