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."""
    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")