1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

import numpy as np
import time
from sklearn.neighbors import NearestNeighbors, LocalOutlierFactor
 
def benchmark_spatial_indices(n_samples: int, n_features: int, k: int = 20):
    """
    Benchmark different spatial indexing strategies for k-NN.
    """
    np.random.seed(42)
    X = np.random.randn(n_samples, n_features)
    
    results = {}
    
    # Brute Force
    start = time.time()
    nn_brute = NearestNeighbors(n_neighbors=k, algorithm='brute')
    nn_brute.fit(X)
    nn_brute.kneighbors(X)
    results['brute'] = time.time() - start
    
    # k-d Tree
    start = time.time()
    nn_kd = NearestNeighbors(n_neighbors=k, algorithm='kd_tree')
    nn_kd.fit(X)
    nn_kd.kneighbors(X)
    results['kd_tree'] = time.time() - start
    
    # Ball Tree
    start = time.time()
    nn_ball = NearestNeighbors(n_neighbors=k, algorithm='ball_tree')
    nn_ball.fit(X)
    nn_ball.kneighbors(X)
    results['ball_tree'] = time.time() - start
    
    return results
 
def benchmark_lof_scalability():
    """
    Benchmark LOF at different scales to show real-world timing.
    """
    print("LOF Scalability Benchmark:")
    print("=" * 60)
    print(f"{'n_samples':<12} {'d=5 (sec)':<12} {'d=20 (sec)':<12} {'d=50 (sec)':<12}")
    print("-" * 60)
    
    for n in [1000, 5000, 10000, 25000]:
        times = []
        for d in [5, 20, 50]:
            np.random.seed(42)
            X = np.random.randn(n, d)
            
            start = time.time()
            lof = LocalOutlierFactor(n_neighbors=20, algorithm='auto')
            lof.fit_predict(X)
            elapsed = time.time() - start
            times.append(elapsed)
        
        print(f"{n:<12} {times[0]:<12.3f} {times[1]:<12.3f} {times[2]:<12.3f}")
 
# Run benchmarks
print("
Spatial Index Benchmark (n=10000, k=20):")
print("-" * 40)
for d in [5, 10, 20, 50]:
    results = benchmark_spatial_indices(10000, d)
    best = min(results, key=results.get)
    print(f"d={d}: Best={best} ({results[best]:.3f}s), "
          f"Speedup over brute: {results['brute']/results[best]:.1f}x")
 
print()
benchmark_lof_scalability()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

import numpy as np
import time
 
# Note: These require additional packages
# pip install annoy hnswlib faiss-cpu
 
def demo_approximate_nn_lof():
    """
    Demonstrate approximate k-NN based LOF for large-scale data.
    """
    from sklearn.neighbors import LocalOutlierFactor
    
    class ApproximateLOF:
        """
        LOF using approximate nearest neighbors for scalability.
        
        Uses ANNOY for fast approximate k-NN queries.
        """
        
        def __init__(self, k: int = 20, n_trees: int = 50):
            self.k = k
            self.n_trees = n_trees
            
        def fit_predict(self, X: np.ndarray) -> np.ndarray:
            """
            Compute approximate LOF scores.
            """
            try:
                from annoy import AnnoyIndex
            except ImportError:
                print("ANNOY not installed, falling back to sklearn")
                lof = LocalOutlierFactor(n_neighbors=self.k)
                lof.fit_predict(X)
                return -lof.negative_outlier_factor_
            
            n_samples, n_features = X.shape
            
            # Build ANNOY index
            index = AnnoyIndex(n_features, 'euclidean')
            for i in range(n_samples):
                index.add_item(i, X[i])
            index.build(self.n_trees)
            
            # Compute approximate k-NN distances
            k_distances = np.zeros(n_samples)
            neighbors = np.zeros((n_samples, self.k), dtype=int)
            distances = np.zeros((n_samples, self.k))
            
            for i in range(n_samples):
                # Get k+1 neighbors (including self)
                nn_idx, nn_dist = index.get_nns_by_item(
                    i, self.k + 1, include_distances=True
                )
                # Exclude self (first neighbor)
                neighbors[i] = nn_idx[1:self.k+1]
                distances[i] = nn_dist[1:self.k+1]
                k_distances[i] = distances[i, -1]
            
            # Compute LRD (Local Reachability Density)
            lrd = np.zeros(n_samples)
            for i in range(n_samples):
                reach_dists = np.maximum(k_distances[neighbors[i]], distances[i])
                lrd[i] = self.k / (reach_dists.sum() + 1e-10)
            
            # Compute LOF
            lof_scores = np.zeros(n_samples)
            for i in range(n_samples):
                neighbor_lrds = lrd[neighbors[i]]
                lof_scores[i] = neighbor_lrds.mean() / (lrd[i] + 1e-10)
            
            return lof_scores
    
    # Benchmark
    np.random.seed(42)
    
    print("Approximate LOF Benchmark:")
    print("=" * 60)
    
    for n in [10000, 50000, 100000]:
        X = np.random.randn(n, 20)
        
        # Exact LOF (only for smaller datasets)
        if n <= 20000:
            start = time.time()
            lof_exact = LocalOutlierFactor(n_neighbors=20)
            lof_exact.fit_predict(X)
            exact_scores = -lof_exact.negative_outlier_factor_
            exact_time = time.time() - start
        else:
            exact_time = float('nan')
        
        # Approximate LOF
        start = time.time()
        approx_lof = ApproximateLOF(k=20, n_trees=50)
        approx_scores = approx_lof.fit_predict(X)
        approx_time = time.time() - start
        
        print(f"n={n}: Exact={exact_time:.2f}s, Approx={approx_time:.2f}s", end="")
        if n <= 20000:
            # Compute correlation between rankings
            from scipy.stats import spearmanr
            corr, _ = spearmanr(exact_scores, approx_scores)
            print(f", Rank Correlation={corr:.3f}")
        else:
            print()
 
# Example: Using FAISS for billion-scale
def faiss_example():
    """
    Example of using FAISS for very large scale k-NN.
    """
    try:
        import faiss
    except ImportError:
        print("FAISS not installed. Install with: pip install faiss-cpu")
        return
    
    np.random.seed(42)
    n, d, k = 100000, 64, 20
    
    # Generate data
    X = np.random.randn(n, d).astype('float32')
    
    # Method 1: Flat (exact, brute force with GPU support)
    start = time.time()
    index_flat = faiss.IndexFlatL2(d)
    index_flat.add(X)
    distances_flat, indices_flat = index_flat.search(X, k + 1)
    flat_time = time.time() - start
    
    # Method 2: IVF (inverted file) - approximate
    start = time.time()
    nlist = 100  # Number of clusters
    quantizer = faiss.IndexFlatL2(d)
    index_ivf = faiss.IndexIVFFlat(quantizer, d, nlist)
    index_ivf.train(X)
    index_ivf.add(X)
    index_ivf.nprobe = 10  # How many clusters to search
    distances_ivf, indices_ivf = index_ivf.search(X, k + 1)
    ivf_time = time.time() - start
    
    print(f"
FAISS Benchmark (n={n}, d={d}):")
    print(f"  Flat (exact): {flat_time:.3f}s")
    print(f"  IVF (approximate): {ivf_time:.3f}s")
    print(f"  Speedup: {flat_time/ivf_time:.1f}x")
 
# Run demos
demo_approximate_nn_lof()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

import numpy as np
from sklearn.neighbors import LocalOutlierFactor, NearestNeighbors
 
class SampledLOF:
    """
    LOF using a sampled reference set for scalability.
    
    Scores all points against a sampled reference.
    """
    
    def __init__(self, k: int = 20, sample_size: int = 10000,
                 sample_method: str = 'random'):
        """
        Parameters:
        -----------
        k : int
            Number of neighbors
        sample_size : int
            Size of reference sample
        sample_method : str
            'random', 'stratified', or 'importance'
        """
        self.k = k
        self.sample_size = sample_size
        self.sample_method = sample_method
        
    def _random_sample(self, X: np.ndarray) -> np.ndarray:
        """Uniform random sampling."""
        n = len(X)
        sample_size = min(self.sample_size, n)
        indices = np.random.choice(n, size=sample_size, replace=False)
        return indices
    
    def _importance_sample(self, X: np.ndarray) -> np.ndarray:
        """Sample inversely proportional to approximate density."""
        n = len(X)
        sample_size = min(self.sample_size, n)
        
        # Quick density estimate using random subset
        pilot_size = min(1000, n)
        pilot_idx = np.random.choice(n, size=pilot_size, replace=False)
        
        nn = NearestNeighbors(n_neighbors=min(10, pilot_size))
        nn.fit(X[pilot_idx])
        dists, _ = nn.kneighbors(X)
        approx_density = 1.0 / (dists[:, -1] + 1e-10)
        
        # Sample inversely proportional to density
        weights = 1.0 / (approx_density + 1e-10)
        weights /= weights.sum()
        
        indices = np.random.choice(n, size=sample_size, replace=False, p=weights)
        return indices
    
    def fit_predict(self, X: np.ndarray) -> np.ndarray:
        """
        Compute LOF scores using sampled reference.
        """ 
        n = len(X)
        
        # Sample reference points
        if self.sample_method == 'importance':
            sample_idx = self._importance_sample(X)
        else:
            sample_idx = self._random_sample(X)
        
        X_sample = X[sample_idx]
        
        # Compute LOF on sample
        lof_sample = LocalOutlierFactor(n_neighbors=self.k, novelty=True)
        lof_sample.fit(X_sample)
        
        # Score all points against sample
        # Note: Uses score_samples which returns negative LOF (higher = more inlier)
        # We negate to get positive scores where higher = more outlier
        scores = -lof_sample.score_samples(X)
        
        return scores
 
 
class MiniBatchLOF:
    """
    LOF computed in mini-batches for memory efficiency.
    """
    
    def __init__(self, k: int = 20, batch_size: int = 10000):
        self.k = k
        self.batch_size = batch_size
        
    def fit_predict(self, X: np.ndarray) -> np.ndarray:
        """
        Compute LOF in mini-batches.
        """
        n = len(X)
        scores = np.zeros(n)
        
        # Build global k-NN index
        nn = NearestNeighbors(n_neighbors=self.k + 1, algorithm='auto')
        nn.fit(X)
        
        # Process in batches
        for start in range(0, n, self.batch_size):
            end = min(start + self.batch_size, n)
            batch = X[start:end]
            
            # Get neighbors from global index
            distances, indices = nn.kneighbors(batch)
            
            # Compute k-distances for batch points
            k_dists = distances[:, -1]  # Last column = k-th neighbor
            
            # For each point in batch, compute LRD
            batch_lrd = np.zeros(end - start)
            for i in range(end - start):
                neighbor_idx = indices[i, 1:]  # Exclude self
                neighbor_kdists = nn.kneighbors(X[neighbor_idx], return_distance=False)
                neighbor_kdists = np.linalg.norm(
                    X[neighbor_idx] - X[neighbor_idx][:, np.newaxis], axis=2
                )[:, -1]
                
                # Simplified LRD computation
                reach_dists = np.maximum(k_dists[0], distances[i, 1:])
                batch_lrd[i] = self.k / (reach_dists.sum() + 1e-10)
            
            # Compute LOF for batch (simplified)
            # In practice, would need neighbor LRDs from previous passes
            scores[start:end] = 1.0 / (batch_lrd + 1e-10)  # Placeholder
        
        return scores
 
 
# Demonstration
np.random.seed(42)
 
# Generate large dataset with known anomalies
n_normal = 50000
n_anomalies = 100
 
X_normal = np.random.randn(n_normal, 10)
X_anomalies = np.random.randn(n_anomalies, 10) + 5  # Shifted anomalies
X = np.vstack([X_normal, X_anomalies])
true_labels = np.array([0] * n_normal + [1] * n_anomalies)
 
print("Sampling-Based Anomaly Detection:")
print("=" * 60)
print(f"Dataset: {len(X)} points, {n_anomalies} anomalies")
 
# Full LOF (for comparison, if feasible)
import time
 
# Sampled LOF variants
for method in ['random', 'importance']:
    start = time.time()
    sampled_lof = SampledLOF(k=20, sample_size=5000, sample_method=method)
    scores = sampled_lof.fit_predict(X)
    elapsed = time.time() - start
    
    # Evaluate: what fraction of top-k are true anomalies?
    top_k = 200
    top_indices = np.argsort(scores)[-top_k:]
    precision_at_k = true_labels[top_indices].mean()
    
    print(f"
{method.capitalize()} Sampling (sample_size=5000):")
    print(f"  Time: {elapsed:.2f}s")
    print(f"  Precision@{top_k}: {precision_at_k:.1%}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

import numpy as np
from collections import deque
from sklearn.neighbors import NearestNeighbors
 
class SlidingWindowLOF:
    """
    LOF with sliding window for streaming anomaly detection.
    """
    
    def __init__(self, k: int = 20, window_size: int = 5000):
        """
        Parameters:
        -----------
        k : int
            Number of neighbors
        window_size : int
            Number of recent points to keep as reference
        """
        self.k = k
        self.window_size = window_size
        self.window = deque(maxlen=window_size)
        self.nn_model = None
        self.rebuild_every = 100  # Rebuild index every N points
        self.points_since_rebuild = 0
        
    def _rebuild_index(self):
        """Rebuild k-NN index from current window."""
        if len(self.window) < self.k + 1:
            self.nn_model = None
            return
            
        X = np.array(list(self.window))
        self.nn_model = NearestNeighbors(n_neighbors=self.k + 1, algorithm='auto')
        self.nn_model.fit(X)
        self.points_since_rebuild = 0
        
    def score(self, point: np.ndarray) -> float:
        """
        Score a single point against current window.
        
        Returns:
        --------
        float : LOF-like anomaly score (higher = more anomalous)
        """
        if len(self.window) < self.k + 1:
            # Not enough history
            return 0.0
            
        if self.nn_model is None or self.points_since_rebuild >= self.rebuild_every:
            self._rebuild_index()
        
        point = np.array(point).reshape(1, -1)
        
        # Get neighbors
        distances, indices = self.nn_model.kneighbors(point)
        distances = distances[0, 1:]  # Exclude potential self-match
        indices = indices[0, 1:]
        
        k_dist_point = distances[-1]
        
        # Get neighbors' k-distances
        X_window = np.array(list(self.window))
        neighbor_kdists = []
        for idx in indices:
            if idx < len(X_window):
                d, _ = self.nn_model.kneighbors(X_window[idx:idx+1])
                neighbor_kdists.append(d[0, -1])
        
        if not neighbor_kdists:
            return 0.0
            
        neighbor_kdists = np.array(neighbor_kdists)
        
        # Compute reachability distances
        reach_dists = np.maximum(neighbor_kdists, distances[:len(neighbor_kdists)])
        
        # LRD of point
        lrd_point = self.k / (reach_dists.sum() + 1e-10)
        
        # Approximate neighbors' LRDs (using their k-distances)
        neighbor_lrds = self.k / (neighbor_kdists * self.k + 1e-10)
        
        # LOF score
        lof = neighbor_lrds.mean() / (lrd_point + 1e-10)
        
        return lof
    
    def update(self, point: np.ndarray) -> float:
        """
        Add point to window and return its anomaly score.
        """
        # Score before adding to window
        score = self.score(point)
        
        # Add to window
        self.window.append(point)
        self.points_since_rebuild += 1
        
        return score
 
 
class ExponentialDecayReference:
    """
    Reference distribution with exponential decay for concept drift.
    """
    
    def __init__(self, decay: float = 0.99, k: int = 20):
        """
        Parameters:
        -----------
        decay : float
            Weight decay factor per observation (0.99 = 1% decay per point)
        k : int
            Number of neighbors
        """
        self.decay = decay
        self.k = k
        self.points = []
        self.weights = []
        self.max_points = 10000
        
    def update(self, point: np.ndarray) -> float:
        """Add point and return its score."""
        # Decay existing weights
        self.weights = [w * self.decay for w in self.weights]
        
        # Prune low-weight points
        if len(self.weights) > 0:
            min_weight = 0.01
            keep_mask = [w > min_weight for w in self.weights]
            self.points = [p for p, k in zip(self.points, keep_mask) if k]
            self.weights = [w for w, k in zip(self.weights, keep_mask) if k]
        
        # Score current point
        if len(self.points) < self.k + 1:
            score = 0.0
        else:
            X_ref = np.array(self.points)
            weights = np.array(self.weights)
            
            # Weighted k-NN query
            dists = np.linalg.norm(X_ref - point.reshape(1, -1), axis=1)
            weighted_dists = dists / (weights + 1e-10)
            nearest = np.argsort(weighted_dists)[:self.k]
            
            # Score based on weighted distance to k-th neighbor
            score = weighted_dists[nearest[-1]]
        
        # Add new point with weight 1
        self.points.append(point)
        self.weights.append(1.0)
        
        # Limit size
        if len(self.points) > self.max_points:
            self.points = self.points[-self.max_points:]
            self.weights = self.weights[-self.max_points:]
        
        return score
 
 
# Demonstration: Streaming with concept drift
np.random.seed(42)
 
def simulate_stream_with_drift(n_points: int = 10000):
    """
    Simulate a data stream with concept drift.
    
    - First half: centered at origin
    - Second half: centered at [2, 2, ...]
    - Anomalies: random spikes throughout
    """
    d = 5
    points = []
    labels = []  # 0 = normal, 1 = anomaly
    
    for i in range(n_points):
        # Concept drift: mean shifts halfway
        if i < n_points // 2:
            mean = np.zeros(d)
        else:
            mean = np.ones(d) * 2
        
        # 1% anomalies
        if np.random.rand() < 0.01:
            point = np.random.randn(d) * 0.5 + mean + np.random.randn(d) * 5
            labels.append(1)
        else:
            point = np.random.randn(d) * 0.5 + mean
            labels.append(0)
        
        points.append(point)
    
    return np.array(points), np.array(labels)
 
# Run simulation
stream_points, stream_labels = simulate_stream_with_drift(5000)
 
detector = SlidingWindowLOF(k=20, window_size=1000)
scores = []
 
for point in stream_points:
    score = detector.update(point)
    scores.append(score)
 
scores = np.array(scores)
 
# Evaluate detection after warmup
warmup = 1000
scores_eval = scores[warmup:]
labels_eval = stream_labels[warmup:]
 
print("Streaming Anomaly Detection Results:")
print("=" * 50)
print(f"Total anomalies in evaluation period: {labels_eval.sum()}")
 
# Precision at top 100
top_100_idx = np.argsort(scores_eval)[-100:]
precision_at_100 = labels_eval[top_100_idx].mean()
print(f"Precision@100: {precision_at_100:.1%}")
 
# Detection after drift (second half)
post_drift_start = (5000 - warmup) // 2
post_drift_scores = scores_eval[post_drift_start:]
post_drift_labels = labels_eval[post_drift_start:]
top_50_post = np.argsort(post_drift_scores)[-50:]
precision_post_drift = post_drift_labels[top_50_post].mean()
print(f"Precision@50 after drift: {precision_post_drift:.1%}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170

# Pseudocode for Spark-based distributed LOF
# Requires: PySpark
 
"""
from pyspark.sql import SparkSession
from pyspark.ml.feature import VectorAssembler
from pyspark.ml.clustering import KMeans
import numpy as np
 
def distributed_lof_approximate(df, k=20, n_partitions=100):
    '''
    Approximate distributed LOF using partitioning.
    
    Strategy: Partition data spatially, compute LOF within partitions,
    then handle boundary effects.
    '''
    spark = SparkSession.builder.getOrCreate()
    
    # Step 1: Create spatial partitions using k-means
    # This groups similar points together
    kmeans = KMeans(k=n_partitions, featuresCol='features')
    model = kmeans.fit(df)
    df_with_partition = model.transform(df)
    
    # Step 2: Compute LOF within each partition
    def partition_lof(partition_df):
        '''Compute LOF for a single partition.'''
        from sklearn.neighbors import LocalOutlierFactor
        
        X = np.array([row['features'].toArray() for row in partition_df])
        if len(X) < k + 1:
            return [(row['id'], 1.0) for row in partition_df]
        
        lof = LocalOutlierFactor(n_neighbors=k)
        lof.fit_predict(X)
        scores = -lof.negative_outlier_factor_
        
        return [(row['id'], float(score)) 
                for row, score in zip(partition_df, scores)]
    
    # Apply LOF per partition
    lof_rdd = (df_with_partition
               .rdd
               .map(lambda row: (row['partition'], row))
               .groupByKey()
               .flatMap(lambda x: partition_lof(list(x[1]))))
    
    # Step 3: Handle boundary points (optional refinement)
    # Points near partition boundaries should be re-scored
    # using neighbors from adjacent partitions
    
    return lof_rdd.toDF(['id', 'lof_score'])
 
 
def distributed_sampling_lof(df, k=20, sample_fraction=0.1):
    '''
    Distributed LOF using sampling strategy.
    
    Sample reference points, broadcast to all partitions,
    score each point against the sample.
    '''
    # Sample reference points
    sample_df = df.sample(fraction=sample_fraction, seed=42)
    reference = sample_df.collect()
    reference_array = np.array([r['features'].toArray() for r in reference])
    
    # Broadcast reference to all executors
    reference_bc = spark.sparkContext.broadcast(reference_array)
    
    def score_against_reference(features):
        '''Score a point against broadcasted reference.'''
        from sklearn.neighbors import NearestNeighbors
        
        ref = reference_bc.value
        nn = NearestNeighbors(n_neighbors=k)
        nn.fit(ref)
        
        dists, _ = nn.kneighbors([features.toArray()])
        return float(dists[0].mean())  # Use mean distance as score
    
    # Apply scoring function
    from pyspark.sql.functions import udf
    from pyspark.sql.types import DoubleType
    
    score_udf = udf(score_against_reference, DoubleType())
    result = df.withColumn('anomaly_score', score_udf(df['features']))
    
    return result
"""
 
# For demonstration without Spark:
# Simulated distributed LOF logic
 
import numpy as np
from sklearn.neighbors import LocalOutlierFactor
from sklearn.cluster import MiniBatchKMeans
 
def simulated_distributed_lof(X: np.ndarray, k: int = 20, n_partitions: int = 10):
    """
    Simulate distributed LOF using partition-based approach.
    
    This demonstrates the partitioning logic that would be distributed.
    """
    n_samples = len(X)
    
    # Partition data using k-means
    kmeans = MiniBatchKMeans(n_clusters=n_partitions, random_state=42)
    partition_labels = kmeans.fit_predict(X)
    
    # Compute LOF within each partition
    scores = np.zeros(n_samples)
    
    for part in range(n_partitions):
        mask = partition_labels == part
        X_part = X[mask]
        
        if len(X_part) < k + 1:
            scores[mask] = 1.0  # Default score for small partitions
            continue
        
        lof = LocalOutlierFactor(n_neighbors=min(k, len(X_part) - 1))
        lof.fit_predict(X_part)
        scores[mask] = -lof.negative_outlier_factor_
    
    return scores, partition_labels
 
# Demonstrate
np.random.seed(42)
n = 50000
d = 10
 
X_normal = np.random.randn(n - 100, d)
X_outliers = np.random.randn(100, d) + 5
 
X = np.vstack([X_normal, X_outliers])
true_labels = np.array([0] * (n - 100) + [1] * 100)
 
import time
 
# Single-machine LOF
start = time.time()
lof_single = LocalOutlierFactor(n_neighbors=20)
lof_single.fit_predict(X)
scores_single = -lof_single.negative_outlier_factor_
single_time = time.time() - start
 
# Simulated distributed
start = time.time()
scores_dist, partitions = simulated_distributed_lof(X, k=20, n_partitions=50)
dist_time = time.time() - start
 
print("Distributed LOF Simulation:")
print("=" * 50)
print(f"Dataset size: {n} points, {d} dimensions")
print(f"
Single-machine LOF: {single_time:.2f}s")
print(f"Partitioned LOF (50 partitions): {dist_time:.2f}s")
print(f"Speedup: {single_time/dist_time:.1f}x")
 
# Evaluate
from scipy.stats import spearmanr
corr, _ = spearmanr(scores_single, scores_dist)
print(f"
Score correlation (single vs distributed): {corr:.3f}")
 
# Top-100 precision
top_100_single = np.argsort(scores_single)[-100:]
top_100_dist = np.argsort(scores_dist)[-100:]
print(f"Top-100 precision (single): {true_labels[top_100_single].mean():.1%}")
print(f"Top-100 precision (distributed): {true_labels[top_100_dist].mean():.1%}")