# src/mpi_examples/matrix_multiplication.py
from mpi4py import MPI
import numpy as np
import time

class ParallelMatrixMultiplier:
    def __init__(self):
        self.comm = MPI.COMM_WORLD
        self.rank = self.comm.Get_rank()
        self.size = self.comm.Get_size()
    
    def multiply_matrices(self, A, B):
        """Perform parallel matrix multiplication C = A * B"""
        rows_A, cols_A = A.shape
        rows_B, cols_B = B.shape
        
        if cols_A != rows_B:
            raise ValueError("Matrix dimensions incompatible for multiplication")
        
        # Distribute rows of matrix A across processes
        rows_per_process = rows_A // self.size
        start_row = self.rank * rows_per_process
        end_row = start_row + rows_per_process
        
        # Handle remaining rows for the last process
        if self.rank == self.size - 1:
            end_row = rows_A
        
        # Local computation
        local_A = A[start_row:end_row, :]
        local_C = np.dot(local_A, B)
        
        # Gather results from all processes
        C = None
        if self.rank == 0:
            C = np.zeros((rows_A, cols_B))
        
        # Gather local results
        self.comm.Gather(local_C, C, root=0)
        
        return C
    
    def benchmark_performance(self, matrix_sizes):
        """Benchmark matrix multiplication performance"""
        results = {}
        
        for size in matrix_sizes:
            if self.rank == 0:
                print(f"Testing matrix size: {size}x{size}")
            
            # Generate random matrices
            if self.rank == 0:
                A = np.random.rand(size, size)
                B = np.random.rand(size, size)
            else:
                A = None
                B = None
            
            # Broadcast matrices to all processes
            A = self.comm.bcast(A, root=0)
            B = self.comm.bcast(B, root=0)
            
            # Measure execution time
            start_time = time.time()
            C = self.multiply_matrices(A, B)
            end_time = time.time()
            
            # Collect timing information
            execution_time = end_time - start_time
            times = self.comm.gather(execution_time, root=0)
            
            if self.rank == 0:
                max_time = max(times)
                results[size] = {
                    'execution_time': max_time,
                    'processes': self.size,
                    'matrix_size': size
                }
        
        return results

def main():
    multiplier = ParallelMatrixMultiplier()
    
    # Test different matrix sizes
    matrix_sizes = [100, 200, 500, 1000]
    
    if multiplier.rank == 0:
        print(f"Running parallel matrix multiplication with {multiplier.size} processes")
    
    results = multiplier.benchmark_performance(matrix_sizes)
    
    if multiplier.rank == 0:
        print("\nPerformance Results:")
        for size, result in results.items():
            print(f"Matrix {size}x{size}: {result['execution_time']:.4f} seconds")

if __name__ == "__main__":
    main()

# src/mpi_examples/parallel_sort.py
from mpi4py import MPI
import numpy as np
import random

class ParallelSorter:
    def __init__(self):
        self.comm = MPI.COMM_WORLD
        self.rank = self.comm.Get_rank()
        self.size = self.comm.Get_size()
    
    def parallel_quicksort(self, data):
        """Implement parallel quicksort using MPI"""
        if self.size == 1:
            return sorted(data)
        
        # Distribute data among processes
        local_data = self.distribute_data(data)
        
        # Sort local data
        local_data.sort()
        
        # Perform parallel merge
        sorted_data = self.parallel_merge(local_data)
        
        return sorted_data
    
    def distribute_data(self, data):
        """Distribute data evenly among processes"""
        if self.rank == 0:
            # Calculate distribution
            data_size = len(data)
            chunk_size = data_size // self.size
            remainder = data_size % self.size
            
            # Distribute chunks
            for i in range(self.size):
                start = i * chunk_size
                end = start + chunk_size
                
                if i == self.size - 1:
                    end += remainder
                
                chunk = data[start:end]
                
                if i == 0:
                    local_data = chunk
                else:
                    self.comm.send(chunk, dest=i)
        else:
            local_data = self.comm.recv(source=0)
        
        return local_data
    
    def parallel_merge(self, local_data):
        """Perform parallel merge of sorted local data"""
        current_data = local_data
        
        # Merge in pairs
        step = 1
        while step < self.size:
            if self.rank % (2 * step) == 0:
                # This process will receive data
                partner_rank = self.rank + step
                if partner_rank < self.size:
                    partner_data = self.comm.recv(source=partner_rank)
                    current_data = self.merge_sorted_arrays(current_data, partner_data)
            else:
                # This process will send data
                partner_rank = self.rank - step
                self.comm.send(current_data, dest=partner_rank)
                return None
            
            step *= 2
        
        return current_data
    
    def merge_sorted_arrays(self, arr1, arr2):
        """Merge two sorted arrays"""
        result = []
        i = j = 0
        
        while i < len(arr1) and j < len(arr2):
            if arr1[i] <= arr2[j]:
                result.append(arr1[i])
                i += 1
            else:
                result.append(arr2[j])
                j += 1
        
        # Add remaining elements
        result.extend(arr1[i:])
        result.extend(arr2[j:])
        
        return result
    
    def benchmark_sorting(self, data_sizes):
        """Benchmark sorting performance"""
        results = {}
        
        for size in data_sizes:
            if self.rank == 0:
                print(f"Testing data size: {size}")
                # Generate random data
                data = [random.randint(1, 1000) for _ in range(size)]
            else:
                data = None
            
            # Broadcast data to all processes
            data = self.comm.bcast(data, root=0)
            
            # Measure sorting time
            start_time = MPI.Wtime()
            sorted_data = self.parallel_quicksort(data)
            end_time = MPI.Wtime()
            
            # Collect timing information
            execution_time = end_time - start_time
            times = self.comm.gather(execution_time, root=0)
            
            if self.rank == 0:
                max_time = max(times)
                results[size] = {
                    'execution_time': max_time,
                    'processes': self.size,
                    'data_size': size
                }
        
        return results

def main():
    sorter = ParallelSorter()
    
    # Test different data sizes
    data_sizes = [1000, 5000, 10000, 50000]
    
    if sorter.rank == 0:
        print(f"Running parallel quicksort with {sorter.size} processes")
    
    results = sorter.benchmark_sorting(data_sizes)
    
    if sorter.rank == 0:
        print("\nSorting Performance Results:")
        for size, result in results.items():
            print(f"Data size {size}: {result['execution_time']:.4f} seconds")

if __name__ == "__main__":
    main()

# src/mpi_examples/monte_carlo.py
from mpi4py import MPI
import numpy as np
import random
import math

class ParallelMonteCarlo:
    def __init__(self):
        self.comm = MPI.COMM_WORLD
        self.rank = self.rank.Get_rank()
        self.size = self.comm.Get_size()
    
    def estimate_pi(self, num_samples):
        """Estimate π using Monte Carlo method"""
        # Distribute samples among processes
        samples_per_process = num_samples // self.size
        remainder = num_samples % self.size
        
        if self.rank < remainder:
            local_samples = samples_per_process + 1
        else:
            local_samples = samples_per_process
        
        # Generate random points
        points_in_circle = 0
        
        for _ in range(local_samples):
            x = random.uniform(-1, 1)
            y = random.uniform(-1, 1)
            
            if x*x + y*y <= 1:
                points_in_circle += 1
        
        # Gather results from all processes
        total_points_in_circle = self.comm.reduce(points_in_circle, op=MPI.SUM, root=0)
        
        if self.rank == 0:
            pi_estimate = 4 * total_points_in_circle / num_samples
            return pi_estimate
        
        return None
    
    def parallel_integration(self, func, a, b, num_intervals):
        """Perform parallel numerical integration"""
        # Distribute intervals among processes
        intervals_per_process = num_intervals // self.size
        remainder = num_intervals % self.size
        
        if self.rank < remainder:
            local_intervals = intervals_per_process + 1
            start_interval = self.rank * (intervals_per_process + 1)
        else:
            local_intervals = intervals_per_process
            start_interval = self.rank * intervals_per_process + remainder
        
        # Calculate local integral
        h = (b - a) / num_intervals
        local_integral = 0
        
        for i in range(local_intervals):
            x = a + (start_interval + i) * h
            local_integral += func(x) * h
        
        # Gather results from all processes
        total_integral = self.comm.reduce(local_integral, op=MPI.SUM, root=0)
        
        if self.rank == 0:
            return total_integral
        
        return None
    
    def parallel_random_walk(self, num_steps, num_walks):
        """Simulate parallel random walks"""
        # Distribute walks among processes
        walks_per_process = num_walks // self.size
        remainder = num_walks % self.size
        
        if self.rank < remainder:
            local_walks = walks_per_process + 1
        else:
            local_walks = walks_per_process
        
        # Perform local random walks
        local_results = []
        
        for _ in range(local_walks):
            position = 0
            for _ in range(num_steps):
                if random.random() < 0.5:
                    position += 1
                else:
                    position -= 1
            local_results.append(position)
        
        # Gather all results
        all_results = self.comm.gather(local_results, root=0)
        
        if self.rank == 0:
            # Flatten results
            all_positions = []
            for result_list in all_results:
                all_positions.extend(result_list)
            
            return {
                'mean_position': np.mean(all_positions),
                'std_position': np.std(all_positions),
                'final_positions': all_positions
            }
        
        return None

def main():
    monte_carlo = ParallelMonteCarlo()
    
    if monte_carlo.rank == 0:
        print(f"Running Monte Carlo simulations with {monte_carlo.size} processes")
    
    # Estimate π
    pi_estimate = monte_carlo.estimate_pi(1000000)
    if monte_carlo.rank == 0:
        print(f"π estimate: {pi_estimate:.6f}")
        print(f"Error: {abs(pi_estimate - math.pi):.6f}")
    
    # Numerical integration
    def func(x):
        return x * x  # x²
    
    integral = monte_carlo.parallel_integration(func, 0, 1, 1000000)
    if monte_carlo.rank == 0:
        print(f"Integral of x² from 0 to 1: {integral:.6f}")
        print(f"Expected: 1/3 = {1/3:.6f}")
    
    # Random walk simulation
    walk_results = monte_carlo.parallel_random_walk(1000, 10000)
    if monte_carlo.rank == 0:
        print(f"Random walk results:")
        print(f"Mean position: {walk_results['mean_position']:.4f}")
        print(f"Standard deviation: {walk_results['std_position']:.4f}")

if __name__ == "__main__":
    main()

# src/mpi_examples/heat_equation.py
from mpi4py import MPI
import numpy as np
import matplotlib.pyplot as plt

class ParallelHeatSolver:
    def __init__(self, nx, ny, nt, alpha=0.1):
        self.comm = MPI.COMM_WORLD
        self.rank = self.comm.Get_rank()
        self.size = self.comm.Get_size()
        
        self.nx = nx
        self.ny = ny
        self.nt = nt
        self.alpha = alpha
        
        # Calculate grid spacing
        self.dx = 1.0 / (nx - 1)
        self.dy = 1.0 / (ny - 1)
        self.dt = 0.01
        
        # Distribute grid points among processes
        self.setup_grid_distribution()
    
    def setup_grid_distribution(self):
        """Setup grid distribution among processes"""
        # Distribute rows among processes
        self.rows_per_process = self.ny // self.size
        self.remainder_rows = self.ny % self.size
        
        if self.rank < self.remainder_rows:
            self.local_rows = self.rows_per_process + 1
            self.start_row = self.rank * (self.rows_per_process + 1)
        else:
            self.local_rows = self.rows_per_process
            self.start_row = self.rank * self.rows_per_process + self.remainder_rows
        
        self.end_row = self.start_row + self.local_rows
        
        # Add ghost cells for communication
        self.local_grid = np.zeros((self.local_rows + 2, self.nx))
    
    def initialize_grid(self):
        """Initialize the temperature grid"""
        # Set initial conditions
        for i in range(self.local_rows + 2):
            for j in range(self.nx):
                x = j * self.dx
                y = (self.start_row + i - 1) * self.dy
                
                # Initial temperature distribution
                if 0.3 <= x <= 0.7 and 0.3 <= y <= 0.7:
                    self.local_grid[i, j] = 100.0
                else:
                    self.local_grid[i, j] = 0.0
    
    def update_boundary_conditions(self):
        """Update boundary conditions"""
        # Fixed boundary conditions
        if self.start_row == 0:
            self.local_grid[0, :] = 0.0  # Bottom boundary
        
        if self.end_row == self.ny:
            self.local_grid[-1, :] = 0.0  # Top boundary
        
        # Left and right boundaries
        self.local_grid[:, 0] = 0.0
        self.local_grid[:, -1] = 0.0
    
    def exchange_ghost_cells(self):
        """Exchange ghost cells between neighboring processes"""
        # Send to upper neighbor
        if self.rank < self.size - 1:
            self.comm.Send(self.local_grid[-2, :], dest=self.rank + 1, tag=0)
        
        # Send to lower neighbor
        if self.rank > 0:
            self.comm.Send(self.local_grid[1, :], dest=self.rank - 1, tag=1)
        
        # Receive from upper neighbor
        if self.rank < self.size - 1:
            self.comm.Recv(self.local_grid[-1, :], source=self.rank + 1, tag=1)
        
        # Receive from lower neighbor
        if self.rank > 0:
            self.comm.Recv(self.local_grid[0, :], source=self.rank - 1, tag=0)
    
    def solve_heat_equation(self):
        """Solve the heat equation using finite differences"""
        self.initialize_grid()
        
        for t in range(self.nt):
            # Exchange ghost cells
            self.exchange_ghost_cells()
            
            # Update interior points
            new_grid = self.local_grid.copy()
            
            for i in range(1, self.local_rows + 1):
                for j in range(1, self.nx - 1):
                    # Finite difference approximation
                    d2u_dx2 = (self.local_grid[i, j+1] - 2*self.local_grid[i, j] + 
                              self.local_grid[i, j-1]) / (self.dx**2)
                    d2u_dy2 = (self.local_grid[i+1, j] - 2*self.local_grid[i, j] + 
                              self.local_grid[i-1, j]) / (self.dy**2)
                    
                    new_grid[i, j] = self.local_grid[i, j] + self.alpha * self.dt * (d2u_dx2 + d2u_dy2)
            
            self.local_grid = new_grid
            
            # Update boundary conditions
            self.update_boundary_conditions()
            
            # Synchronize all processes
            self.comm.Barrier()
        
        return self.local_grid[1:-1, :]  # Return without ghost cells
    
    def gather_solution(self):
        """Gather solution from all processes"""
        if self.rank == 0:
            full_solution = np.zeros((self.ny, self.nx))
        else:
            full_solution = None
        
        # Gather local solutions
        self.comm.Gather(self.local_grid[1:-1, :], full_solution, root=0)
        
        return full_solution
    
    def visualize_solution(self, solution):
        """Visualize the temperature distribution"""
        if self.rank == 0:
            plt.figure(figsize=(10, 8))
            plt.imshow(solution, cmap='hot', origin='lower')
            plt.colorbar(label='Temperature')
            plt.title('Temperature Distribution')
            plt.xlabel('X')
            plt.ylabel('Y')
            plt.show()

def main():
    # Problem parameters
    nx, ny = 100, 100  # Grid size
    nt = 1000  # Time steps
    
    solver = ParallelHeatSolver(nx, ny, nt)
    
    if solver.rank == 0:
        print(f"Solving heat equation with {solver.size} processes")
        print(f"Grid size: {nx}x{ny}")
        print(f"Time steps: {nt}")
    
    # Solve the heat equation
    start_time = MPI.Wtime()
    solution = solver.solve_heat_equation()
    end_time = MPI.Wtime()
    
    if solver.rank == 0:
        print(f"Solution completed in {end_time - start_time:.4f} seconds")
    
    # Gather and visualize solution
    full_solution = solver.gather_solution()
    solver.visualize_solution(full_solution)

if __name__ == "__main__":
    main()

# src/benchmarks/scalability_tests.py
from mpi4py import MPI
import numpy as np
import time
import matplotlib.pyplot as plt

class ScalabilityAnalyzer:
    def __init__(self):
        self.comm = MPI.COMM_WORLD
        self.rank = self.comm.Get_rank()
        self.size = self.comm.Get_size()
    
    def test_matrix_multiplication_scalability(self, matrix_sizes, process_counts):
        """Test matrix multiplication scalability"""
        results = {}
        
        for size in matrix_sizes:
            results[size] = {}
            
            for processes in process_counts:
                if self.size == processes:
                    # Generate matrices
                    if self.rank == 0:
                        A = np.random.rand(size, size)
                        B = np.random.rand(size, size)
                    else:
                        A = None
                        B = None
                    
                    # Broadcast matrices
                    A = self.comm.bcast(A, root=0)
                    B = self.comm.bcast(B, root=0)
                    
                    # Measure execution time
                    start_time = MPI.Wtime()
                    
                    # Perform parallel multiplication
                    rows_per_process = size // processes
                    start_row = self.rank * rows_per_process
                    end_row = start_row + rows_per_process
                    
                    if self.rank == processes - 1:
                        end_row = size
                    
                    local_A = A[start_row:end_row, :]
                    local_C = np.dot(local_A, B)
                    
                    # Gather results
                    C = None
                    if self.rank == 0:
                        C = np.zeros((size, size))
                    
                    self.comm.Gather(local_C, C, root=0)
                    
                    end_time = MPI.Wtime()
                    
                    if self.rank == 0:
                        execution_time = end_time - start_time
                        results[size][processes] = execution_time
        
        return results
    
    def calculate_speedup(self, results):
        """Calculate speedup for different process counts"""
        speedup_results = {}
        
        for size, process_times in results.items():
            speedup_results[size] = {}
            baseline_time = process_times[1]  # Single process time
            
            for processes, execution_time in process_times.items():
                speedup = baseline_time / execution_time
                speedup_results[size][processes] = speedup
        
        return speedup_results
    
    def calculate_efficiency(self, speedup_results):
        """Calculate parallel efficiency"""
        efficiency_results = {}
        
        for size, speedups in speedup_results.items():
            efficiency_results[size] = {}
            
            for processes, speedup in speedups.items():
                efficiency = speedup / processes
                efficiency_results[size][processes] = efficiency
        
        return efficiency_results
    
    def plot_scalability_results(self, results, output_file='scalability_results.png'):
        """Plot scalability results"""
        if self.rank == 0:
            fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
            
            # Plot execution times
            for size, process_times in results.items():
                processes = list(process_times.keys())
                times = list(process_times.values())
                ax1.plot(processes, times, marker='o', label=f'Matrix {size}x{size}')
            
            ax1.set_xlabel('Number of Processes')
            ax1.set_ylabel('Execution Time (seconds)')
            ax1.set_title('Execution Time vs Number of Processes')
            ax1.legend()
            ax1.grid(True)
            
            # Calculate and plot speedup
            speedup_results = self.calculate_speedup(results)
            
            for size, speedups in speedup_results.items():
                processes = list(speedups.keys())
                speedups_list = list(speedups.values())
                ax2.plot(processes, speedups_list, marker='s', label=f'Matrix {size}x{size}')
            
            # Plot ideal speedup line
            max_processes = max(max(process_times.keys()) for process_times in results.values())
            ideal_speedup = list(range(1, max_processes + 1))
            ax2.plot(ideal_speedup, ideal_speedup, 'k--', label='Ideal Speedup')
            
            ax2.set_xlabel('Number of Processes')
            ax2.set_ylabel('Speedup')
            ax2.set_title('Speedup vs Number of Processes')
            ax2.legend()
            ax2.grid(True)
            
            plt.tight_layout()
            plt.savefig(output_file, dpi=300, bbox_inches='tight')
            plt.show()

def main():
    analyzer = ScalabilityAnalyzer()
    
    if analyzer.rank == 0:
        print("Running scalability analysis...")
    
    # Test parameters
    matrix_sizes = [100, 200, 500]
    process_counts = [1, 2, 4, 8]
    
    # Run scalability tests
    results = analyzer.test_matrix_multiplication_scalability(matrix_sizes, process_counts)
    
    # Plot results
    analyzer.plot_scalability_results(results)
    
    if analyzer.rank == 0:
        print("Scalability analysis completed!")

if __name__ == "__main__":
    main()

# src/utils/communication_patterns.py
from mpi4py import MPI
import numpy as np

class CommunicationPatterns:
    def __init__(self):
        self.comm = MPI.COMM_WORLD
        self.rank = self.rank.Get_rank()
        self.size = self.comm.Get_size()
    
    def ring_allreduce(self, data):
        """Implement ring allreduce algorithm"""
        local_data = data.copy()
        step_size = 1
        
        while step_size < self.size:
            # Calculate partner rank
            partner_rank = (self.rank + step_size) % self.size
            
            # Send data to partner
            self.comm.Send(local_data, dest=partner_rank, tag=0)
            
            # Receive data from partner
            received_data = np.zeros_like(local_data)
            self.comm.Recv(received_data, source=partner_rank, tag=0)
            
            # Add received data to local data
            local_data += received_data
            
            step_size *= 2
        
        return local_data
    
    def butterfly_allreduce(self, data):
        """Implement butterfly allreduce algorithm"""
        local_data = data.copy()
        step_size = self.size // 2
        
        while step_size > 0:
            # Calculate partner rank
            partner_rank = self.rank ^ step_size
            
            # Send data to partner
            self.comm.Send(local_data, dest=partner_rank, tag=0)
            
            # Receive data from partner
            received_data = np.zeros_like(local_data)
            self.comm.Recv(received_data, source=partner_rank, tag=0)
            
            # Add received data to local data
            local_data += received_data
            
            step_size //= 2
        
        return local_data
    
    def scatter_gather_pattern(self, data):
        """Demonstrate scatter-gather pattern"""
        if self.rank == 0:
            # Scatter data to all processes
            chunk_size = len(data) // self.size
            scattered_data = []
            
            for i in range(self.size):
                start = i * chunk_size
                end = start + chunk_size
                if i == self.size - 1:
                    end = len(data)
                
                chunk = data[start:end]
                scattered_data.append(chunk)
        else:
            scattered_data = None
        
        # Receive scattered data
        local_data = self.comm.scatter(scattered_data, root=0)
        
        # Process local data (example: square each element)
        processed_data = [x**2 for x in local_data]
        
        # Gather processed data
        gathered_data = self.comm.gather(processed_data, root=0)
        
        if self.rank == 0:
            # Flatten gathered data
            result = []
            for chunk in gathered_data:
                result.extend(chunk)
            return result
        
        return None
    
    def pipeline_pattern(self, data, stages):
        """Demonstrate pipeline pattern"""
        if self.rank == 0:
            # First stage: data generation
            processed_data = data
        else:
            processed_data = None
        
        # Pass data through pipeline stages
        for stage in range(stages):
            if self.rank == stage:
                # Process data at this stage
                if processed_data is not None:
                    processed_data = [x * 2 for x in processed_data]  # Example processing
                
                # Send to next stage
                if stage < self.size - 1:
                    self.comm.send(processed_data, dest=stage + 1, tag=stage)
            
            if self.rank == stage + 1:
                # Receive from previous stage
                processed_data = self.comm.recv(source=stage, tag=stage)
        
        return processed_data

def main():
    patterns = CommunicationPatterns()
    
    if patterns.rank == 0:
        print(f"Testing communication patterns with {patterns.size} processes")
    
    # Test data
    test_data = np.array([1, 2, 3, 4, 5, 6, 7, 8])
    
    # Test ring allreduce
    if patterns.rank == 0:
        print("Testing ring allreduce...")
    
    ring_result = patterns.ring_allreduce(test_data)
    
    if patterns.rank == 0:
        print(f"Ring allreduce result: {ring_result}")
    
    # Test butterfly allreduce
    if patterns.rank == 0:
        print("Testing butterfly allreduce...")
    
    butterfly_result = patterns.butterfly_allreduce(test_data)
    
    if patterns.rank == 0:
        print(f"Butterfly allreduce result: {butterfly_result}")
    
    # Test scatter-gather pattern
    if patterns.rank == 0:
        print("Testing scatter-gather pattern...")
        scatter_data = list(range(16))
    else:
        scatter_data = None
    
    scatter_result = patterns.scatter_gather_pattern(scatter_data)
    
    if patterns.rank == 0:
        print(f"Scatter-gather result: {scatter_result}")

if __name__ == "__main__":
    main()