24  Lab 24: AI as a Microscope for the Brain

Open In Colab

24.1 Learning Objectives

1. Understand core concepts from Chapter 24
2. Implement key algorithms and techniques
3. Apply methods to practical problems
4. Analyze results and draw insights
5. Connect to broader NeuroAI themes

24.2 Prerequisites

• Reading: Chapter 24: AI for Neuroscience Discovery
• Libraries: NumPy, Matplotlib, relevant frameworks
• Concepts: ML for neural data

24.3 Setup

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
plt.rcParams['figure.figsize'] = (12, 8)

24.4 Part 1: Fundamentals

Core concepts and theory from AI for Neuroscience Discovery.

## Example implementation
def example_function():
    '''Placeholder for main implementation.'''
    pass

24.5 Part 2: Hands-On Implementation

Practical exercises exploring ML for neural data.

24.6 Part 3: Analysis

Analyze and visualize results.

24.7 Exercises

24.7.1 Exercise 1

Implement the basic algorithm.

24.7.2 Exercise 2

Explore parameter variations.

24.7.3 Exercise 3

Apply to real or simulated data.

24.7.4 Exercise 4

Compare different approaches.

24.7.5 Exercise 5

Discuss NeuroAI connections.

---

24.8 Challenge Problems

24.8.1 Challenge 1

Advanced implementation.

24.8.2 Challenge 2

Reproduce research findings.

24.8.3 Challenge 3

Novel application.

---

Open In Colab

24.9 Discussion Questions

24.9.1 Question 1

How does this relate to neuroscience?

24.9.2 Question 2

What are the AI implications?

24.9.3 Question 3

Future research directions?

---

24.10 Summary

Explored ML for neural data through theory and practice.

Key Takeaways: - Core concepts mastered - Practical implementation skills developed - NeuroAI connections understood

---

Open In Colab

24.11 Code Examples from Chapter 24 (Part 1)

24.11.1 Code Block 1: Neural Speech Decoder (CNN+LSTM for ECoG Data)

# Neural speech decoder using CNN+LSTM for ECoG data
import torch
import torch.nn as nn

class SpeechDecoder(nn.Module):
    def __init__(self, n_channels, n_classes):
        super().__init__()
        # A convolutional and recurrent architecture to process spatio-temporal ECoG data
        self.conv1 = nn.Conv1d(n_channels, 128, kernel_size=3)
        self.lstm = nn.LSTM(128, 256, batch_first=True)
        self.fc = nn.Linear(256, n_classes) # n_classes = vocabulary size

    def forward(self, ecog_data):
        # ecog_data shape: (batch, channels, time)
        x = torch.relu(self.conv1(ecog_data))
        x = x.transpose(1, 2) # Switch to (batch, time, features) for LSTM
        lstm_out, _ = self.lstm(x)
        # Take the output of the last time step
        last_step_out = lstm_out[:, -1, :]
        logits = self.fc(last_step_out)
        return logits

24.11.2 Code Block 2: AlphaFold-style Contact Prediction Map

# AlphaFold-style contact prediction map visualization
import numpy as np
import matplotlib.pyplot as plt

# Simulate a contact map showing predicted amino acid interactions
np.random.seed(42)
sequence_length = 100

# Create a realistic-looking contact map
contact_map = np.zeros((sequence_length, sequence_length))

# Add diagonal (sequential contacts)
for i in range(sequence_length-1):
    contact_map[i, i+1] = 0.9
    contact_map[i+1, i] = 0.9

# Add secondary structure patterns (alpha helices and beta sheets)
# Alpha helix pattern (i, i+3 and i, i+4 contacts)
for i in range(20, 40):
    if i+3 < sequence_length:
        contact_map[i, i+3] = 0.7
        contact_map[i+3, i] = 0.7
    if i+4 < sequence_length:
        contact_map[i, i+4] = 0.6
        contact_map[i+4, i] = 0.6

# Beta sheet pattern (long-range contacts)
for i in range(60, 70):
    for j in range(80, 90):
        if abs(i-j) > 10:
            contact_map[i, j] = 0.5 + 0.2 * np.random.rand()
            contact_map[j, i] = contact_map[i, j]

# Add some noise
contact_map += np.random.rand(sequence_length, sequence_length) * 0.1
contact_map = np.clip(contact_map, 0, 1)

fig, ax = plt.subplots(1, 1, figsize=(8, 8))
im = ax.imshow(contact_map, cmap='RdPu', origin='lower')
ax.set_xlabel('Amino Acid Position', fontsize=12)
ax.set_ylabel('Amino Acid Position', fontsize=12)
ax.set_title('AlphaFold-style Contact Prediction Map', fontsize=14, fontweight='bold')
plt.colorbar(im, ax=ax, label='Contact Probability')
plt.tight_layout()
plt.show()

print("Contact map shows predicted spatial proximity between amino acids")
print("Diagonal: Sequential contacts | Patterns: Secondary structures (helices, sheets)")

24.11.3 Code Block 3: AlphaFold Impact on Structural Biology

# AlphaFold impact on structural biology - growth visualization
import pandas as pd

# Historical data on protein structures
years = np.array([2000, 2005, 2010, 2015, 2018, 2020, 2021, 2022, 2023])
experimental = np.array([15000, 30000, 65000, 110000, 145000, 170000, 180000, 190000, 200000])
predicted = np.array([0, 0, 0, 0, 0, 0, 365000, 550000, 800000])  # AlphaFold launched 2020

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(years, experimental, 'o-', linewidth=2, markersize=8,
        label='Experimental Structures (PDB)', color='#0066cc')
ax.plot(years, predicted, 's-', linewidth=2, markersize=8,
        label='AI-Predicted Structures', color='#cc0000')
ax.axvline(x=2020, color='gray', linestyle='--', alpha=0.7, label='AlphaFold 2 Release')
ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Number of Protein Structures', fontsize=12)
ax.set_title('The AlphaFold Revolution in Structural Biology', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print("AlphaFold expanded known protein structures by orders of magnitude")
print("Impact: From ~170K experimental to ~800K+ total structures in 3 years")

24.11.4 Code Block 4: Calcium Imaging Cell Detection (Suite2p)

# Calcium imaging cell detection and ROI segmentation (Suite2p)
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
from matplotlib.patches import Ellipse

np.random.seed(42)

# Create synthetic calcium imaging frame
img_size = 256
n_cells = 45

# Generate background (neuropil)
background = np.random.rand(img_size, img_size) * 0.3

# Add cell bodies
cell_image = background.copy()
cell_centers = []

for i in range(n_cells):
    x = np.random.randint(15, img_size-15)
    y = np.random.randint(15, img_size-15)
    cell_centers.append((x, y))

    # Create Gaussian cell body
    y_grid, x_grid = np.ogrid[:img_size, :img_size]
    radius = np.random.uniform(3, 6)
    intensity = np.random.uniform(0.7, 1.0)
    mask = ((x_grid - x)**2 + (y_grid - y)**2) <= radius**2
    cell_image[mask] += intensity

# Add noise
cell_image += np.random.randn(img_size, img_size) * 0.1
cell_image = np.clip(cell_image, 0, 1.5)

# Create "detected" ROIs
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Panel 1: Raw calcium imaging data
axes[0].imshow(cell_image, cmap='gray', interpolation='nearest')
axes[0].set_title('Raw Calcium Imaging Frame', fontsize=12, fontweight='bold')
axes[0].axis('off')

# Panel 2: Detected cells (Suite2p style)
axes[1].imshow(cell_image, cmap='gray', interpolation='nearest')
for x, y in cell_centers:
    circle = plt.Circle((x, y), 5, color='#cc0000', fill=False, linewidth=1.5, alpha=0.7)
    axes[1].add_patch(circle)
axes[1].set_title('Suite2p: Automated Cell Detection', fontsize=12, fontweight='bold')
axes[1].axis('off')

# Panel 3: Extracted ROIs colored by identity
roi_image = np.zeros((img_size, img_size, 3))
roi_image[:, :, 0] = cell_image  # Red channel: raw image
colors = plt.cm.Set3(np.linspace(0, 1, n_cells))

for idx, (x, y) in enumerate(cell_centers):
    y_grid, x_grid = np.ogrid[:img_size, :img_size]
    radius = 5
    mask = ((x_grid - x)**2 + (y_grid - y)**2) <= radius**2
    for c in range(3):
        roi_image[mask, c] = colors[idx, c]

axes[2].imshow(roi_image)
axes[2].set_title('Segmented ROIs (Color-coded)', fontsize=12, fontweight='bold')
axes[2].axis('off')

plt.tight_layout()
plt.show()

print(f"Detected {n_cells} neurons automatically")
print("Suite2p reduces analysis time from days to minutes")

24.11.5 Code Block 5: Spike Inference from Calcium Traces

# Spike inference from calcium traces with performance metrics
time = np.linspace(0, 10, 1000)  # 10 seconds
dt = time[1] - time[0]

# Generate ground truth spikes (Poisson process)
np.random.seed(42)
spike_rate = 5  # Hz
spike_prob = spike_rate * dt
spikes = np.random.rand(len(time)) < spike_prob

# Simulate calcium dynamics (convolution with exponential kernel)
tau_rise = 0.05  # 50 ms rise
tau_decay = 0.3  # 300 ms decay
t_kernel = np.arange(0, 1, dt)
kernel = (1 - np.exp(-t_kernel/tau_rise)) * np.exp(-t_kernel/tau_decay)
kernel /= kernel.sum()

calcium = np.convolve(spikes.astype(float), kernel, mode='same')
calcium += np.random.randn(len(time)) * 0.05  # Add noise

# Simulate deep learning spike inference
# Simple threshold-based inference for demonstration
calcium_smooth = ndimage.gaussian_filter1d(calcium, sigma=5)
inferred_spikes = np.zeros_like(spikes)
threshold = np.percentile(calcium_smooth, 75)
inferred_spikes[calcium_smooth > threshold] = 1

fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

# Ground truth spikes
axes[0].eventplot([time[spikes]], lineoffsets=0.5, linelengths=0.8, colors='black', linewidths=2)
axes[0].set_ylabel('Spikes', fontsize=11)
axes[0].set_title('Ground Truth: Action Potentials', fontsize=12, fontweight='bold')
axes[0].set_ylim([0, 1])
axes[0].set_yticks([])

# Calcium trace
axes[1].plot(time, calcium, linewidth=1.5, color='#0066cc')
axes[1].set_ylabel('ΔF/F', fontsize=11)
axes[1].set_title('Measured: Calcium Fluorescence', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

# Inferred spikes
axes[2].eventplot([time[inferred_spikes.astype(bool)]], lineoffsets=0.5,
                   linelengths=0.8, colors='#cc0000', linewidths=2)
axes[2].set_ylabel('Inferred', fontsize=11)
axes[2].set_xlabel('Time (s)', fontsize=11)
axes[2].set_title('Deep Learning Output: Inferred Spikes', fontsize=12, fontweight='bold')
axes[2].set_ylim([0, 1])
axes[2].set_yticks([])

plt.tight_layout()
plt.show()

# Calculate performance metrics
true_positive = np.sum(spikes & inferred_spikes.astype(bool))
false_positive = np.sum(~spikes & inferred_spikes.astype(bool))
false_negative = np.sum(spikes & ~inferred_spikes.astype(bool))
precision = true_positive / (true_positive + false_positive) if (true_positive + false_positive) > 0 else 0
recall = true_positive / (true_positive + false_negative) if (true_positive + false_negative) > 0 else 0

print(f"Spike Inference Performance:")
print(f"  Precision: {precision:.2f} | Recall: {recall:.2f}")
print(f"  True Positive: {true_positive} | False Positive: {false_positive} | False Negative: {false_negative}")

---

24.12 Code Examples from Chapter 24 (Part 3: Lines 1000-1500)

24.12.1 Code Block 1: CNN Decoder Training Visualization

# CNN decoder training visualization - loss curves, accuracy comparison, confusion matrix
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: Training curves
ax = axes[0, 0]
ax.plot(train_losses, label='Training Loss', linewidth=2, color='#cc0000')
ax.set_xlabel('Epoch', fontsize=11)
ax.set_ylabel('Loss', fontsize=11)
ax.set_title('CNN Training Loss', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 2: Accuracy curves
ax = axes[0, 1]
ax.plot(train_accs, label='Train Accuracy', linewidth=2, color='#0066cc')
ax.plot(test_accs, label='Test Accuracy', linewidth=2, color='#cc0000')
ax.axhline(y=linear_accuracy, color='gray', linestyle='--', label='Linear Baseline', linewidth=2)
ax.set_xlabel('Epoch', fontsize=11)
ax.set_ylabel('Accuracy', fontsize=11)
ax.set_title('CNN Learning Curves', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 3: CNN confusion matrix
model.eval()
with torch.no_grad():
    outputs = model(X_test_pt)
    _, y_pred_cnn = torch.max(outputs.data, 1)
    y_pred_cnn = y_pred_cnn.numpy()

cm_cnn = confusion_matrix(y_test, y_pred_cnn)
ax = axes[1, 0]
sns.heatmap(cm_cnn, annot=True, fmt='d', cmap='Greens', ax=ax, cbar_kws={'label': 'Count'})
ax.set_xlabel('Predicted Stimulus', fontsize=11)
ax.set_ylabel('True Stimulus', fontsize=11)
ax.set_title(f'CNN Decoder Confusion Matrix - Accuracy: {test_accs[-1]:.3f}',
             fontsize=12, fontweight='bold')

# Panel 4: Comparison
ax = axes[1, 1]
methods = ['Linear - Decoder', 'CNN - Decoder']
accuracies = [linear_accuracy, test_accs[-1]]
colors = ['#0066cc', '#cc0000']
bars = ax.bar(methods, accuracies, color=colors, alpha=0.7, edgecolor='black', linewidth=2)

# Add value labels on bars
for bar, acc in zip(bars, accuracies):
    height = bar.get_height()
    ax.text(bar.get_x() + bar.get_width()/2., height,
            f'{acc:.3f}',
            ha='center', va='bottom', fontsize=12, fontweight='bold')

ax.set_ylabel('Decoding Accuracy', fontsize=11)
ax.set_title('Decoder Performance Comparison', fontsize=12, fontweight='bold')
ax.set_ylim([0, 1.0])
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print(); print("Key Findings:")
print(f"• Linear decoder: {linear_accuracy:.1%} accuracy")
print(f"• CNN decoder: {test_accs[-1]:.1%} accuracy")
print(f"• Deep learning captures nonlinear temporal dynamics")
print(f"• Improvement is most pronounced for confusable stimuli")

24.12.2 Code Block 2: CNN Feature Visualization

# CNN feature visualization - learned filters, activation patterns, neuron importance, latent space
# Extract learned convolutional filters
conv1_weights = model.conv1.weight.data.cpu().numpy()  # Shape: (64, n_neurons, 3)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: First layer filters (sample)
ax = axes[0, 0]
n_filters_show = 16
filters_to_show = conv1_weights[:n_filters_show, :, 1]  # Middle of kernel
im = ax.imshow(filters_to_show, aspect='auto', cmap='RdBu_r', interpolation='nearest')
ax.set_xlabel('Neuron #', fontsize=11)
ax.set_ylabel('Filter #', fontsize=11)
ax.set_title('Learned Conv1 Filters (Spatial Patterns)', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax, label='Weight')

# Panel 2: Activation patterns for different stimuli
ax = axes[0, 1]
model.eval()
# Get activations for example stimuli
example_stimuli = [0, 2, 4, 6]
with torch.no_grad():
    for stim in example_stimuli:
        stim_data = X_test_pt[y_test_pt == stim][:5]  # First 5 trials
        # Forward through first conv layer
        activations = torch.relu(model.bn1(model.conv1(stim_data)))
        mean_activation = activations.mean(dim=(0, 2)).numpy()  # Average over trials and time
        ax.plot(mean_activation, label=f'Stimulus {stim}', linewidth=2)
ax.set_xlabel('Filter #', fontsize=11)
ax.set_ylabel('Mean Activation', fontsize=11)
ax.set_title('Filter Activations by Stimulus', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 3: Decoding confidence across neurons
ax = axes[1, 0]
# Measure which neurons contribute most to decoding
neuron_importance = np.zeros(n_neurons)
for neuron in range(n_neurons):
    # Shuffle this neuron's activity
    X_test_shuffled = X_test_pt.clone()
    shuffle_idx = torch.randperm(X_test_shuffled.shape[0])
    X_test_shuffled[:, neuron, :] = X_test_shuffled[shuffle_idx, neuron, :]

    # Measure accuracy drop
    with torch.no_grad():
        outputs_shuffled = model(X_test_shuffled)
        _, pred_shuffled = torch.max(outputs_shuffled.data, 1)
        acc_shuffled = (pred_shuffled == y_test_pt).sum().item() / len(y_test_pt)

    neuron_importance[neuron] = test_accs[-1] - acc_shuffled  # Drop in accuracy

ax.bar(range(n_neurons), neuron_importance, color='#9966cc', alpha=0.7)
ax.set_xlabel('Neuron #', fontsize=11)
ax.set_ylabel('Importance - (Accuracy Drop When Shuffled)', fontsize=11)
ax.set_title('Neuron Importance for Decoding', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3, axis='y')

# Panel 4: Decision boundary in latent space
ax = axes[1, 1]
# Get latent representations (before final classification layer)
model.eval()
with torch.no_grad():
    # Forward to second-to-last layer
    x = X_test_pt
    x = torch.relu(model.bn1(model.conv1(x)))
    x = model.pool(torch.relu(model.bn2(model.conv2(x))))
    x = x.view(x.size(0), -1)
    latent = torch.relu(model.fc1(x)).numpy()

# PCA on latent space
pca_latent = PCA(n_components=2)
latent_pca = pca_latent.fit_transform(latent)

colors = plt.cm.rainbow(np.linspace(0, 1, n_stimuli))
for stim in range(n_stimuli):
    mask = y_test == stim
    ax.scatter(latent_pca[mask, 0], latent_pca[mask, 1], c=[colors[stim]],
               label=f'Stim {stim}', alpha=0.6, s=30)
ax.set_xlabel(f'PC1 ({pca_latent.explained_variance_ratio_[0]:.1%})', fontsize=11)
ax.set_ylabel(f'PC2 ({pca_latent.explained_variance_ratio_[1]:.1%})', fontsize=11)
ax.set_title('CNN Latent Space (PCA)', fontsize=12, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(); print("Visualization Insights:")
print("• Conv filters learn to detect specific neuron-time patterns")
print("• Different stimuli activate different filter combinations")
print("• Neurons with broad tuning are most important for decoding")
print("• Latent space shows better stimulus separation than raw activity")

24.12.3 Code Block 3: Simulating Neural Trajectories for Decision-Making

# Simulating neural trajectories for decision-making task with latent dynamics
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(42)

# Simulate a decision-making task with neural trajectories
n_neurons = 80
n_trials = 200
n_timepoints = 50
n_conditions = 2  # Two choice conditions

def generate_decision_trajectories(n_neurons, n_trials_per_condition, n_timepoints, condition):
    """
    Generate neural trajectories during a decision-making task.
    Neural activity follows a low-dimensional trajectory from stimulus onset to choice.
    """
    trajectories = np.zeros((n_trials_per_condition, n_neurons, n_timepoints))

    # Define a low-dimensional trajectory in "latent space"
    # Three dimensions: stimulus encoding, decision process, motor preparation
    for trial in range(n_trials_per_condition):
        # Latent trajectory parameters
        stimulus_strength = np.random.uniform(0.5, 1.5)
        decision_speed = np.random.uniform(0.8, 1.2)
        motor_noise = np.random.randn()

        # Time-varying latent states
        t_normalized = np.linspace(0, 1, n_timepoints)

        # Dimension 1: Stimulus encoding (rises then plateaus)
        latent_1 = stimulus_strength * (1 - np.exp(-5 * t_normalized))

        # Dimension 2: Decision variable (S-shaped accumulation, condition-dependent)
        if condition == 0:
            latent_2 = 1.0 / (1 + np.exp(-10 * (t_normalized - 0.5)))
        else:
            latent_2 = -1.0 / (1 + np.exp(-10 * (t_normalized - 0.5)))

        # Dimension 3: Motor preparation (late ramp)
        latent_3 = np.maximum(0, decision_speed * (t_normalized - 0.6)) + 0.1 * motor_noise

        # Create random projection from latent space to neural space
        # Each neuron is a random linear combination of latent dimensions plus noise
        if trial == 0:  # Create mixing matrix once per condition
            mixing_matrix = np.random.randn(n_neurons, 3)
            # Normalize
            mixing_matrix = mixing_matrix / np.linalg.norm(mixing_matrix, axis=1, keepdims=True)

        # Project latent trajectory to neural space
        latent_trajectory = np.vstack([latent_1, latent_2, latent_3])  # (3, n_timepoints)
        neural_trajectory = mixing_matrix @ latent_trajectory  # (n_neurons, n_timepoints)

        # Add noise and ensure non-negativity (firing rates)
        neural_trajectory += np.random.randn(n_neurons, n_timepoints) * 0.3
        neural_trajectory = np.maximum(neural_trajectory, 0)

        # Add baseline firing rate
        neural_trajectory += 2.0

        trajectories[trial, :, :] = neural_trajectory

    return trajectories

# Generate data for both conditions
print("Generating neural trajectories for decision-making task...")
trials_per_condition = n_trials // n_conditions

trajectories_cond0 = generate_decision_trajectories(n_neurons, trials_per_condition, n_timepoints, condition=0)
trajectories_cond1 = generate_decision_trajectories(n_neurons, trials_per_condition, n_timepoints, condition=1)

# Combine
X_traj = np.vstack([trajectories_cond0, trajectories_cond1])
conditions = np.array([0]*trials_per_condition + [1]*trials_per_condition)

print(f"Data shape: {X_traj.shape} (trials, neurons, time)")

# Visualize raw neural data
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: Single trial activity
ax = axes[0, 0]
trial_idx = 0
im = ax.imshow(X_traj[trial_idx, :, :], aspect='auto', cmap='hot', interpolation='nearest')
ax.set_xlabel('Time Step', fontsize=11)
ax.set_ylabel('Neuron #', fontsize=11)
ax.set_title(f'Single Trial Neural Activity (Condition {conditions[trial_idx]})',
             fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax, label='Firing Rate')

# Panel 2: Average activity across trials for each condition
ax = axes[0, 1]
avg_cond0 = trajectories_cond0.mean(axis=0)
avg_cond1 = trajectories_cond1.mean(axis=0)
im = ax.imshow(avg_cond0 - avg_cond1, aspect='auto', cmap='RdBu_r',
               interpolation='nearest', vmin=-1, vmax=1)
ax.set_xlabel('Time Step', fontsize=11)
ax.set_ylabel('Neuron #', fontsize=11)
ax.set_title('Differential Activity: Condition 0 - Condition 1', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax, label='ΔFiring Rate')

# Panel 3: Example neurons over time
ax = axes[1, 0]
example_neurons = [5, 20, 35, 50]
time = np.arange(n_timepoints)
for neuron in example_neurons:
    mean_0 = trajectories_cond0[:, neuron, :].mean(axis=0)
    mean_1 = trajectories_cond1[:, neuron, :].mean(axis=0)
    ax.plot(time, mean_0, label=f'Neuron {neuron} (Cond 0)', linewidth=2)
    ax.plot(time, mean_1, '--', label=f'Neuron {neuron} (Cond 1)', linewidth=2)
ax.set_xlabel('Time Step', fontsize=11)
ax.set_ylabel('Firing Rate', fontsize=11)
ax.set_title('Example Neuron Trajectories', fontsize=12, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)

# Panel 4: Population variance over time
ax = axes[1, 1]
variance_cond0 = trajectories_cond0.var(axis=0).mean(axis=0)  # Avg variance across neurons
variance_cond1 = trajectories_cond1.var(axis=0).mean(axis=0)
ax.plot(time, variance_cond0, label='Condition 0', linewidth=2, color='#0066cc')
ax.plot(time, variance_cond1, label='Condition 1', linewidth=2, color='#cc0000')
ax.set_xlabel('Time Step', fontsize=11)
ax.set_ylabel('Population Variance', fontsize=11)
ax.set_title('Neural Variability Over Time', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

24.12.4 Code Block 4: PCA for Neural Manifolds

# PCA for neural manifolds - explained variance, trajectories in PC space, 3D manifold
# Reshape data for PCA: (trials * time, neurons)
X_flat = X_traj.reshape(-1, n_neurons)

# Apply PCA
pca = PCA(n_components=10)
X_pca = pca.fit_transform(X_flat)

# Reshape back to (trials, time, components)
X_pca = X_pca.reshape(n_trials, n_timepoints, -1)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: Explained variance
ax = axes[0, 0]
explained_var = pca.explained_variance_ratio_
cumulative_var = np.cumsum(explained_var)
ax.bar(range(1, 11), explained_var, alpha=0.7, color='#0066cc', label='Individual')
ax.plot(range(1, 11), cumulative_var, 'o-', color='#cc0000', linewidth=2,
        markersize=8, label='Cumulative')
ax.set_xlabel('Principal Component', fontsize=11)
ax.set_ylabel('Explained Variance Ratio', fontsize=11)
ax.set_title('PCA: Explained Variance', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

print(f"Top 3 PCs explain {cumulative_var[2]:.1%} of variance")

# Panel 2: Neural trajectories in PC space (2D)
ax = axes[0, 1]
colors_cond0 = plt.cm.Blues(np.linspace(0.3, 1, n_timepoints))
colors_cond1 = plt.cm.Reds(np.linspace(0.3, 1, n_timepoints))

# Plot a few example trials
n_trials_plot = 5
for trial in range(n_trials_plot):
    if conditions[trial] == 0:
        ax.plot(X_pca[trial, :, 0], X_pca[trial, :, 1], 'o-',
                color='#0066cc', alpha=0.4, linewidth=1, markersize=2)
    else:
        trial_idx = trials_per_condition + trial
        ax.plot(X_pca[trial_idx, :, 0], X_pca[trial_idx, :, 1], 's-',
                color='#cc0000', alpha=0.4, linewidth=1, markersize=2)

# Mark start and end
for cond in [0, 1]:
    mask = conditions == cond
    start_mean = X_pca[mask, 0, :2].mean(axis=0)
    end_mean = X_pca[mask, -1, :2].mean(axis=0)
    color = '#0066cc' if cond == 0 else '#cc0000'
    ax.plot(start_mean[0], start_mean[1], 'o', markersize=15, color=color,
            markeredgecolor='black', markeredgewidth=2, label=f'Cond {cond} Start')
    ax.plot(end_mean[0], end_mean[1], '*', markersize=20, color=color,
            markeredgecolor='black', markeredgewidth=2, label=f'Cond {cond} End')

ax.set_xlabel(f'PC1 ({explained_var[0]:.1%})', fontsize=11)
ax.set_ylabel(f'PC2 ({explained_var[1]:.1%})', fontsize=11)
ax.set_title('Neural Trajectories in PC Space', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Panel 3: 3D trajectories
ax = fig.add_subplot(2, 2, 3, projection='3d')
for trial in range(n_trials_plot):
    if conditions[trial] == 0:
        ax.plot(X_pca[trial, :, 0], X_pca[trial, :, 1], X_pca[trial, :, 2],
                color='#0066cc', alpha=0.4, linewidth=1.5)
    else:
        trial_idx = trials_per_condition + trial
        ax.plot(X_pca[trial_idx, :, 0], X_pca[trial_idx, :, 1], X_pca[trial_idx, :, 2],
                color='#cc0000', alpha=0.4, linewidth=1.5)

ax.set_xlabel(f'PC1 ({explained_var[0]:.1%})', fontsize=10)
ax.set_ylabel(f'PC2 ({explained_var[1]:.1%})', fontsize=10)
ax.set_zlabel(f'PC3 ({explained_var[2]:.1%})', fontsize=10)
ax.set_title('3D Neural Manifold (PCA)', fontsize=12, fontweight='bold')

# Panel 4: PC loadings (which neurons contribute to each PC)
ax = axes[1, 1]
loadings = pca.components_[:3, :]  # Top 3 PCs
im = ax.imshow(loadings, aspect='auto', cmap='RdBu_r', interpolation='nearest')
ax.set_xlabel('Neuron #', fontsize=11)
ax.set_ylabel('Principal Component', fontsize=11)
ax.set_title('PC Loadings: Neuron Contributions', fontsize=12, fontweight='bold')
ax.set_yticks([0, 1, 2])
ax.set_yticklabels(['PC1', 'PC2', 'PC3'])
plt.colorbar(im, ax=ax, label='Loading')

plt.tight_layout()
plt.show()

---

Open In Colab

24.13 Code Examples from Chapter 24 (Part 2: Lines 500-1000)

24.13.1 Code Block 1: Connectomics Project Scale and AI Contribution Visualization

# Visualize the scale and accuracy of connectomics projects
import matplotlib.pyplot as plt
import numpy as np

projects = ['C. elegans (1986)', 'Larval Fly (2020)', 'Adult Fly (2020)',
            'FlyWire (2022)', 'Mouse Retina', 'Mouse Cortex']
neurons = [302, 3000, 25000, 130000, 80000, 100000]
synapses = [7000, 500000, 20000000, 54500000, 30000000, 1000000000]
year_completed = [1986, 2020, 2020, 2022, 2019, 2025]
ai_percentage = [0, 50, 90, 95, 85, 98]  # Estimated AI contribution

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: Neurons and synapses mapped
ax1 = axes[0, 0]
x_pos = np.arange(len(projects))
width = 0.35
ax1.bar(x_pos - width/2, np.array(neurons)/1000, width, label='Neurons (thousands)', color='#0066cc')
ax1_twin = ax1.twinx()
ax1_twin.bar(x_pos + width/2, np.array(synapses)/1000000, width, label='Synapses (millions)', color='#cc0000')
ax1.set_xlabel('Project', fontsize=11)
ax1.set_ylabel('Neurons (thousands)', fontsize=11, color='#0066cc')
ax1_twin.set_ylabel('Synapses (millions)', fontsize=11, color='#cc0000')
ax1.set_title('Scale of Connectomics Projects', fontsize=12, fontweight='bold')
ax1.set_xticks(x_pos)
ax1.set_xticklabels(projects, rotation=45, ha='right')
ax1.set_yscale('log')
ax1_twin.set_yscale('log')
ax1.grid(True, alpha=0.3, axis='y')

# Panel 2: AI contribution over time
ax2 = axes[0, 1]
ax2.plot(year_completed, ai_percentage, 'o-', linewidth=2, markersize=10, color='#9966cc')
for i, project in enumerate(projects):
    ax2.annotate(project.split(' - ')[0], (year_completed[i], ai_percentage[i]),
                 textcoords="offset points", xytext=(0,10), ha='center', fontsize=8)
ax2.set_xlabel('Year', fontsize=11)
ax2.set_ylabel('AI Contribution (%)', fontsize=11)
ax2.set_title('Increasing Role of AI in Connectomics', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_ylim([0, 105])

# Panel 3: Time to complete (estimated)
ax3 = axes[1, 0]
manual_years = [30, 500, 3000, 15000, 2000, 100000]  # Estimated if done manually
actual_years = [15, 3, 3, 2, 4, 5]  # Actual project duration
speedup = np.array(manual_years) / np.array(actual_years)

x_pos = np.arange(len(projects))
bars = ax3.barh(x_pos, speedup, color='#cc9900', alpha=0.7)
ax3.set_yticks(x_pos)
ax3.set_yticklabels(projects)
ax3.set_xlabel('Speedup Factor (log scale)', fontsize=11)
ax3.set_title('AI Acceleration of Connectomics', fontsize=12, fontweight='bold')
ax3.set_xscale('log')
ax3.grid(True, alpha=0.3, axis='x')

# Add annotations
for i, (bar, speed) in enumerate(zip(bars, speedup)):
    ax3.text(speed, bar.get_y() + bar.get_height()/2, f'{speed:.0f}x',
             va='center', ha='left', fontweight='bold', fontsize=9)

# Panel 4: Data size and processing time
ax4 = axes[1, 1]
data_sizes = [0.001, 0.1, 26, 95, 10, 1000]  # TB
processing_times = [5, 2, 1.5, 1, 1.2, 0.8]  # years

scatter = ax4.scatter(data_sizes, processing_times, s=np.array(neurons)/200,
                      c=ai_percentage, cmap='RdYlGn', alpha=0.6, edgecolors='black', linewidth=1.5)
for i, project in enumerate(projects):
    ax4.annotate(project.split(' - ')[0], (data_sizes[i], processing_times[i]),
                 textcoords="offset points", xytext=(5,5), ha='left', fontsize=8)

ax4.set_xlabel('Dataset Size (TB, log scale)', fontsize=11)
ax4.set_ylabel('Processing Time (years)', fontsize=11)
ax4.set_title('Efficiency Gains: Data Size vs Processing Time', fontsize=12, fontweight='bold')
ax4.set_xscale('log')
ax4.grid(True, alpha=0.3)
cbar = plt.colorbar(scatter, ax=ax4, label='AI Contribution (%)')

plt.tight_layout()
plt.show()

print("Key Insights:")
print("- AI enables 100-10,000x speedup in connectome reconstruction")
print("- Modern projects process PB-scale data in years, not centuries")
print("- Hybrid AI+human approaches achieve highest accuracy")

24.13.2 Code Block 2: Simulating Neural Population Responses with Tuning Curves

import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, confusion_matrix
import seaborn as sns

np.random.seed(42)
torch.manual_seed(42)

# Simulation parameters
n_neurons = 100  # Population size
n_stimuli = 8    # Number of different stimuli (e.g., oriented gratings)
n_trials = 50    # Trials per stimulus
n_timepoints = 20  # Time bins per trial

# Create tuning curves for neurons
# Each neuron has a preferred stimulus orientation
preferred_stimuli = np.random.randint(0, n_stimuli, n_neurons)
tuning_width = 2.0  # Width of tuning curves

def generate_neural_responses(stimulus_id, n_trials, n_neurons, n_timepoints):
    """Generate realistic neural population responses to a stimulus."""
    responses = np.zeros((n_trials, n_neurons, n_timepoints))

    for trial in range(n_trials):
        for neuron in range(n_neurons):
            # Tuning curve: Gaussian centered on preferred stimulus
            pref = preferred_stimuli[neuron]
            # Circular distance on stimulus space
            dist = min(abs(stimulus_id - pref), n_stimuli - abs(stimulus_id - pref))
            tuning_response = np.exp(-(dist**2) / (2 * tuning_width**2))

            # Base firing rate + tuned response + noise
            base_rate = 5.0  # Hz
            max_response = 30.0  # Hz
            mean_rate = base_rate + max_response * tuning_response

            # Generate Poisson spike counts with temporal dynamics
            for t in range(n_timepoints):
                # Add temporal modulation (onset transient)
                temporal_mod = 1.0 + 0.5 * np.exp(-t / 5.0)
                rate = mean_rate * temporal_mod
                responses[trial, neuron, t] = np.random.poisson(rate * 0.05)  # 50ms bins

    return responses

# Generate dataset
X = []  # Neural activity
y = []  # Stimulus labels

print("Generating neural population data...")
for stim in range(n_stimuli):
    responses = generate_neural_responses(stim, n_trials, n_neurons, n_timepoints)
    X.append(responses)
    y.extend([stim] * n_trials)

X = np.vstack(X)  # Shape: (n_stimuli * n_trials, n_neurons, n_timepoints)
y = np.array(y)

print(f"Dataset shape: {X.shape}")
print(f"Labels shape: {y.shape}")

# Visualize population response to one stimulus
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: Raster plot for one trial
trial_idx = 0
ax = axes[0, 0]
for neuron in range(min(50, n_neurons)):
    spike_times = np.where(X[trial_idx, neuron, :] > 0)[0]
    spike_counts = X[trial_idx, neuron, spike_times]
    for t, count in zip(spike_times, spike_counts):
        ax.plot([t]*int(count), [neuron]*int(count), 'k.', markersize=2)
ax.set_xlabel('Time Bin', fontsize=11)
ax.set_ylabel('Neuron #', fontsize=11)
ax.set_title(f'Raster Plot: Stimulus {y[trial_idx]}', fontsize=12, fontweight='bold')

# Panel 2: Population response heatmap
ax = axes[0, 1]
trial_avg = X[y == 0].mean(axis=0)  # Average across trials for stimulus 0
im = ax.imshow(trial_avg, aspect='auto', cmap='hot', interpolation='nearest')
ax.set_xlabel('Time Bin', fontsize=11)
ax.set_ylabel('Neuron #', fontsize=11)
ax.set_title('Population Activity Heatmap (Stimulus 0)', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax, label='Firing Rate (Hz)')

# Panel 3: Tuning curves of example neurons
ax = axes[1, 0]
example_neurons = [0, 25, 50, 75]
for neuron in example_neurons:
    mean_responses = []
    for stim in range(n_stimuli):
        stim_trials = X[y == stim, neuron, :].mean(axis=1)  # Average across time
        mean_responses.append(stim_trials.mean())  # Average across trials
    ax.plot(range(n_stimuli), mean_responses, 'o-', label=f'Neuron {neuron}', linewidth=2)
ax.set_xlabel('Stimulus ID', fontsize=11)
ax.set_ylabel('Mean Firing Rate (Hz)', fontsize=11)
ax.set_title('Tuning Curves: Example Neurons', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Panel 4: Population decoding potential
ax = axes[1, 1]
# Show how well each stimulus is separated in neural space (first 2 PCs)
from sklearn.decomposition import PCA
X_flat = X.reshape(X.shape[0], -1)  # Flatten for PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_flat)

colors = plt.cm.rainbow(np.linspace(0, 1, n_stimuli))
for stim in range(n_stimuli):
    mask = y == stim
    ax.scatter(X_pca[mask, 0], X_pca[mask, 1], c=[colors[stim]],
               label=f'Stim {stim}', alpha=0.6, s=30)
ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%} var)', fontsize=11)
ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%} var)', fontsize=11)
ax.set_title('Neural Population Space (PCA)', fontsize=12, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"PC1+PC2 explain {(pca.explained_variance_ratio_[:2].sum()):.1%} of variance")

24.13.3 Code Block 3: Linear Decoder Baseline Using Logistic Regression

# Prepare data for sklearn
X_flat = X.reshape(X.shape[0], -1)  # Flatten time and neurons
X_train, X_test, y_train, y_test = train_test_split(X_flat, y, test_size=0.3, random_state=42)

# Train linear decoder
print("Training linear decoder (Logistic Regression)...")
linear_decoder = LogisticRegression(max_iter=1000, multi_class='multinomial', random_state=42)
linear_decoder.fit(X_train, y_train)

# Evaluate
y_pred_linear = linear_decoder.predict(X_test)
linear_accuracy = accuracy_score(y_test, y_pred_linear)
print(f"Linear Decoder Accuracy: {linear_accuracy:.3f}")

# Confusion matrix
cm_linear = confusion_matrix(y_test, y_pred_linear)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot confusion matrix
ax = axes[0]
sns.heatmap(cm_linear, annot=True, fmt='d', cmap='Blues', ax=ax, cbar_kws={'label': 'Count'})
ax.set_xlabel('Predicted Stimulus', fontsize=11)
ax.set_ylabel('True Stimulus', fontsize=11)
ax.set_title(f'Linear Decoder Confusion Matrix - Accuracy: {linear_accuracy:.3f}',
             fontsize=12, fontweight='bold')

# Plot decoder weights
ax = axes[1]
weights = linear_decoder.coef_  # Shape: (n_stimuli, n_features)
weights_reshaped = weights.reshape(n_stimuli, n_neurons, n_timepoints)
mean_weights = np.abs(weights_reshaped).mean(axis=2)  # Average across time

im = ax.imshow(mean_weights, aspect='auto', cmap='RdBu_r', interpolation='nearest')
ax.set_xlabel('Neuron #', fontsize=11)
ax.set_ylabel('Stimulus #', fontsize=11)
ax.set_title('Linear Decoder Weights - (Absolute Mean Across Time)', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax, label='Weight Magnitude')

plt.tight_layout()
plt.show()

24.13.4 Code Block 4: CNN-Based Deep Learning Decoder for Neural Data

# Prepare data for PyTorch
X_train_pt = torch.FloatTensor(X_train.reshape(-1, n_neurons, n_timepoints))
X_test_pt = torch.FloatTensor(X_test.reshape(-1, n_neurons, n_timepoints))
y_train_pt = torch.LongTensor(y_train)
y_test_pt = torch.LongTensor(y_test)

# Create data loaders
from torch.utils.data import TensorDataset, DataLoader

train_dataset = TensorDataset(X_train_pt, y_train_pt)
test_dataset = TensorDataset(X_test_pt, y_test_pt)
train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)
test_loader = DataLoader(test_dataset, batch_size=32, shuffle=False)

# Define CNN decoder architecture
class CNNDecoder(nn.Module):
    def __init__(self, n_neurons, n_timepoints, n_classes):
        super(CNNDecoder, self).__init__()

        # Convolutional layers to process spatiotemporal patterns
        self.conv1 = nn.Conv1d(n_neurons, 64, kernel_size=3, padding=1)
        self.bn1 = nn.BatchNorm1d(64)
        self.conv2 = nn.Conv1d(64, 128, kernel_size=3, padding=1)
        self.bn2 = nn.BatchNorm1d(128)
        self.pool = nn.MaxPool1d(2)

        # Calculate size after convolutions
        conv_output_size = 128 * (n_timepoints // 2)

        # Fully connected layers
        self.fc1 = nn.Linear(conv_output_size, 256)
        self.dropout = nn.Dropout(0.5)
        self.fc2 = nn.Linear(256, n_classes)

    def forward(self, x):
        # x shape: (batch, neurons, time)
        x = torch.relu(self.bn1(self.conv1(x)))
        x = self.pool(torch.relu(self.bn2(self.conv2(x))))

        # Flatten
        x = x.view(x.size(0), -1)

        # Fully connected
        x = torch.relu(self.fc1(x))
        x = self.dropout(x)
        x = self.fc2(x)

        return x

# Initialize model
model = CNNDecoder(n_neurons, n_timepoints, n_stimuli)
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=0.001)

# Train the model
print(); print("Training CNN decoder...")
n_epochs = 50
train_losses = []
train_accs = []
test_accs = []

for epoch in range(n_epochs):
    model.train()
    epoch_loss = 0
    correct = 0
    total = 0

    for batch_X, batch_y in train_loader:
        optimizer.zero_grad()
        outputs = model(batch_X)
        loss = criterion(outputs, batch_y)
        loss.backward()
        optimizer.step()

        epoch_loss += loss.item()
        _, predicted = torch.max(outputs.data, 1)
        total += batch_y.size(0)
        correct += (predicted == batch_y).sum().item()

    train_loss = epoch_loss / len(train_loader)
    train_acc = correct / total
    train_losses.append(train_loss)
    train_accs.append(train_acc)

    # Evaluate on test set
    model.eval()
    correct = 0
    total = 0
    with torch.no_grad():
        for batch_X, batch_y in test_loader:
            outputs = model(batch_X)
            _, predicted = torch.max(outputs.data, 1)
            total += batch_y.size(0)
            correct += (predicted == batch_y).sum().item()

    test_acc = correct / total
    test_accs.append(test_acc)

    if (epoch + 1) % 10 == 0:
        print(f'Epoch [{epoch+1}/{n_epochs}], Loss: {train_loss:.4f}, '
              f'Train Acc: {train_acc:.3f}, Test Acc: {test_acc:.3f}')

print(f"Final CNN Decoder Accuracy: {test_accs[-1]:.3f}")
print(f"Improvement over linear: {(test_accs[-1] - linear_accuracy):.3f} "
      f"({((test_accs[-1] - linear_accuracy)/linear_accuracy * 100):.1f}%)")

24.13.5 Code Block 5: t-SNE Nonlinear Dimensionality Reduction for Neural Manifolds

# Apply t-SNE to time-averaged activity per trial
X_trial_avg = X_traj.mean(axis=2)  # Average over time: (trials, neurons)

print("Running t-SNE (this may take a minute)...")
tsne = TSNE(n_components=2, random_state=42, perplexity=30)
X_tsne = tsne.fit_transform(X_trial_avg)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Panel 1: t-SNE embedding colored by condition
ax = axes[0]
scatter0 = ax.scatter(X_tsne[conditions==0, 0], X_tsne[conditions==0, 1],
                       c='#0066cc', s=50, alpha=0.6, label='Condition 0', edgecolors='black', linewidth=0.5)
scatter1 = ax.scatter(X_tsne[conditions==1, 0], X_tsne[conditions==1, 1],
                       c='#cc0000', s=50, alpha=0.6, label='Condition 1', edgecolors='black', linewidth=0.5)
ax.set_xlabel('t-SNE Dimension 1', fontsize=11)
ax.set_ylabel('t-SNE Dimension 2', fontsize=11)
ax.set_title('t-SNE Embedding: Trial-Averaged Activity', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 2: Compare PCA vs t-SNE separation
ax = axes[1]
# Calculate silhouette score (measure of cluster separation)
from sklearn.metrics import silhouette_score

# PCA separation (using first 2 PCs of trial-averaged data)
X_trial_avg_pca = pca.transform(X_trial_avg)[:, :2]
silhouette_pca = silhouette_score(X_trial_avg_pca, conditions)

# t-SNE separation
silhouette_tsne = silhouette_score(X_tsne, conditions)

methods = ['PCA', 't-SNE']
silhouettes = [silhouette_pca, silhouette_tsne]
colors = ['#0066cc', '#cc0000']
bars = ax.bar(methods, silhouettes, color=colors, alpha=0.7, edgecolor='black', linewidth=2)

for bar, score in zip(bars, silhouettes):
    height = bar.get_height()
    ax.text(bar.get_x() + bar.get_width()/2., height,
            f'{score:.3f}',
            ha='center', va='bottom', fontsize=12, fontweight='bold')

ax.set_ylabel('Silhouette Score - (Condition Separation)', fontsize=11)
ax.set_title('PCA vs t-SNE: Cluster Separation', fontsize=12, fontweight='bold')
ax.set_ylim([0, 1])
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print(f"PCA Silhouette Score: {silhouette_pca:.3f}")
print(f"t-SNE Silhouette Score: {silhouette_tsne:.3f}")
print(f"t-SNE provides {((silhouette_tsne - silhouette_pca)/silhouette_pca * 100):.1f}% better separation")

24.13.6 Code Block 6: UMAP Dimensionality Reduction and Method Comparison

# Note: UMAP requires installation: pip install umap-learn
try:
    import umap

    print("Running UMAP...")
    umap_model = umap.UMAP(n_components=2, random_state=42, n_neighbors=15)
    X_umap = umap_model.fit_transform(X_trial_avg)

    fig, axes = plt.subplots(1, 3, figsize=(18, 5))

    # Panel 1: UMAP embedding
    ax = axes[0]
    ax.scatter(X_umap[conditions==0, 0], X_umap[conditions==0, 1],
               c='#0066cc', s=50, alpha=0.6, label='Condition 0', edgecolors='black', linewidth=0.5)
    ax.scatter(X_umap[conditions==1, 0], X_umap[conditions==1, 1],
               c='#cc0000', s=50, alpha=0.6, label='Condition 1', edgecolors='black', linewidth=0.5)
    ax.set_xlabel('UMAP Dimension 1', fontsize=11)
    ax.set_ylabel('UMAP Dimension 2', fontsize=11)
    ax.set_title('UMAP Embedding', fontsize=12, fontweight='bold')
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)

    # Panel 2: All three methods side by side
    ax = axes[1]
    ax.scatter(X_trial_avg_pca[conditions==0, 0], X_trial_avg_pca[conditions==0, 1],
               c='#0066cc', s=30, alpha=0.5, label='Condition 0', marker='o')
    ax.scatter(X_trial_avg_pca[conditions==1, 0], X_trial_avg_pca[conditions==1, 1],
               c='#cc0000', s=30, alpha=0.5, label='Condition 1', marker='s')
    ax.set_xlabel('Dimension 1', fontsize=11)
    ax.set_ylabel('Dimension 2', fontsize=11)
    ax.set_title('PCA (Linear)', fontsize=12, fontweight='bold')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    # Panel 3: Comparison metrics
    ax = axes[2]
    silhouette_umap = silhouette_score(X_umap, conditions)

    methods = ['PCA', 't-SNE', 'UMAP']
    silhouettes = [silhouette_pca, silhouette_tsne, silhouette_umap]
    colors = ['#0066cc', '#cc0000', '#9966cc']
    bars = ax.bar(methods, silhouettes, color=colors, alpha=0.7, edgecolor='black', linewidth=2)

    for bar, score in zip(bars, silhouettes):
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2., height,
                f'{score:.3f}',
                ha='center', va='bottom', fontsize=11, fontweight='bold')

    ax.set_ylabel('Silhouette Score', fontsize=11)
    ax.set_title('Method Comparison', fontsize=12, fontweight='bold')
    ax.set_ylim([0, 1])
    ax.grid(True, alpha=0.3, axis='y')

    plt.tight_layout()
    plt.show()

    print(f"UMAP Silhouette Score: {silhouette_umap:.3f}")
    print("UMAP often provides best balance of local and global structure")

except ImportError:
    print("UMAP not installed. Install with: pip install umap-learn")
    print("Skipping UMAP analysis...")

---

24.14 Code Examples from Chapter 24 (Part 4: Lines 1500-End)

24.14.1 Code Block 1: CEBRA-Style Contrastive Learning for Neural Embeddings

# Conceptual illustration of CEBRA-style contrastive learning
import numpy as np
import matplotlib.pyplot as plt
from sklearn.manifold import MDS

np.random.seed(42)

# Simulate neural activity during 4 distinct behavioral states
n_neurons = 50
n_samples_per_state = 30
n_states = 4

# Generate neural patterns for each state
neural_data = []
labels = []

for state in range(n_states):
    # Each state has a characteristic population pattern
    base_pattern = np.random.randn(n_neurons) * 2
    for sample in range(n_samples_per_state):
        # Add trial-to-trial variability
        pattern = base_pattern + np.random.randn(n_neurons) * 0.5
        neural_data.append(pattern)
        labels.append(state)

neural_data = np.array(neural_data)
labels = np.array(labels)

# Simulate CEBRA embedding (using MDS as approximation)
mds = MDS(n_components=2, random_state=42)
embeddings = mds.fit_transform(neural_data)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Panel 1: High-dimensional neural space (projected to 2D for visualization)
ax = axes[0]
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
neural_pca = pca.fit_transform(neural_data)

colors = ['#0066cc', '#cc0000', '#00cc66', '#cc9900']
for state in range(n_states):
    mask = labels == state
    ax.scatter(neural_pca[mask, 0], neural_pca[mask, 1],
               c=colors[state], label=f'Behavior {state+1}',
               alpha=0.6, s=50, edgecolors='black', linewidth=0.5)
ax.set_xlabel('PC1', fontsize=11)
ax.set_ylabel('PC2', fontsize=11)
ax.set_title('Raw Neural Activity (PCA)', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 2: CEBRA embedding (better separation)
ax = axes[1]
for state in range(n_states):
    mask = labels == state
    ax.scatter(embeddings[mask, 0], embeddings[mask, 1],
               c=colors[state], label=f'Behavior {state+1}',
               alpha=0.6, s=50, edgecolors='black', linewidth=0.5)
ax.set_xlabel('CEBRA Dimension 1', fontsize=11)
ax.set_ylabel('CEBRA Dimension 2', fontsize=11)
ax.set_title('CEBRA Embedding (Contrastive Learning)', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 3: Quantify separation improvement
from sklearn.metrics import silhouette_score
silhouette_pca = silhouette_score(neural_pca, labels)
silhouette_cebra = silhouette_score(embeddings, labels)

ax = axes[2]
methods = ['PCA', 'CEBRA']
scores = [silhouette_pca, silhouette_cebra]
bars = ax.bar(methods, scores, color=['#0066cc', '#cc0000'],
              alpha=0.7, edgecolor='black', linewidth=2)

for bar, score in zip(bars, scores):
    height = bar.get_height()
    ax.text(bar.get_x() + bar.get_width()/2., height,
            f'{score:.3f}',
            ha='center', va='bottom', fontsize=12, fontweight='bold')

ax.set_ylabel('Behavioral Clustering - (Silhouette Score)', fontsize=11)
ax.set_title('Separation Quality', fontsize=12, fontweight='bold')
ax.set_ylim([0, 1])
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print(f"PCA Silhouette: {silhouette_pca:.3f}")
print(f"CEBRA Silhouette: {silhouette_cebra:.3f}")
print(f"Improvement: {((silhouette_cebra - silhouette_pca) / silhouette_pca * 100):.1f}%")

24.14.2 Code Block 2: Brain-Score Benchmarking

# Visualize Brain-Score concept: comparing model representations to neural data
import numpy as np
import matplotlib.pyplot as plt

# Simulate neural responses and model predictions
n_images = 100
n_neurons = 50

np.random.seed(42)

# True neural responses to images
neural_responses = np.random.randn(n_images, n_neurons)
neural_responses = np.cumsum(neural_responses, axis=0)  # Add structure

# Simulate different model predictions
model_predictions = {
    'Random CNN': neural_responses + np.random.randn(n_images, n_neurons) * 3,
    'ImageNet CNN': neural_responses + np.random.randn(n_images, n_neurons) * 1.5,
    'Self-Supervised': neural_responses + np.random.randn(n_images, n_neurons) * 0.8,
    'Vision Transformer': neural_responses + np.random.randn(n_images, n_neurons) * 0.5,
}

# Calculate correlation for each model
correlations = {}
for model_name, predictions in model_predictions.items():
    corrs = []
    for neuron in range(n_neurons):
        corr = np.corrcoef(neural_responses[:, neuron], predictions[:, neuron])[0, 1]
        corrs.append(corr)
    correlations[model_name] = np.mean(corrs)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Panel 1: Example predictions vs actual
ax = axes[0]
example_neuron = 10
ax.scatter(neural_responses[:, example_neuron],
           model_predictions['Random CNN'][:, example_neuron],
           alpha=0.5, s=30, label='Random CNN', color='#cccccc')
ax.scatter(neural_responses[:, example_neuron],
           model_predictions['Vision Transformer'][:, example_neuron],
           alpha=0.5, s=30, label='Vision Transformer', color='#cc0000')
ax.plot([neural_responses[:, example_neuron].min(), neural_responses[:, example_neuron].max()],
        [neural_responses[:, example_neuron].min(), neural_responses[:, example_neuron].max()],
        'k--', linewidth=2, label='Perfect Match')
ax.set_xlabel('True Neural Response', fontsize=11)
ax.set_ylabel('Model Prediction', fontsize=11)
ax.set_title(f'Prediction vs Reality: Example Neuron #{example_neuron}', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 2: Brain-Score comparison
ax = axes[1]
models = list(correlations.keys())
scores = list(correlations.values())
colors = ['#cccccc', '#0066cc', '#9966cc', '#cc0000']

bars = ax.barh(models, scores, color=colors, alpha=0.7, edgecolor='black', linewidth=1.5)

for bar, score in zip(bars, scores):
    width = bar.get_width()
    ax.text(width, bar.get_y() + bar.get_height()/2.,
            f'{score:.3f}',
            ha='left', va='center', fontsize=11, fontweight='bold', color='black')

ax.set_xlabel('Neural Predictivity (Correlation)', fontsize=11)
ax.set_title('Brain-Score: Model Performance', fontsize=12, fontweight='bold')
ax.set_xlim([0, 1])
ax.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()

print(); print("Brain-Score Rankings:")
for i, (model, score) in enumerate(sorted(correlations.items(), key=lambda x: x[1], reverse=True), 1):
    print(f"{i}. {model}: {score:.3f}")

24.14.3 Code Block 3: Transfer Learning Efficiency

# Illustrate transfer learning benefit
import numpy as np
import matplotlib.pyplot as plt

# Simulate data requirements and performance
n_samples = np.array([10, 20, 50, 100, 200, 500, 1000, 2000])

# Performance with training from scratch
perf_scratch = 0.5 + 0.4 * (1 - np.exp(-n_samples / 500))

# Performance with transfer learning (pre-trained model)
perf_transfer = 0.65 + 0.3 * (1 - np.exp(-n_samples / 100))

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Panel 1: Learning curves
ax = axes[0]
ax.plot(n_samples, perf_scratch, 'o-', linewidth=2, markersize=8,
        color='#0066cc', label='Train from Scratch')
ax.plot(n_samples, perf_transfer, 's-', linewidth=2, markersize=8,
        color='#cc0000', label='Transfer Learning')
ax.axhline(y=0.9, color='gray', linestyle='--', linewidth=1.5, alpha=0.7, label='Target Performance')

# Highlight data efficiency
target_perf = 0.85
scratch_needed = n_samples[np.argmin(np.abs(perf_scratch - target_perf))]
transfer_needed = n_samples[np.argmin(np.abs(perf_transfer - target_perf))]

ax.axvline(x=scratch_needed, color='#0066cc', linestyle=':', alpha=0.5)
ax.axvline(x=transfer_needed, color='#cc0000', linestyle=':', alpha=0.5)

ax.set_xlabel('Training Samples', fontsize=11)
ax.set_ylabel('Decoding Performance', fontsize=11)
ax.set_title('Transfer Learning Efficiency', fontsize=12, fontweight='bold')
ax.set_xscale('log')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 2: Data efficiency gain
ax = axes[1]
data_reduction = []
for target in np.linspace(0.6, 0.9, 10):
    if np.max(perf_scratch) >= target and np.max(perf_transfer) >= target:
        scratch_n = n_samples[np.argmin(np.abs(perf_scratch - target))]
        transfer_n = n_samples[np.argmin(np.abs(perf_transfer - target))]
        reduction = scratch_n / transfer_n
        data_reduction.append((target, reduction))

targets, reductions = zip(*data_reduction)
ax.plot(targets, reductions, 'o-', linewidth=2, markersize=8, color='#9966cc')
ax.set_xlabel('Target Performance', fontsize=11)
ax.set_ylabel('Data Efficiency - (Fold Reduction)', fontsize=11)
ax.set_title('Transfer Learning Advantage', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)

# Annotate typical operating point
ax.annotate(f'{reductions[5]:.1f}x fewer - samples needed',
            xy=(targets[5], reductions[5]),
            xytext=(targets[5] - 0.05, reductions[5] + 1),
            arrowprops=dict(arrowstyle='->', color='black', lw=1.5),
            fontsize=10, fontweight='bold',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

plt.tight_layout()
plt.show()

print(f" - To reach 85% performance:")
print(f"  Training from scratch: {scratch_needed} samples")
print(f"  Transfer learning: {transfer_needed} samples")
print(f"  Efficiency gain: {scratch_needed/transfer_needed:.1f}x")

24.14.4 Code Block 4: Active Learning for Receptive Field Mapping

# Simulate active learning for receptive field mapping
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# True receptive field (unknown to experimenter)
true_rf_center = (30, 40)
true_rf_sigma = 15

def neuron_response(x, y):
    """Simulate neuron with Gaussian receptive field."""
    dist = np.sqrt((x - true_rf_center[0])**2 + (y - true_rf_center[1])**2)
    response = np.exp(-(dist**2) / (2 * true_rf_sigma**2))
    # Add noise
    response += np.random.randn() * 0.1
    return np.clip(response, 0, 1)

# Active learning: choose next stimulus location based on uncertainty
def random_sampling(n_samples, grid_size=80):
    """Baseline: random stimulus locations."""
    samples = []
    for _ in range(n_samples):
        x = np.random.randint(0, grid_size)
        y = np.random.randint(0, grid_size)
        response = neuron_response(x, y)
        samples.append((x, y, response))
    return samples

def active_sampling(n_samples, grid_size=80):
    """Active learning: sample where uncertainty is highest."""
    samples = []
    # Start with a few random samples
    for _ in range(5):
        x = np.random.randint(0, grid_size)
        y = np.random.randint(0, grid_size)
        response = neuron_response(x, y)
        samples.append((x, y, response))

    # Then sample near high-response regions (exploitation)
    # and unexplored regions (exploration)
    for _ in range(n_samples - 5):
        if np.random.rand() < 0.7:  # Exploitation
            # Sample near previously strong responses
            strong_samples = [s for s in samples if s[2] > 0.5]
            if strong_samples:
                base_x, base_y, _ = strong_samples[np.random.randint(len(strong_samples))]
                x = int(np.clip(base_x + np.random.randn() * 10, 0, grid_size-1))
                y = int(np.clip(base_y + np.random.randn() * 10, 0, grid_size-1))
            else:
                x = np.random.randint(0, grid_size)
                y = np.random.randint(0, grid_size)
        else:  # Exploration
            x = np.random.randint(0, grid_size)
            y = np.random.randint(0, grid_size)

        response = neuron_response(x, y)
        samples.append((x, y, response))

    return samples

# Compare strategies
n_samples_list = [10, 20, 50, 100]
grid_size = 80

fig, axes = plt.subplots(2, 4, figsize=(18, 9))

for idx, n_samples in enumerate(n_samples_list):
    # Random sampling
    random_samples = random_sampling(n_samples, grid_size)
    x_rand, y_rand, r_rand = zip(*random_samples)

    # Active sampling
    active_samples = active_sampling(n_samples, grid_size)
    x_active, y_active, r_active = zip(*active_samples)

    # Plot random sampling
    ax = axes[0, idx]
    scatter = ax.scatter(x_rand, y_rand, c=r_rand, s=50, cmap='hot',
                         vmin=0, vmax=1, edgecolors='black', linewidth=0.5)
    # Show true RF center
    circle = plt.Circle(true_rf_center, true_rf_sigma, color='cyan',
                        fill=False, linewidth=2, linestyle='--', label='True RF')
    ax.add_patch(circle)
    ax.set_xlim([0, grid_size])
    ax.set_ylim([0, grid_size])
    ax.set_aspect('equal')
    ax.set_title(f'Random Sampling (N={n_samples})', fontsize=11, fontweight='bold')
    if idx == 0:
        ax.set_ylabel('Y Position', fontsize=10)

    # Plot active sampling
    ax = axes[1, idx]
    scatter = ax.scatter(x_active, y_active, c=r_active, s=50, cmap='hot',
                         vmin=0, vmax=1, edgecolors='black', linewidth=0.5)
    circle = plt.Circle(true_rf_center, true_rf_sigma, color='cyan',
                        fill=False, linewidth=2, linestyle='--', label='True RF')
    ax.add_patch(circle)
    ax.set_xlim([0, grid_size])
    ax.set_ylim([0, grid_size])
    ax.set_aspect('equal')
    ax.set_title(f'Active Learning (N={n_samples})', fontsize=11, fontweight='bold')
    ax.set_xlabel('X Position', fontsize=10)
    if idx == 0:
        ax.set_ylabel('Y Position', fontsize=10)

# Add colorbar
cbar = plt.colorbar(scatter, ax=axes, orientation='vertical', fraction=0.02, pad=0.02)
cbar.set_label('Neural Response', fontsize=11)

plt.tight_layout()
plt.show()

print("Active learning concentrates samples in informative regions")
print("Result: More accurate receptive field map with fewer samples")

24.14.5 Code Block 5: Bayesian Optimization for Parameter Tuning

# Simulate Bayesian optimization for parameter tuning
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

np.random.seed(42)

# True (unknown) objective function
def true_objective(x):
    """Simulate response to a parameter (e.g., stimulation intensity)."""
    return np.sin(x * 3) * np.exp(-x/10) + 0.5

# Bayesian optimization
class SimpleBayesianOptimization:
    def __init__(self, bounds, n_initial=3):
        self.bounds = bounds
        self.X_sampled = []
        self.y_sampled = []

        # Initial random samples
        for _ in range(n_initial):
            x = np.random.uniform(bounds[0], bounds[1])
            y = true_objective(x) + np.random.randn() * 0.05
            self.X_sampled.append(x)
            self.y_sampled.append(y)

    def acquisition_function(self, X):
        """Expected Improvement acquisition function (simplified)."""
        if len(self.y_sampled) == 0:
            return np.ones_like(X)

        # Simple heuristic: balance exploration (distance from samples) and exploitation (predicted value)
        exploration = np.ones_like(X)
        for x_samp in self.X_sampled:
            exploration *= np.abs(X - x_samp)

        # Exploitation: interpolate between known points
        exploitation = np.zeros_like(X)
        for i, x in enumerate(X):
            weights = np.exp(-np.array([(x - xs)**2 / 2 for xs in self.X_sampled]))
            weights /= weights.sum() if weights.sum() > 0 else 1
            exploitation[i] = np.dot(weights, self.y_sampled)

        # Combine
        acquisition = 0.5 * exploration / exploration.max() + 0.5 * exploitation
        return acquisition

    def next_sample(self):
        """Choose next point to sample."""
        X_candidate = np.linspace(self.bounds[0], self.bounds[1], 1000)
        acquisition_values = self.acquisition_function(X_candidate)
        next_x = X_candidate[np.argmax(acquisition_values)]

        # Sample the objective
        next_y = true_objective(next_x) + np.random.randn() * 0.05
        self.X_sampled.append(next_x)
        self.y_sampled.append(next_y)

        return next_x, next_y

# Compare Bayesian optimization vs random search
bounds = (0, 10)
n_iterations = 15

# Bayesian optimization
bayes_opt = SimpleBayesianOptimization(bounds, n_initial=3)
bayes_history = [(x, y) for x, y in zip(bayes_opt.X_sampled, bayes_opt.y_sampled)]

for _ in range(n_iterations - 3):
    x, y = bayes_opt.next_sample()
    bayes_history.append((x, y))

# Random search
random_history = []
for _ in range(n_iterations):
    x = np.random.uniform(bounds[0], bounds[1])
    y = true_objective(x) + np.random.randn() * 0.05
    random_history.append((x, y))

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Panel 1: Bayesian optimization progress
ax = axes[0, 0]
X_true = np.linspace(bounds[0], bounds[1], 200)
y_true = true_objective(X_true)
ax.plot(X_true, y_true, 'k-', linewidth=2, label='True Function', alpha=0.7)

x_bayes, y_bayes = zip(*bayes_history)
ax.scatter(x_bayes, y_bayes, c=range(len(x_bayes)), cmap='cool',
           s=100, edgecolors='black', linewidth=1.5, label='Sampled Points', zorder=5)

# Show best found
best_idx = np.argmax(y_bayes)
ax.scatter([x_bayes[best_idx]], [y_bayes[best_idx]], c='red', s=300,
           marker='*', edgecolors='black', linewidth=2, label='Best Found', zorder=10)

ax.set_xlabel('Parameter Value', fontsize=11)
ax.set_ylabel('Objective (Response)', fontsize=11)
ax.set_title('Bayesian Optimization Trajectory', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Panel 2: Random search progress
ax = axes[0, 1]
ax.plot(X_true, y_true, 'k-', linewidth=2, label='True Function', alpha=0.7)

x_random, y_random = zip(*random_history)
ax.scatter(x_random, y_random, c=range(len(x_random)), cmap='cool',
           s=100, edgecolors='black', linewidth=1.5, label='Sampled Points', zorder=5)

best_idx_random = np.argmax(y_random)
ax.scatter([x_random[best_idx_random]], [y_random[best_idx_random]], c='red', s=300,
           marker='*', edgecolors='black', linewidth=2, label='Best Found', zorder=10)

ax.set_xlabel('Parameter Value', fontsize=11)
ax.set_ylabel('Objective (Response)', fontsize=11)
ax.set_title('Random Search Trajectory', fontsize=12, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Panel 3: Convergence comparison
ax = axes[1, 0]
bayes_best = [max(y_bayes[:i+1]) for i in range(len(y_bayes))]
random_best = [max(y_random[:i+1]) for i in range(len(y_random))]
true_optimum = np.max(y_true)

ax.plot(bayes_best, 'o-', linewidth=2, markersize=8, color='#0066cc',
        label='Bayesian Optimization')
ax.plot(random_best, 's-', linewidth=2, markersize=8, color='#cc0000',
        label='Random Search')
ax.axhline(y=true_optimum, color='green', linestyle='--', linewidth=2,
           label='True Optimum', alpha=0.7)

ax.set_xlabel('Iteration', fontsize=11)
ax.set_ylabel('Best Value Found', fontsize=11)
ax.set_title('Convergence Comparison', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# Panel 4: Efficiency gain
ax = axes[1, 1]
# Calculate how many iterations needed to reach 90% of optimum
target = 0.9 * true_optimum

bayes_to_target = next((i for i, val in enumerate(bayes_best) if val >= target), len(bayes_best))
random_to_target = next((i for i, val in enumerate(random_best) if val >= target), len(random_best))

categories = ['Bayesian - Optimization', 'Random - Search']
iterations = [bayes_to_target + 1, random_to_target + 1]
colors = ['#0066cc', '#cc0000']

bars = ax.bar(categories, iterations, color=colors, alpha=0.7,
              edgecolor='black', linewidth=2)

for bar, iters in zip(bars, iterations):
    height = bar.get_height()
    ax.text(bar.get_x() + bar.get_width()/2., height,
            f'{iters}',
            ha='center', va='bottom', fontsize=12, fontweight='bold')

ax.set_ylabel('Iterations to 90% Optimum', fontsize=11)
ax.set_title('Efficiency: Iterations Required', fontsize=12, fontweight='bold')
ax.set_ylim([0, max(iterations) * 1.2])
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print(f"Bayesian Optimization reached target in {bayes_to_target + 1} iterations")
print(f"Random Search reached target in {random_to_target + 1} iterations")
print(f"Efficiency gain: {random_to_target / bayes_to_target:.1f}x")