Kernel centres semi-stratified selection without replacement

densratio/RuLSIF.py

@@ -11,11 +11,10 @@
 """
 
 from numpy import array, asarray, asmatrix, diag, diagflat, empty, exp, inf, log, matrix, multiply, ones, power, sum
-from numpy.random import randint
 from numpy.linalg import solve
 from warnings import warn
 from .density_ratio import DensityRatio, KernelInfo
-from .helpers import guvectorize_compute, np_float, to_ndarray
+from .helpers import guvectorize_compute, np_float, semi_stratified_sample, to_ndarray
 
 
 def RuLSIF(x, y, alpha, sigma_range, lambda_range, kernel_num=100, verbose=True):
@@ -46,7 +45,7 @@ def RuLSIF(x, y, alpha, sigma_range, lambda_range, kernel_num=100, verbose=True)
     kernel_num = min(kernel_num, nx)
 
     # Randomly take a subset of x, to identify centers for the kernels.
-    centers = x[randint(nx, size=kernel_num)]
+    centers = x[semi_stratified_sample(x, size=kernel_num)]
 
     if verbose:
         print("RuLSIF starting...")

densratio/helpers.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 from numpy import array, matrix, ndarray, result_type
 
 
@@ -33,3 +35,58 @@ def to_ndarray(x):
         raise ValueError("Cannot transform to numpy.matrix.")
     else:
         return to_ndarray(array(x))
+
+
+def semi_stratified_sample(data: ndarray, size: int) -> ndarray:
+    if data.ndim > 2:
+        raise ValueError('Only single and 2d arrays are supported.')
+    if not size:
+        return np.empty(0)
+
+    data_length = data.shape[0]
+    result = np.arange(data_length, dtype=int)
+    if size == data_length:
+        np.random.shuffle(result)
+        return result
+    if size < 0:
+        raise ValueError('Sample size must be a non-negative integer number.')
+    if size > data_length:
+        raise ValueError('Sample size cannot exceed the shape of input data.')
+
+    dims = data.shape[1]
+    indexed = np.column_stack((data, result))
+    result = np.empty(0, dtype=indexed.dtype)
+
+    samples_no = size // dims
+    if samples_no:
+        percentiles = np.linspace(0., 100., num=samples_no, endpoint=False)[1:]
+
+        for d in range(dims):
+            column = indexed[..., d]
+            quantiles = np.append(column.min(), np.percentile(column, percentiles))
+            indices = []
+            i, sample_size = 0, 1
+
+            while i < samples_no:
+                left = quantiles[i]
+                i += 1
+                right = np.Inf if i == samples_no else quantiles[i]
+                try:
+                    indices.extend(np.random.choice(
+                        indexed[(left <= column) & (column < right), dims],
+                        size=sample_size,
+                        replace=False))
+                except ValueError:
+                    sample_size += 1
+                    continue
+                else:
+                    sample_size = 1
+
+            indexed = indexed[~np.isin(indexed[:, dims], indices)]
+            result = np.append(result, indices)
+
+    result = np.append(
+        result,
+        np.random.choice(indexed[..., dims], size=size - result.size, replace=False)).astype(int)
+    np.random.shuffle(result)
+    return result