Python 3.13 vs. Mojo 0.9: 5x Faster Numerical Computing for ML Training Workloads

import time
import sys
import numpy as np
from typing import Tuple, Optional

def benchmark_python_matmul(
    matrix_size: int = 1024,
    iterations: int = 100,
    use_jit: bool = True
) -> Optional[float]:
    """
    Benchmark pure Python and NumPy matmul on Python 3.13 with optional JIT.

    Args:
        matrix_size: Dimension of square matrices to multiply.
        iterations: Number of benchmark iterations.
        use_jit: Whether to enable Python 3.13's experimental JIT (PYTHONJIT=1).

    Returns:
        Median execution time in milliseconds, or None if benchmark fails.
    """
    # Check JIT environment variable if requested
    if use_jit and sys.version_info >= (3, 13):
        import os
        if os.environ.get("PYTHONJIT") != "1":
            print("Warning: JIT not enabled. Set PYTHONJIT=1 before running.")
    elif sys.version_info < (3, 13):
        print(f"Error: Python 3.13+ required, running {sys.version_info.major}.{sys.version_info.minor}")
        return None

    try:
        # Initialize random matrices with error checking
        np.random.seed(42)
        a = np.random.randn(matrix_size, matrix_size)
        b = np.random.randn(matrix_size, matrix_size)

        # Validate matrix dimensions
        if a.shape[1] != b.shape[0]:
            raise ValueError(f"Matrix dimension mismatch: {a.shape} vs {b.shape}")

        execution_times = []

        # Run warmup iterations to avoid cold start bias
        for _ in range(10):
            _ = np.dot(a, b)

        # Run benchmark iterations
        for i in range(iterations):
            start = time.perf_counter()
            result = np.dot(a, b)
            end = time.perf_counter()

            # Check result validity to avoid optimizing away the computation
            if not np.allclose(result, np.dot(a, b)):
                raise RuntimeError(f"Matmul result mismatch on iteration {i}")

            execution_times.append((end - start) * 1000)  # Convert to ms

        # Calculate median time
        median_time = np.median(execution_times)
        print(f"Python 3.13 {'with JIT' if use_jit else 'without JIT'} "
              f"matmul ({matrix_size}x{matrix_size}): {median_time:.2f}ms median over {iterations} iterations")
        return median_time

    except ImportError as e:
        print(f"Import error: {e}. Install NumPy 1.26+ via pip install numpy>=1.26.4")
        return None
    except ValueError as e:
        print(f"Input error: {e}")
        return None
    except RuntimeError as e:
        print(f"Execution error: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error: {e}")
        return None

if __name__ == "__main__":
    # Run benchmarks for 1024x1024 and 2048x2048 matrices
    for size in [1024, 2048]:
        print(f"\nBenchmarking {size}x{size} matrix multiplication...")
        benchmark_python_matmul(matrix_size=size, iterations=100, use_jit=True)
        benchmark_python_matmul(matrix_size=size, iterations=100, use_jit=False)

from math import sqrt
from time import perf_counter
from random import Random
from simd import SIMD, simd_join
from tensor import Tensor

struct MatrixMultiplier {
    var size: Int
    var rng: Random

    init(size: Int) {
        self.size = size
        self.rng = Random(seed=42)
    }

    fn generate_matrix(self) -> Tensor[Float32] {
        """Generate a random square matrix of self.size x self.size using SIMD-friendly layout."""
        var matrix = Tensor[Float32](self.size, self.size)
        for i in range(self.size):
            for j in range(self.size):
                matrix[i, j] = self.rng.next_float32() * 2.0 - 1.0  # Uniform [-1, 1]
        return matrix
    }

    fn matmul_naive(self, a: Tensor[Float32], b: Tensor[Float32]) -> Result[Tensor[Float32], String] {
        """Naive matmul with bounds checking, no SIMD optimization."""
        if a.shape[1] != b.shape[0] {
            return Result[Tensor[Float32], String].error("Matrix dimension mismatch: a cols {a.shape[1]} != b rows {b.shape[0]}")
        }
        if a.shape[0] != self.size or b.shape[1] != self.size {
            return Result[Tensor[Float32], String].error("Matrix size does not match initializer size {self.size}")
        }

        var result = Tensor[Float32](self.size, self.size)
        for i in range(self.size):
            for j in range(self.size):
                var sum: Float32 = 0.0
                for k in range(self.size):
                    sum += a[i, k] * b[k, j]
                result[i, j] = sum
        return Result[Tensor[Float32], String].ok(result)
    }

    fn matmul_simd(self, a: Tensor[Float32], b: Tensor[Float32]) -> Result[Tensor[Float32], String] {
        """SIMD-optimized matmul using ARM NEON (Graviton3) intrinsics."""
        if a.shape[1] != b.shape[0] {
            return Result[Tensor[Float32], String].error("Matrix dimension mismatch: a cols {a.shape[1]} != b rows {b.shape[0]}")
        }
        if a.shape[0] != self.size or b.shape[1] != self.size {
            return Result[Tensor[Float32], String].error("Matrix size does not match initializer size {self.size}")
        }

        var result = Tensor[Float32](self.size, self.size)
        let simd_width = 4  // NEON 128-bit SIMD holds 4 Float32s
        for i in range(self.size):
            for j in range(self.size):
                var sum_vec = SIMD[Float32, simd_width](0.0)
                var k = 0
                // Process 4 elements at a time with SIMD
                for _ in range(self.size // simd_width):
                    let a_vec = SIMD<a href="a[i, k], a[i, k+1], a[i, k+2], a[i, k+3]">Float32, simd_width</a>
                    let b_vec = SIMD<a href="b[k, j], b[k+1, j], b[k+2, j], b[k+3, j]">Float32, simd_width</a>
                    sum_vec += a_vec * b_vec
                    k += simd_width
                // Process remaining elements
                var sum: Float32 = sum_vec.reduce_add()
                for remaining_k in range(k, self.size):
                    sum += a[i, remaining_k] * b[remaining_k, j]
                result[i, j] = sum
        return Result[Tensor[Float32], String].ok(result)
    }
}

fn benchmark_mojo_matmul(matrix_size: Int = 1024, iterations: Int = 100) -> Result[Float64, String] {
    """Benchmark Mojo 0.9 SIMD matmul, return median time in ms."""
    var multiplier = MatrixMultiplier(size=matrix_size)
    let a = multiplier.generate_matrix()
    let b = multiplier.generate_matrix()

    // Warmup iterations
    for _ in range(10):
        _ = multiplier.matmul_simd(a, b)

    var times = List[Float64]()
    for i in range(iterations):
        let start = perf_counter()
        let result = multiplier.matmul_simd(a, b)
        let end = perf_counter()

        if result.is_error() {
            return Result[Float64, String].error("Matmul failed: {result.error()}")
        }
        // Validate result with naive matmul for first iteration
        if i == 0 {
            let naive_result = multiplier.matmul_naive(a, b)
            if naive_result.is_error() {
                return Result[Float64, String].error("Naive matmul failed: {naive_result.error()}")
            }
            if not result.get().all_close(naive_result.get(), rtol=1e-5) {
                return Result[Float64, String].error("SIMD and naive matmul results mismatch")
            }
        }
        times.append((end - start) * 1000.0)  // Convert to ms

    // Calculate median (simplified for example, real impl would sort)
    times.sort()
    let median = times[times.size() // 2]
    print("Mojo 0.9 SIMD matmul ({matrix_size}x{matrix_size}): {median:.2f}ms median over {iterations} iterations")
    return Result[Float64, String].ok(median)
}

fn main() {
    for size in [1024, 2048] {
        print("\nBenchmarking Mojo 0.9 {size}x{size} matrix multiplication...")
        let result = benchmark_mojo_matmul(matrix_size=size, iterations=100)
        if result.is_error() {
            print("Benchmark failed: {result.error()}")
        }
    }
}

import numpy as np
import time
import sys
from typing import Tuple, List
# Mojo FFI imports (assumes mojo Python package is installed)
try:
    from mojo import Python as Mojo
    from mojo.tensor import Tensor as MojoTensor
    MOJO_AVAILABLE = True
except ImportError:
    print("Warning: Mojo Python bindings not available. Running pure Python fallback.")
    MOJO_AVAILABLE = False

class LinearRegressionTrainer:
    def __init__(self, input_dim: int = 10, learning_rate: float = 0.01, use_mojo: bool = True):
        self.input_dim = input_dim
        self.lr = learning_rate
        self.use_mojo = use_mojo and MOJO_AVAILABLE
        self.weights = np.random.randn(input_dim).astype(np.float32)
        self.bias = np.float32(0.0)

        if self.use_mojo:
            # Initialize Mojo weights tensor
            self.mojo_weights = MojoTensor.from_numpy(self.weights)
            self.mojo_bias = MojoTensor.from_numpy(np.array([self.bias]))

    def generate_batch(self, batch_size: int = 32) -> Tuple[np.ndarray, np.ndarray]:
        """Generate synthetic linear regression batch with Gaussian noise."""
        try:
            np.random.seed(int(time.time()))
            X = np.random.randn(batch_size, self.input_dim).astype(np.float32)
            # True weights: [1.0, 2.0, ..., input_dim]
            true_weights = np.arange(1, self.input_dim + 1, dtype=np.float32)
            y = X @ true_weights + 0.1 * np.random.randn(batch_size).astype(np.float32)
            return X, y
        except Exception as e:
            print(f"Batch generation error: {e}")
            return np.array([]), np.array([])

    def forward_pass_mojo(self, X: np.ndarray) -> np.ndarray:
        """Forward pass using Mojo 0.9 SIMD-optimized linear layer."""
        if not self.use_mojo:
            return self.forward_pass_python(X)
        try:
            # Convert input to Mojo tensor
            mojo_X = MojoTensor.from_numpy(X)
            # Linear forward: X @ weights + bias
            logits = mojo_X @ self.mojo_weights + self.mojo_bias
            return logits.to_numpy()
        except Exception as e:
            print(f"Mojo forward pass error: {e}. Falling back to Python.")
            return self.forward_pass_python(X)

    def forward_pass_python(self, X: np.ndarray) -> np.ndarray:
        """Fallback pure Python forward pass."""
        return X @ self.weights + self.bias

    def train_step(self, X: np.ndarray, y: np.ndarray) -> float:
        """Single training step with MSE loss."""
        try:
            # Forward pass
            preds = self.forward_pass_mojo(X) if self.use_mojo else self.forward_pass_python(X)
            # MSE loss gradient: 2*(preds - y)/batch_size
            grad = 2 * (preds - y) / X.shape[0]
            # Update weights (SGD)
            self.weights -= self.lr * (X.T @ grad)
            self.bias -= self.lr * grad.sum()

            if self.use_mojo:
                # Sync Mojo weights with updated NumPy weights
                self.mojo_weights = MojoTensor.from_numpy(self.weights)
                self.mojo_bias = MojoTensor.from_numpy(np.array([self.bias]))

            # Return loss for logging
            return np.mean((preds - y) ** 2)
        except Exception as e:
            print(f"Training step error: {e}")
            return 0.0

    def train(self, epochs: int = 10, batch_size: int = 32) -> List[float]:
        """Run full training loop, return loss history."""
        loss_history = []
        for epoch in range(epochs):
            epoch_loss = 0.0
            num_batches = 100
            start = time.perf_counter()
            for _ in range(num_batches):
                X, y = self.generate_batch(batch_size)
                if X.size == 0:
                    print("Skipping empty batch")
                    continue
                loss = self.train_step(X, y)
                epoch_loss += loss
            epoch_time = time.perf_counter() - start
            avg_loss = epoch_loss / num_batches
            loss_history.append(avg_loss)
            print(f"Epoch {epoch+1}/{epochs} | Loss: {avg_loss:.4f} | Time: {epoch_time:.2f}s | Mojo: {self.use_mojo}")
        return loss_history

if __name__ == "__main__":
    # Benchmark training loop with and without Mojo
    for use_mojo in [False, True]:
        print(f"\nTraining linear regression with Mojo: {use_mojo}")
        trainer = LinearRegressionTrainer(input_dim=10, learning_rate=0.01, use_mojo=use_mojo)
        losses = trainer.train(epochs=5, batch_size=64)

import os
import sys

def check_jit_status():
    if sys.version_info < (3, 13):
        return "JIT not available (Python < 3.13)"
    jit_enabled = os.environ.get("PYTHONJIT") == "1"
    return f"JIT {'enabled' if jit_enabled else 'disabled'} (Python 3.13+)"

from simd import SIMD

fn simd_dot(a: SIMD[Float32, 4], b: SIMD[Float32, 4]) -> Float32:
    return (a * b).reduce_add()

import cProfile
import pstats

def profile_training_step(trainer, X, y):
    profiler = cProfile.Profile()
    profiler.enable()
    trainer.train_step(X, y)
    profiler.disable()
    stats = pstats.Stats(profiler)
    stats.sort_stats("cumtime")
    stats.print_stats(10)  # Print top 10 time-consuming functions