import numpy as np

class MatrixFactorization:
    """
    Collaborative filtering using gradient descent.
    This is similar to what Netflix uses!
    """
    
    def __init__(self, n_users, n_items, n_factors=20, lr=0.01, reg=0.02):
        self.n_factors = n_factors
        self.lr = lr
        self.reg = reg
        
        # Initialize parameters
        self.user_factors = np.random.randn(n_users, n_factors) * 0.1
        self.item_factors = np.random.randn(n_items, n_factors) * 0.1
        self.user_bias = np.zeros(n_users)
        self.item_bias = np.zeros(n_items)
        self.global_bias = 0
    
    def predict(self, user, item):
        """Predict rating for user-item pair."""
        return (self.global_bias + 
                self.user_bias[user] + 
                self.item_bias[item] + 
                np.dot(self.user_factors[user], self.item_factors[item]))
    
    def train(self, ratings, n_epochs=50):
        """
        Train using stochastic gradient descent.
        
        ratings: list of (user, item, rating) tuples
        """
        self.global_bias = np.mean([r for _, _, r in ratings])
        
        for epoch in range(n_epochs):
            np.random.shuffle(ratings)
            total_loss = 0
            
            for user, item, rating in ratings:
                # Prediction
                pred = self.predict(user, item)
                error = rating - pred
                
                # Gradients (derivatives of loss w.r.t. parameters)
                # ∂L/∂b_u = -2*error + 2*λ*b_u
                # ∂L/∂b_i = -2*error + 2*λ*b_i
                # ∂L/∂p_u = -2*error*q_i + 2*λ*p_u
                # ∂L/∂q_i = -2*error*p_u + 2*λ*q_i
                
                # Update biases
                self.user_bias[user] += self.lr * (error - self.reg * self.user_bias[user])
                self.item_bias[item] += self.lr * (error - self.reg * self.item_bias[item])
                
                # Update latent factors
                user_factor_old = self.user_factors[user].copy()
                self.user_factors[user] += self.lr * (error * self.item_factors[item] - 
                                                       self.reg * self.user_factors[user])
                self.item_factors[item] += self.lr * (error * user_factor_old - 
                                                       self.reg * self.item_factors[item])
                
                total_loss += error**2
            
            rmse = np.sqrt(total_loss / len(ratings))
            if epoch % 10 == 0:
                print(f"Epoch {epoch}: RMSE = {rmse:.4f}")
        
        return self

# Demo with synthetic movie ratings
np.random.seed(42)
n_users, n_items = 100, 50

# Generate ratings
ratings = []
for u in range(n_users):
    for i in range(n_items):
        if np.random.random() < 0.1:  # 10% density
            # True rating based on hidden factors
            true_rating = 3 + np.random.randn() * 0.5
            ratings.append((u, i, np.clip(true_rating, 1, 5)))

print(f"Training on {len(ratings)} ratings...")
model = MatrixFactorization(n_users, n_items, n_factors=10)
model.train(ratings, n_epochs=50)

# Make predictions
print("\nSample predictions:")
for u in [0, 10, 50]:
    preds = [(i, model.predict(u, i)) for i in range(n_items)]
    preds.sort(key=lambda x: x[1], reverse=True)
    print(f"User {u} top recommendations: items {[p[0] for p in preds[:5]]}")

import numpy as np

def softmax(x, axis=-1):
    """
    Numerically stable softmax.
    
    The key trick: subtract max(x) before exponentiating.
    
    Why this matters: if x contains a value like 1000, then exp(1000)
    overflows to infinity. But exp(1000 - 1000) = exp(0) = 1, which is fine.
    Mathematically, subtracting a constant from all inputs does not change
    the softmax output (the constant cancels in numerator and denominator),
    but it prevents overflow.
    
    This is one of the most important numerical stability tricks in ML.
    Every production softmax implementation does this.
    """
    x_max = np.max(x, axis=axis, keepdims=True)
    exp_x = np.exp(x - x_max)
    return exp_x / np.sum(exp_x, axis=axis, keepdims=True)

def softmax_gradient(softmax_output):
    """
    Compute the Jacobian of softmax.
    
    For each output s_i, we need ∂s_i/∂z_j for all j.
    Jacobian[i,j] = s_i * (δ_ij - s_j)
    """
    s = softmax_output.reshape(-1, 1)
    return np.diagflat(s) - np.dot(s, s.T)

def scaled_dot_product_attention(Q, K, V):
    """
    Single-head attention mechanism.
    
    Q: (seq_len, d_k) query vectors
    K: (seq_len, d_k) key vectors
    V: (seq_len, d_v) value vectors
    """
    d_k = Q.shape[-1]
    
    # Attention scores
    scores = Q @ K.T / np.sqrt(d_k)
    
    # Softmax to get attention weights
    attention_weights = softmax(scores, axis=-1)
    
    # Weighted sum of values
    output = attention_weights @ V
    
    return output, attention_weights

# Demo
np.random.seed(42)
seq_len, d_k, d_v = 5, 8, 8

Q = np.random.randn(seq_len, d_k)
K = np.random.randn(seq_len, d_k)
V = np.random.randn(seq_len, d_v)

output, weights = scaled_dot_product_attention(Q, K, V)

print("Attention weights (each row sums to 1):")
print(weights.round(3))
print(f"\nRow sums: {weights.sum(axis=1).round(3)}")

# Visualize attention
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 6))
plt.imshow(weights, cmap='Blues')
plt.colorbar(label='Attention Weight')
plt.xlabel('Key Position')
plt.ylabel('Query Position')
plt.title('Attention Pattern (Self-Attention)')
for i in range(seq_len):
    for j in range(seq_len):
        plt.text(j, i, f'{weights[i,j]:.2f}', ha='center', va='center')
plt.show()

import numpy as np
from scipy import signal

def conv2d_forward(x, kernel, stride=1, padding=0):
    """Forward pass of 2D convolution."""
    if padding > 0:
        x = np.pad(x, padding, mode='constant')
    
    h, w = x.shape
    kh, kw = kernel.shape
    out_h = (h - kh) // stride + 1
    out_w = (w - kw) // stride + 1
    
    output = np.zeros((out_h, out_w))
    
    for i in range(out_h):
        for j in range(out_w):
            region = x[i*stride:i*stride+kh, j*stride:j*stride+kw]
            output[i, j] = np.sum(region * kernel)
    
    return output

def conv2d_backward(grad_output, x, kernel, stride=1, padding=0):
    """
    Backward pass: compute gradients w.r.t. input and kernel.
    
    This is where calculus meets convolution!
    """
    if padding > 0:
        x = np.pad(x, padding, mode='constant')
    
    kh, kw = kernel.shape
    
    # Gradient w.r.t. kernel
    # Each weight sees a portion of the input
    grad_kernel = np.zeros_like(kernel)
    for i in range(kh):
        for j in range(kw):
            # Sum over all positions where this weight was used
            grad_kernel[i, j] = np.sum(
                x[i:i+grad_output.shape[0]*stride:stride, 
                  j:j+grad_output.shape[1]*stride:stride] * grad_output
            )
    
    # Gradient w.r.t. input (full convolution)
    rotated_kernel = np.rot90(kernel, 2)
    grad_input = signal.convolve2d(grad_output, rotated_kernel, mode='full')
    
    return grad_input, grad_kernel

# Demo
np.random.seed(42)

# Simple image
image = np.random.randn(8, 8)

# Edge detection kernel
kernel = np.array([[-1, 0, 1],
                   [-2, 0, 2],
                   [-1, 0, 1]])

# Forward pass
output = conv2d_forward(image, kernel, padding=1)

# Backward pass (assuming unit gradient from next layer)
grad_output = np.ones_like(output)
grad_input, grad_kernel = conv2d_backward(grad_output, image, kernel, padding=1)

# Visualize
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].imshow(image, cmap='gray')
axes[0, 0].set_title('Input Image')
axes[0, 0].axis('off')

axes[0, 1].imshow(kernel, cmap='RdBu', vmin=-2, vmax=2)
axes[0, 1].set_title('Kernel (Sobel Edge)')
axes[0, 1].axis('off')

axes[0, 2].imshow(output, cmap='gray')
axes[0, 2].set_title('Output (Edges Detected)')
axes[0, 2].axis('off')

axes[1, 0].imshow(grad_input[1:-1, 1:-1], cmap='RdBu')
axes[1, 0].set_title('Gradient w.r.t. Input')
axes[1, 0].axis('off')

axes[1, 1].imshow(grad_kernel, cmap='RdBu')
axes[1, 1].set_title('Gradient w.r.t. Kernel')
axes[1, 1].axis('off')

# Show gradient flow
axes[1, 2].text(0.5, 0.5, 
    "Gradient flows from\nloss → output → input\n\n"
    "∂L/∂kernel = ∑ input * ∂L/∂output\n"
    "∂L/∂input = conv(∂L/∂output, rot180(kernel))",
    ha='center', va='center', fontsize=12,
    transform=axes[1, 2].transAxes)
axes[1, 2].axis('off')

plt.tight_layout()
plt.show()

import numpy as np

class PolicyGradientAgent:
    """
    Simple REINFORCE agent for cart-pole style problems.
    
    This is how AI learns to play games!
    """
    
    def __init__(self, state_dim, action_dim, hidden_dim=32, lr=0.01):
        self.lr = lr
        
        # Policy network weights
        self.W1 = np.random.randn(state_dim, hidden_dim) * 0.1
        self.b1 = np.zeros(hidden_dim)
        self.W2 = np.random.randn(hidden_dim, action_dim) * 0.1
        self.b2 = np.zeros(action_dim)
        
        # For storing episode data
        self.saved_log_probs = []
        self.rewards = []
    
    def forward(self, state):
        """Forward pass through policy network."""
        self.h = np.maximum(0, state @ self.W1 + self.b1)  # ReLU
        logits = self.h @ self.W2 + self.b2
        probs = self._softmax(logits)
        return probs
    
    def _softmax(self, x):
        exp_x = np.exp(x - np.max(x))
        return exp_x / exp_x.sum()
    
    def select_action(self, state):
        """Sample action from policy."""
        probs = self.forward(state)
        action = np.random.choice(len(probs), p=probs)
        
        # Store log probability for gradient computation
        self.saved_log_probs.append(np.log(probs[action]))
        self.last_state = state
        self.last_action = action
        self.last_probs = probs
        
        return action
    
    def update(self):
        """
        Update policy using REINFORCE.
        
        The key equation:
        ∇_θ J(θ) = Σ ∇_θ log π(a|s) * R
        """
        # Compute discounted returns
        gamma = 0.99
        R = 0
        returns = []
        for r in self.rewards[::-1]:
            R = r + gamma * R
            returns.insert(0, R)
        
        returns = np.array(returns)
        returns = (returns - returns.mean()) / (returns.std() + 1e-8)  # Normalize
        
        # Compute policy gradient and update
        # This is gradient descent on the negative expected reward
        for log_prob, R in zip(self.saved_log_probs, returns):
            # In practice, we'd compute full gradients here
            # For demo, we just show the gradient direction
            gradient_direction = log_prob * R
            
        # Clear episode data
        self.saved_log_probs = []
        self.rewards = []
        
        return np.mean(returns)

# Simulate training
agent = PolicyGradientAgent(state_dim=4, action_dim=2)

print("Policy Gradient Training (Simulated)")
print("=" * 50)

for episode in range(10):
    state = np.random.randn(4)  # Random state
    total_reward = 0
    
    for t in range(100):
        action = agent.select_action(state)
        reward = 1.0  # Simplified reward
        agent.rewards.append(reward)
        total_reward += reward
        state = np.random.randn(4)  # Next state
    
    avg_return = agent.update()
    print(f"Episode {episode}: Total Reward = {total_reward}, Avg Return = {avg_return:.4f}")

import numpy as np
import matplotlib.pyplot as plt

class ExtendedKalmanFilter:
    """
    EKF for vehicle tracking.
    State: [x, y, velocity, heading]
    """
    
    def __init__(self):
        # State: [x, y, v, θ]
        self.x = np.zeros(4)
        
        # Covariance matrix
        self.P = np.eye(4) * 1.0
        
        # Process noise
        self.Q = np.diag([0.1, 0.1, 0.01, 0.01])
        
        # Measurement noise (for position only)
        self.R = np.diag([0.5, 0.5])
    
    def predict(self, dt):
        """
        Predict step using nonlinear motion model.
        
        x_new = x + v*cos(θ)*dt
        y_new = y + v*sin(θ)*dt
        v_new = v
        θ_new = θ
        """
        x, y, v, theta = self.x
        
        # Predicted state
        x_pred = np.array([
            x + v * np.cos(theta) * dt,
            y + v * np.sin(theta) * dt,
            v,
            theta
        ])
        
        # Jacobian of motion model (THE CALCULUS!)
        # ∂f/∂x
        F = np.array([
            [1, 0, np.cos(theta)*dt, -v*np.sin(theta)*dt],
            [0, 1, np.sin(theta)*dt,  v*np.cos(theta)*dt],
            [0, 0, 1, 0],
            [0, 0, 0, 1]
        ])
        
        # Predict covariance
        P_pred = F @ self.P @ F.T + self.Q
        
        self.x = x_pred
        self.P = P_pred
        
        return x_pred
    
    def update(self, z):
        """
        Update step using position measurement.
        z = [x_measured, y_measured]
        """
        # Measurement model: H extracts position from state
        H = np.array([
            [1, 0, 0, 0],
            [0, 1, 0, 0]
        ])
        
        # Innovation
        y = z - H @ self.x
        
        # Innovation covariance
        S = H @ self.P @ H.T + self.R
        
        # Kalman gain
        K = self.P @ H.T @ np.linalg.inv(S)
        
        # Update state and covariance
        self.x = self.x + K @ y
        self.P = (np.eye(4) - K @ H) @ self.P
        
        return self.x

# Simulate vehicle tracking
np.random.seed(42)

ekf = ExtendedKalmanFilter()
ekf.x = np.array([0, 0, 5.0, np.pi/6])  # Start: origin, 5 m/s, 30° heading

true_positions = []
measured_positions = []
estimated_positions = []

dt = 0.1
n_steps = 100

true_x, true_y = 0, 0
true_v, true_theta = 5.0, np.pi/6

for t in range(n_steps):
    # True motion
    true_x += true_v * np.cos(true_theta) * dt
    true_y += true_v * np.sin(true_theta) * dt
    true_positions.append([true_x, true_y])
    
    # Noisy measurement
    z = np.array([true_x, true_y]) + np.random.randn(2) * 0.5
    measured_positions.append(z)
    
    # EKF predict and update
    ekf.predict(dt)
    est = ekf.update(z)
    estimated_positions.append(est[:2])

true_positions = np.array(true_positions)
measured_positions = np.array(measured_positions)
estimated_positions = np.array(estimated_positions)

# Visualize
plt.figure(figsize=(12, 6))
plt.plot(true_positions[:, 0], true_positions[:, 1], 'g-', linewidth=2, label='True Path')
plt.scatter(measured_positions[:, 0], measured_positions[:, 1], c='red', s=10, alpha=0.5, label='Noisy Measurements')
plt.plot(estimated_positions[:, 0], estimated_positions[:, 1], 'b-', linewidth=2, label='EKF Estimate')

plt.xlabel('X Position (m)')
plt.ylabel('Y Position (m)')
plt.title('Extended Kalman Filter: Vehicle Tracking')
plt.legend()
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.show()

# Error analysis
errors = np.sqrt(np.sum((estimated_positions - true_positions)**2, axis=1))
print(f"Average tracking error: {errors.mean():.3f} m")
print(f"Final position error: {errors[-1]:.3f} m")