Random Forests

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

import numpy as np
import time
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
import joblib
import multiprocessing
 
def analyze_parallel_scaling(X, y, n_estimators=200):
    """
    Analyze how training time scales with number of cores.
    """
    n_cores = multiprocessing.cpu_count()
    print(f"System has {n_cores} CPU cores")
    print(f"Dataset: {X.shape[0]} samples, {X.shape[1]} features")
    print(f"Training {n_estimators} trees")
    print("=" * 60)
    
    results = []
    
    for n_jobs in [1, 2, 4, 8, n_cores // 2, n_cores]:
        if n_jobs > n_cores:
            continue
            
        # Multiple runs for reliable timing
        times = []
        for _ in range(3):
            rf = RandomForestClassifier(
                n_estimators=n_estimators,
                n_jobs=n_jobs,
                random_state=42
            )
            
            start = time.time()
            rf.fit(X, y)
            times.append(time.time() - start)
        
        avg_time = np.mean(times)
        std_time = np.std(times)
        
        if n_jobs == 1:
            base_time = avg_time
            speedup = 1.0
        else:
            speedup = base_time / avg_time
        
        efficiency = speedup / n_jobs * 100
        
        results.append({
            'n_jobs': n_jobs,
            'time': avg_time,
            'speedup': speedup,
            'efficiency': efficiency
        })
        
        print(f"n_jobs={n_jobs:2} | Time: {avg_time:.2f}s (±{std_time:.2f}) | "
              f"Speedup: {speedup:.2f}x | Efficiency: {efficiency:.0f}%")
    
    return results
 
 
def optimize_joblib_backend(X, y):
    """
    Compare different joblib backends for RF training.
    """
    backends = ['loky', 'threading', 'multiprocessing']
    
    print("\nJoblib Backend Comparison")
    print("=" * 60)
    
    for backend in backends:
        try:
            with joblib.parallel_backend(backend):
                rf = RandomForestClassifier(
                    n_estimators=100,
                    n_jobs=-1,
                    random_state=42
                )
                
                start = time.time()
                rf.fit(X, y)
                elapsed = time.time() - start
                
                print(f"{backend:15} | Time: {elapsed:.2f}s")
        except Exception as e:
            print(f"{backend:15} | Error: {str(e)[:50]}")
 
 
# Example
X, y = make_classification(
    n_samples=10000, n_features=100,
    n_informative=50, random_state=42
)
 
analyze_parallel_scaling(X, y)
optimize_joblib_backend(X, y)

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

import numpy as np
import time
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
 
def benchmark_prediction_parallelism(rf, X_test):
    """
    Benchmark prediction with different parallelization settings.
    """
    print("\nPrediction Parallelism Benchmark")
    print("=" * 60)
    print(f"Test samples: {X_test.shape[0]}, Trees: {len(rf.estimators_)}")
    
    for n_jobs in [1, 2, 4, -1]:
        # Warm-up
        _ = rf.predict(X_test[:100])
        
        times = []
        for _ in range(5):
            start = time.time()
            # Note: n_jobs affects .predict() when set during fit
            predictions = rf.predict(X_test)
            times.append(time.time() - start)
        
        avg_time = np.mean(times)
        std_time = np.std(times)
        throughput = len(X_test) / avg_time
        
        jobs_str = "all" if n_jobs == -1 else str(n_jobs)
        print(f"n_jobs={jobs_str:3} | Time: {avg_time*1000:.1f}ms (±{std_time*1000:.1f}) | "
              f"Throughput: {throughput:,.0f} samples/sec")
 
 
def optimize_batch_prediction(rf, X_test, batch_sizes=[100, 1000, 10000]):
    """
    Compare batch vs single prediction performance.
    For production, batch predictions are almost always faster.
    """
    print("\nBatch Size Impact on Prediction")
    print("=" * 60)
    
    n_samples = len(X_test)
    
    for batch_size in batch_sizes:
        n_batches = (n_samples + batch_size - 1) // batch_size
        
        start = time.time()
        predictions = []
        for i in range(0, n_samples, batch_size):
            batch = X_test[i:i+batch_size]
            predictions.append(rf.predict(batch))
        predictions = np.concatenate(predictions)
        elapsed = time.time() - start
        
        throughput = n_samples / elapsed
        print(f"Batch size {batch_size:6} | {n_batches:4} batches | "
              f"Time: {elapsed*1000:.1f}ms | Throughput: {throughput:,.0f}/sec")
 
 
# Example
X, y = make_classification(
    n_samples=50000, n_features=50,
    n_informative=25, random_state=42
)
 
X_train, X_test = X[:40000], X[40000:]
y_train = y[:40000]
 
# Train with parallelism
rf = RandomForestClassifier(n_estimators=200, n_jobs=-1, random_state=42)
rf.fit(X_train, y_train)
 
benchmark_prediction_parallelism(rf, X_test)
optimize_batch_prediction(rf, X_test)

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

from pyspark.sql import SparkSession
from pyspark.ml.feature import VectorAssembler
from pyspark.ml.classification import RandomForestClassifier
from pyspark.ml.evaluation import MulticlassClassificationEvaluator
 
# Initialize Spark session
spark = SparkSession.builder \
    .appName("DistributedRandomForest") \
    .config("spark.executor.memory", "4g") \
    .config("spark.executor.cores", "4") \
    .config("spark.executor.instances", "10") \
    .getOrCreate()
 
def train_distributed_rf(train_df, feature_cols, label_col):
    """
    Train a Random Forest using Spark MLlib.
    
    Parameters:
    -----------
    train_df : Spark DataFrame
        Training data
    feature_cols : list
        Names of feature columns
    label_col : str
        Name of label column
        
    Returns:
    --------
    Fitted RandomForestModel
    """
    # Assemble features into a vector column
    assembler = VectorAssembler(
        inputCols=feature_cols,
        outputCol="features"
    )
    train_df = assembler.transform(train_df)
    
    # Configure Random Forest
    rf = RandomForestClassifier(
        labelCol=label_col,
        featuresCol="features",
        numTrees=200,            # Equivalent to n_estimators
        maxDepth=20,             # Tree depth limit
        maxBins=32,              # Max bins for discretizing features
        featureSubsetStrategy="sqrt",  # Equivalent to max_features='sqrt'
        impurity="gini",
        seed=42
    )
    
    # Train model
    model = rf.fit(train_df)
    
    print(f"Trained Random Forest with {model.getNumTrees} trees")
    print(f"Feature importances: {model.featureImportances}")
    
    return model
 
 
def distributed_prediction(model, test_df, feature_cols):
    """
    Make predictions using the distributed model.
    """
    # Assemble features
    assembler = VectorAssembler(
        inputCols=feature_cols,
        outputCol="features"
    )
    test_df = assembler.transform(test_df)
    
    # Predict
    predictions = model.transform(test_df)
    
    return predictions
 
 
def evaluate_model(predictions, label_col):
    """
    Evaluate the model using Spark evaluators.
    """
    evaluator = MulticlassClassificationEvaluator(
        labelCol=label_col,
        predictionCol="prediction",
        metricName="accuracy"
    )
    
    accuracy = evaluator.evaluate(predictions)
    print(f"Test Accuracy: {accuracy:.4f}")
    
    return accuracy
 
 
# Example usage (assuming data is loaded into Spark DataFrame)
"""
# Load data
train_df = spark.read.parquet("s3://bucket/train_data.parquet")
test_df = spark.read.parquet("s3://bucket/test_data.parquet")
 
# Define columns
feature_cols = [f"feature_{i}" for i in range(100)]
label_col = "label"
 
# Train
model = train_distributed_rf(train_df, feature_cols, label_col)
 
# Predict
predictions = distributed_prediction(model, test_df, feature_cols)
 
# Evaluate
evaluate_model(predictions, label_col)
 
# Save model
model.write().overwrite().save("s3://bucket/rf_model")
"""

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

import dask.array as da
import dask.dataframe as dd
from dask.distributed import Client
from dask_ml.wrappers import ParallelPostFit
from sklearn.ensemble import RandomForestClassifier
import numpy as np
 
def setup_dask_cluster():
    """
    Set up a Dask distributed cluster.
    For local testing, this creates a local cluster.
    For production, connect to an existing cluster.
    """
    # Local cluster (uses all cores)
    client = Client(n_workers=4, threads_per_worker=2)
    
    # For existing cluster:
    # client = Client("scheduler_address:8786")
    
    print(f"Dashboard: {client.dashboard_link}")
    return client
 
 
def train_with_parallel_prediction(X_train, y_train, X_test):
    """
    Train sklearn RF normally, use Dask for parallel prediction.
    Best for: Large prediction sets, normal training data size.
    """
    # Train with regular sklearn (it's already parallel)
    rf = RandomForestClassifier(
        n_estimators=200,
        n_jobs=-1,
        random_state=42
    )
    rf.fit(X_train, y_train)
    
    # Wrap for parallel prediction across Dask chunks
    parallel_rf = ParallelPostFit(rf)
    
    # Convert test data to Dask array (chunked)
    X_test_da = da.from_array(X_test, chunks=(10000, -1))
    
    # Predictions run in parallel across chunks
    predictions = parallel_rf.predict(X_test_da)
    
    # Compute (triggers execution)
    return predictions.compute()
 
 
def distributed_rf_training(X_da, y_da):
    """
    Fully distributed RF training using Dask.
    Uses ensemble of RF models trained on different chunks.
    
    Note: This trains DIFFERENT models on different data chunks,
    then averages predictions. Different from true distributed RF.
    """
    from dask_ml.ensemble import BlockwiseVotingClassifier
    
    # Create a blockwise ensemble
    # Each block trains its own RF, predictions are averaged
    ensemble = BlockwiseVotingClassifier(
        estimator=RandomForestClassifier(n_estimators=50, random_state=42),
        classes=np.array([0, 1]),
        voting='soft'
    )
    
    # Train on Dask arrays (each chunk trains independently)
    ensemble.fit(X_da, y_da)
    
    return ensemble
 
 
def incremental_rf_training(X_chunks, y_chunks):
    """
    Train RF incrementally on data chunks.
    Uses warm_start to add trees as more data is seen.
    """
    rf = RandomForestClassifier(
        n_estimators=50,  # Trees per chunk
        warm_start=True,
        random_state=42
    )
    
    trees_per_chunk = 50
    
    for i, (X_chunk, y_chunk) in enumerate(zip(X_chunks, y_chunks)):
        print(f"Processing chunk {i+1}...")
        
        # Increase target number of trees
        rf.n_estimators = (i + 1) * trees_per_chunk
        
        # Train (adds new trees, keeps old ones)
        rf.fit(X_chunk, y_chunk)
        
        print(f"  Trees so far: {len(rf.estimators_)}")
    
    return rf
 
 
# Example usage
"""
# Set up Dask cluster
client = setup_dask_cluster()
 
# Load data as Dask array/dataframe
X_train = da.from_zarr("data/X_train.zarr")
y_train = da.from_zarr("data/y_train.zarr")
X_test = da.from_zarr("data/X_test.zarr")
 
# Option 1: Train locally, predict in parallel
predictions = train_with_parallel_prediction(
    X_train.compute(), y_train.compute(), X_test.compute()
)
 
# Option 2: Fully distributed (ensemble of RFs)
ensemble = distributed_rf_training(X_train, y_train)
predictions = ensemble.predict(X_test).compute()
 
# Clean up
client.close()
"""

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

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import pickle
import time
 
def find_optimal_tree_count(X, y, max_trees=500, tolerance=0.001):
    """
    Find minimum number of trees that achieves near-optimal accuracy.
    
    tolerance: acceptable accuracy drop from maximum
    """
    tree_counts = [10, 25, 50, 100, 150, 200, 300, 400, 500]
    tree_counts = [t for t in tree_counts if t <= max_trees]
    
    results = []
    for n_trees in tree_counts:
        rf = RandomForestClassifier(n_estimators=n_trees, random_state=42)
        score = cross_val_score(rf, X, y, cv=5).mean()
        results.append((n_trees, score))
        print(f"Trees: {n_trees:3} | CV Accuracy: {score:.4f}")
    
    max_score = max(r[1] for r in results)
    threshold = max_score - tolerance
    
    # Find minimum trees meeting threshold
    for n_trees, score in results:
        if score >= threshold:
            print(f"\nOptimal: {n_trees} trees (achieves {score:.4f}, "
                  f"within {tolerance} of max {max_score:.4f})")
            return n_trees
    
    return tree_counts[-1]
 
 
def prune_trees_by_importance(rf, keep_fraction=0.5):
    """
    Keep only the most important trees based on OOB contribution.
    
    Note: This is a simple heuristic. Trees that perform best on
    OOB samples are kept.
    """
    n_trees = len(rf.estimators_)
    n_keep = max(1, int(n_trees * keep_fraction))
    
    # Estimate tree quality (would need OOB predictions stored)
    # This is a placeholder - actual implementation needs OOB tracking
    print(f"Pruning from {n_trees} to {n_keep} trees")
    
    # Simple approach: keep first n_keep trees (assumes random = uncorrelated)
    pruned_estimators = rf.estimators_[:n_keep]
    
    # Create new RF with subset
    from sklearn.base import clone
    pruned_rf = clone(rf)
    pruned_rf.n_estimators = n_keep
    pruned_rf.estimators_ = pruned_estimators
    
    return pruned_rf
 
 
def quantize_thresholds(rf, precision=6):
    """
    Reduce memory by quantizing split thresholds.
    Reduces model size with minimal accuracy impact.
    """
    for tree in rf.estimators_:
        tree_struct = tree.tree_
        # Round thresholds to fewer decimal places
        tree_struct.threshold[:] = np.round(
            tree_struct.threshold, decimals=precision
        )
    
    return rf
 
 
def measure_model_size(rf):
    """
    Measure serialized model size.
    """
    serialized = pickle.dumps(rf)
    size_mb = len(serialized) / (1024 * 1024)
    print(f"Model size: {size_mb:.2f} MB")
    print(f"Trees: {len(rf.estimators_)}")
    avg_depth = np.mean([t.tree_.max_depth for t in rf.estimators_])
    print(f"Average tree depth: {avg_depth:.1f}")
    return size_mb
 
 
def benchmark_inference_speed(rf, X_test, n_runs=100):
    """
    Benchmark prediction latency.
    """
    # Warm up
    _ = rf.predict(X_test[:10])
    
    # Single sample latency
    times = []
    for i in range(n_runs):
        start = time.time()
        _ = rf.predict(X_test[i:i+1])
        times.append(time.time() - start)
    
    print(f"Single sample latency:")
    print(f"  Mean: {np.mean(times)*1000:.2f}ms")
    print(f"  p50: {np.percentile(times, 50)*1000:.2f}ms")
    print(f"  p99: {np.percentile(times, 99)*1000:.2f}ms")
    
    # Batch latency
    batch_start = time.time()
    _ = rf.predict(X_test)
    batch_time = time.time() - batch_start
    
    print(f"\nBatch prediction ({len(X_test)} samples):")
    print(f"  Total: {batch_time*1000:.2f}ms")
    print(f"  Per sample: {batch_time/len(X_test)*1000:.3f}ms")
 
 
# Example
from sklearn.datasets import make_classification
 
X, y = make_classification(
    n_samples=5000, n_features=50,
    n_informative=25, random_state=42
)
 
X_train, X_test = X[:4000], X[4000:]
y_train = y[:4000]
 
# Train full model
rf = RandomForestClassifier(n_estimators=200, random_state=42)
rf.fit(X_train, y_train)
 
print("Full Model:")
measure_model_size(rf)
benchmark_inference_speed(rf, X_test)
 
# Find optimal tree count
print("\nOptimizing tree count:")
optimal_trees = find_optimal_tree_count(X_train, y_train)

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

from fastapi import FastAPI, HTTPException
from pydantic import BaseModel
import numpy as np
import pickle
from typing import List
import asyncio
from concurrent.futures import ThreadPoolExecutor
 
# Load model at startup
with open("rf_model.pkl", "rb") as f:
    MODEL = pickle.load(f)
 
# Thread pool for CPU-bound prediction
EXECUTOR = ThreadPoolExecutor(max_workers=4)
 
app = FastAPI(title="Random Forest Inference Service")
 
 
class PredictionRequest(BaseModel):
    features: List[List[float]]  # Batch of feature vectors
 
 
class PredictionResponse(BaseModel):
    predictions: List[int]
    probabilities: List[List[float]]
 
 
def predict_sync(features: np.ndarray):
    """Synchronous prediction (runs in thread pool)."""
    predictions = MODEL.predict(features)
    probabilities = MODEL.predict_proba(features)
    return predictions, probabilities
 
 
@app.post("/predict", response_model=PredictionResponse)
async def predict(request: PredictionRequest):
    """
    Make predictions on a batch of samples.
    
    Runs prediction in thread pool to avoid blocking async event loop.
    """
    try:
        features = np.array(request.features)
        
        # Validate input shape
        expected_features = MODEL.n_features_in_
        if features.shape[1] != expected_features:
            raise HTTPException(
                status_code=400,
                detail=f"Expected {expected_features} features, got {features.shape[1]}"
            )
        
        # Run prediction in thread pool (CPU-bound)
        loop = asyncio.get_event_loop()
        predictions, probabilities = await loop.run_in_executor(
            EXECUTOR, predict_sync, features
        )
        
        return PredictionResponse(
            predictions=predictions.tolist(),
            probabilities=probabilities.tolist()
        )
        
    except Exception as e:
        raise HTTPException(status_code=500, detail=str(e))
 
 
@app.get("/health")
async def health():
    """Health check endpoint."""
    return {
        "status": "healthy",
        "model_trees": len(MODEL.estimators_),
        "model_features": MODEL.n_features_in_
    }
 
 
@app.get("/model/info")
async def model_info():
    """Return model metadata."""
    return {
        "n_estimators": len(MODEL.estimators_),
        "n_features": MODEL.n_features_in_,
        "n_classes": len(MODEL.classes_),
        "classes": MODEL.classes_.tolist(),
        "max_depth": max(t.tree_.max_depth for t in MODEL.estimators_)
    }
 
 
# Run with: uvicorn production_serving:app --host 0.0.0.0 --port 8000