import numpy as np
import matplotlib.pyplot as plt

def generate_random_walks_2d_discrete(num_walks, num_steps):
    """
    Generates a specified number of 2D discrete random walks.

    Each walk starts at (0,0). Steps are of size 1 and can be in
    one of four directions: right, left, up, or down.

    Args:
        num_walks (int): The number of random walks to generate.
        num_steps (int): The number of steps each walk should take.

    Returns:
        list: A list of numpy arrays, where each array represents a
              walk and contains (x, y) coordinates for each step.
    """
    walks = []
    # Define possible steps: right, left, up, down
    # Each row is a [dx, dy]
    step_directions = np.array([[1, 0], [-1, 0], [0, 1], [0, -1]])

    for _ in range(num_walks):
        # Initialize walk at (0,0)
        # Shape: (num_steps + 1, 2) to include the starting point
        path = np.zeros((num_steps + 1, 2))

        for i in range(num_steps):
            # Randomly choose one of the four directions
            step = step_directions[np.random.randint(0, 4)]
            # Update position
            path[i+1] = path[i] + step
        walks.append(path)
    return walks

def plot_random_walks(walks, title="2D Discrete Random Walks"):
    """
    Plots a collection of 2D random walks on a single Matplotlib plot.

    Args:
        walks (list): A list of numpy arrays, where each array represents a
                      walk (output from generate_random_walks_2d_discrete).
        title (str): The title for the plot.
    """
    if not walks:
        print("No walks to plot.")
        return

    num_walks = len(walks)
    # Create a colormap to assign a different color to each walk
    colors = plt.cm.viridis(np.linspace(0, 1, num_walks))

    plt.figure(figsize=(10, 10)) # Adjust figure size for better visibility

    for i, walk_path in enumerate(walks):
        # Plot x and y coordinates of the walk
        plt.plot(walk_path[:, 0], walk_path[:, 1], color=colors[i], alpha=0.7, linewidth=0.8)

    plt.title(title)
    plt.xlabel("X position")
    plt.ylabel("Y position")
    plt.grid(True, linestyle='--', alpha=0.6) # Add a grid for readability
    plt.axis('equal') # Ensure steps have equal visual length in x and y
    plt.show()

# --- Main execution ---
if __name__ == '__main__':
    # Parameters for the random walks
    number_of_walks = 1000
    steps_per_walk = 100  # You can adjust the number of steps per walk

    # 1. Generate the random walks
    all_walks = generate_random_walks_2d_discrete(number_of_walks, steps_per_walk)

    # 2. Visualize the walks
    plot_title = f"{number_of_walks} Two-Dimensional Discrete Random Walks ({steps_per_walk} steps each)"
    plot_random_walks(all_walks, title=plot_title)

import numpy as np
import matplotlib.pyplot as plt

def generate_random_walks_2d_continuous(num_walks, num_steps):
    """
    Generates a specified number of 2D continuous random walks.

    Each walk starts at (0,0). Steps have a random direction (0 to 2*pi)
    and a random length (uniformly distributed between 0 and 1).

    Args:
        num_walks (int): The number of random walks to generate.
        num_steps (int): The number of steps each walk should take.

    Returns:
        list: A list of numpy arrays, where each array represents a
              walk and contains (x, y) coordinates for each step.
    """
    walks = []
    for _ in range(num_walks):
        # Initialize walk at (0,0)
        # Shape: (num_steps + 1, 2) to include the starting point
        path = np.zeros((num_steps + 1, 2))

        for i in range(num_steps):
            # 1. Generate a random angle (direction) between 0 and 2*pi
            angle = np.random.uniform(0, 2 * np.pi)

            # 2. Generate a random step size between 0 and 1
            step_size = np.random.uniform(0, 1)

            # 3. Calculate the change in x and y (dx, dy)
            dx = step_size * np.cos(angle)
            dy = step_size * np.sin(angle)

            # 4. Update position
            path[i+1] = path[i] + [dx, dy]
        walks.append(path)
    return walks

def plot_random_walks(walks, title="2D Random Walks"):
    """
    Plots a collection of 2D random walks on a single Matplotlib plot.

    Args:
        walks (list): A list of numpy arrays, where each array represents a
                      walk (output from a generate_random_walks function).
        title (str): The title for the plot.
    """
    if not walks:
        print("No walks to plot.")
        return

    num_walks = len(walks)
    # Create a colormap to assign a different color to each walk
    # Using plt.cm.get_cmap for newer Matplotlib versions if cm.viridis is an issue
    try:
        colors = plt.cm.viridis(np.linspace(0, 1, num_walks))
    except AttributeError: # Fallback for older versions or specific environments
        cmap = plt.get_cmap('magma')
        colors = cmap(np.linspace(0, 1, num_walks))


    plt.figure(figsize=(10, 10)) # Adjust figure size for better visibility

    for i, walk_path in enumerate(walks):
        # Plot x and y coordinates of the walk
        plt.plot(walk_path[:, 0], walk_path[:, 1], color=colors[i], alpha=0.6, linewidth=0.5)

    plt.title(title)
    plt.xlabel("X position")
    plt.ylabel("Y position")
    plt.grid(True, linestyle='--', alpha=0.6) # Add a grid for readability
    plt.axis('equal') # Ensure steps have equal visual length in x and y initially,
                      # though step sizes vary. This maintains aspect ratio.
    plt.show()

# --- Main execution ---
if __name__ == '__main__':
    # Parameters for the random walks
    number_of_walks = 1000
    steps_per_walk = 100  # You can adjust the number of steps per walk

    # 1. Generate the continuous random walks
    all_walks_continuous = generate_random_walks_2d_continuous(number_of_walks, steps_per_walk)

    # 2. Visualize the walks
    plot_title_continuous = (
        f"{number_of_walks} Continuous 2D Random Walks\n"
        f"({steps_per_walk} steps each, random direction & step size [0,1])"
    )
    plot_random_walks(all_walks_continuous, title=plot_title_continuous)

import numpy as np
import matplotlib.pyplot as plt

def estimate_pi_monte_carlo(num_points=1000):
    """
    Estimates Pi using the Monte Carlo method by simulating points
    in a unit square and checking if they fall within a quarter circle.

    Args:
        num_points (int): The total number of random points to generate.

    Returns:
        float: The estimated value of Pi.
        numpy.ndarray: x-coordinates of points inside the quarter circle.
        numpy.ndarray: y-coordinates of points inside the quarter circle.
        numpy.ndarray: x-coordinates of points outside the quarter circle.
        numpy.ndarray: y-coordinates of points outside the quarter circle.
        int: Number of points inside the quarter circle.
    """
    # Generate random points (x, y) within the unit square [0,1] x [0,1]
    # Each point has x and y coordinates between 0 and 1.
    x_coords = np.random.uniform(0, 1, num_points)
    y_coords = np.random.uniform(0, 1, num_points)

    # Calculate the distance of each point from the origin (0,0)
    # A point is inside the quarter circle if x^2 + y^2 <= 1^2 (radius = 1)
    distance_squared = x_coords**2 + y_coords**2
    is_inside_circle = distance_squared <= 1.0

    # Separate points for plotting and counting
    points_inside_x = x_coords[is_inside_circle]
    points_inside_y = y_coords[is_inside_circle]
    points_outside_x = x_coords[~is_inside_circle]
    points_outside_y = y_coords[~is_inside_circle]

    # Count points inside the quarter circle
    num_inside_circle = np.sum(is_inside_circle)

    # Estimate Pi
    # The ratio (num_inside_circle / num_points) approximates (Area of quarter circle / Area of square)
    # Area of quarter circle = (pi * r^2) / 4. Here r=1. So, Area = pi / 4.
    # Area of unit square = 1 * 1 = 1.
    # So, (num_inside_circle / num_points) approx pi / 4.
    # Therefore, pi approx 4 * (num_inside_circle / num_points).
    pi_approximation = 4 * (num_inside_circle / num_points)

    return (pi_approximation,
            points_inside_x, points_inside_y,
            points_outside_x, points_outside_y,
            num_inside_circle)

def plot_pi_estimation(pi_estimate, num_total_points,
                       points_inside_x, points_inside_y,
                       points_outside_x, points_outside_y,
                       num_inside_circle):
    """
    Plots the visualization for the Monte Carlo Pi estimation.
    """
    num_outside_circle = num_total_points - num_inside_circle

    fig, ax = plt.subplots(figsize=(7, 7)) # Create a figure and an axes

    # Plot points inside the quarter circle (green)
    ax.scatter(points_inside_x, points_inside_y, color='green',
               label=f'Inside Circle ({num_inside_circle})', s=10, alpha=0.7)

    # Plot points outside the quarter circle but inside the square (red)
    ax.scatter(points_outside_x, points_outside_y, color='red',
               label=f'Outside Circle ({num_outside_circle})', s=10, alpha=0.7)

    # Draw the arc of the quarter circle (radius 1, from 0 to pi/2)
    theta = np.linspace(0, np.pi/2, 200) # Angles for the arc
    x_circle_arc = np.cos(theta)         # x-coordinates of the arc
    y_circle_arc = np.sin(theta)         # y-coordinates of the arc
    ax.plot(x_circle_arc, y_circle_arc, color='blue', linewidth=2, label='Quarter Circle (r=1)')

    # Set plot properties
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.set_aspect('equal', adjustable='box') # Crucial for correct visual representation
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    title_str = (
        f'Monte Carlo Estimation of $\pi$\n'
        f'{num_total_points} random points\n'
        f'Estimated $\pi \\approx {pi_estimate:.4f}$' # Using LaTeX for Pi
    )
    ax.set_title(title_str)
    ax.legend(loc='upper right')
    ax.grid(True, linestyle='--', alpha=0.5)

    plt.show()

# --- Main execution ---
if __name__ == '__main__':
    # Parameters
    number_of_points = 1000

    # 1. Perform Monte Carlo estimation
    (pi_value_estimate,
     px_in, py_in,
     px_out, py_out,
     n_inside) = estimate_pi_monte_carlo(num_points=number_of_points)

    # 2. Output the results
    print(f"--- Monte Carlo Pi Estimation ---")
    print(f"Total random points generated: {number_of_points}")
    print(f"Points inside the quarter circle: {n_inside}")
    print(f"Points outside the quarter circle (but in square): {number_of_points - n_inside}")
    print(f"Estimated value of Pi: {pi_value_estimate:.5f}")
    print(f"Value of math.pi (for comparison): {np.pi:.5f}")
    print(f"Absolute error: {abs(np.pi - pi_value_estimate):.5f}")


    # 3. Plot the visualization
    plot_pi_estimation(pi_value_estimate, number_of_points,
                       px_in, py_in, px_out, py_out, n_inside)

import numpy as np
import matplotlib.pyplot as plt

# --- Part 1: Core Pi Estimation Function (simplified) ---
def get_pi_estimate_single_run(num_points):
    """
    Estimates Pi using the Monte Carlo method for a single run.
    This function does not produce any plots.

    Args:
        num_points (int): The total number of random points to generate for this run.

    Returns:
        float: The estimated value of Pi for this run.
    """
    # Ensure num_points is positive to avoid division by zero
    if num_points <= 0:
        return 0.0 # Or raise an error, or return np.nan

    # Generate random points (x, y) within the unit square [0,1] x [0,1]
    # Each point has x and y coordinates between 0 and 1.
    x_coords = np.random.uniform(0, 1, num_points)
    y_coords = np.random.uniform(0, 1, num_points)

    # Calculate the distance squared of each point from the origin (0,0)
    # A point is inside the quarter circle if x^2 + y^2 <= 1^2 (radius = 1)
    distance_squared = x_coords**2 + y_coords**2
    is_inside_circle = distance_squared <= 1.0

    # Count points inside the quarter circle
    num_inside_circle = np.sum(is_inside_circle)

    # Estimate Pi
    # The ratio (num_inside_circle / num_points) approximates (Area of quarter circle / Area of square)
    # Area of quarter circle = (pi * r^2) / 4. Here r=1. So, Area = pi / 4.
    # Area of unit square = 1 * 1 = 1.
    # So, (num_inside_circle / num_points) approx pi / 4.
    # Therefore, pi approx 4 * (num_inside_circle / num_points).
    pi_approximation = 4 * (num_inside_circle / num_points)
    return pi_approximation

# --- Part 2: Function to run multiple samples and get stats ---
def run_pi_simulation_series(num_samples, num_points_per_sample):
    """
    Runs a series of Pi estimation simulations and calculates statistics
    on the collected Pi estimates.

    Args:
        num_samples (int): The number of individual Pi estimation runs (samples).
        num_points_per_sample (int): The number of random points used in each sample.

    Returns:
        tuple: (mean_pi_estimate, std_dev_pi_estimate)
               The mean and standard deviation of the Pi estimates from all samples.
    """
    pi_estimates_for_series = []
    for _ in range(num_samples):
        estimate = get_pi_estimate_single_run(num_points_per_sample)
        pi_estimates_for_series.append(estimate)

    # Convert list to NumPy array for easier statistical calculations
    pi_estimates_array = np.array(pi_estimates_for_series)

    # Calculate the mean of the Pi estimates
    mean_pi = np.mean(pi_estimates_array)

    # Calculate the standard deviation of the Pi estimates
    # This represents the spread of the Pi estimates obtained from the samples.
    std_dev_pi = np.std(pi_estimates_array) # ddof=0 by default for np.std

    return mean_pi, std_dev_pi

# --- Main execution block ---
if __name__ == '__main__':
    # Define the different simulation configurations
    # Each dictionary contains a label for plotting, number of samples, and points per sample.
    configurations = [
        {"label": "100 samples\n100 pts/sample", "samples": 100, "points": 100},
        {"label": "1000 samples\n100 pts/sample", "samples": 1000, "points": 100},
        {"label": "100 samples\n1000 pts/sample", "samples": 100, "points": 1000},
        {"label": "1000 samples\n1000 pts/sample", "samples": 1000, "points": 1000},
        {"label": "10000 samples\n1000 pts/sample", "samples": 10000, "points": 1000},
        {"label": "10000 samples\n10000 pts/sample", "samples": 10000, "points": 10000},
        {"label": "1000 samples\n10000 pts/sample", "samples": 1000, "points": 10000},
        {"label": "100 samples\n100000 pts/sample", "samples": 100, "points": 100000},
    ]

    results_means = []      # To store the mean Pi estimate for each configuration
    results_std_devs = []   # To store the std dev of Pi estimates for each configuration
    config_labels_for_plot = [] # To store labels for the plot's x-axis

    print("--- Monte Carlo Pi Estimation: Series Results ---")
    # Loop through each configuration, run simulations, and store results
    for config in configurations:
        mean_pi, std_dev_pi = run_pi_simulation_series(config["samples"], config["points"])

        results_means.append(mean_pi)
        results_std_devs.append(std_dev_pi)
        config_labels_for_plot.append(config["label"])

        # Print results for the current configuration
        # Replacing newline in label for cleaner single-line console output
        print(f"\nConfiguration: {config['label'].replace(chr(10), ' ')}")
        print(f"  Number of Samples: {config['samples']}")
        print(f"  Points per Sample: {config['points']}")
        print(f"  Mean of Pi Estimates (Mean of Sample Means): {mean_pi:.5f}")
        print(f"  Std Dev of Pi Estimates (Sigma of Sample Means): {std_dev_pi:.5f}")

    # --- Part 3: Plotting the summary results ---
    plt.figure(figsize=(12, 7)) # Adjusted figure size for better label display
    x_positions = np.arange(len(config_labels_for_plot)) # X-coordinates for the bars/points

    # Create the error bar plot
    # fmt='o' plots points at the mean values
    # capsize adds caps to the error bars
    plt.errorbar(x_positions, results_means, yerr=results_std_devs,
                 fmt='o', color='dodgerblue', markersize=8,
                 capsize=5, capthick=1, elinewidth=1,
                 label='Mean Estimated $\pi$ (with Std Dev of Sample Means)')

    # Add a horizontal dashed line for the actual value of Pi
    plt.axhline(y=np.pi, color='red', linestyle='--', linewidth=1.5,
                label=f'Actual $\pi \\approx {np.pi:.5f}$')

    # Set plot properties
    plt.xticks(x_positions, config_labels_for_plot, rotation=0, ha='center') # Center labels
    plt.ylabel('Estimated Value of $\pi$')
    plt.xlabel('Simulation Configuration (Samples x Points per Sample)')
    plt.title('Monte Carlo $\pi$ Estimation: Impact of Sample Size and Points per Sample')
    plt.legend(loc='best')
    plt.grid(True, linestyle=':', alpha=0.6, axis='y') # Lighter grid on y-axis
    plt.tight_layout() # Adjust layout to make sure everything fits without overlapping

    plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- Part 1: Core Pi Estimation Function using MC Integration ---
def estimate_pi_mc_integration_single_run(num_points_N):
    """
    Estimates Pi using Monte Carlo integration of sqrt(1-x^2) from 0 to 1.
    The integral of sqrt(1-x^2) dx from 0 to 1 is pi/4.

    Args:
        num_points_N (int): The number of random points (x_i values) to generate.

    Returns:
        float: The estimated value of Pi for this run.
    """
    if num_points_N <= 0:
        # Return NaN or raise error for invalid input
        return np.nan

    # Generate N random numbers uniformly distributed between 0 and 1
    # x_values are in [0, 1), so 1 - x_values**2 is in (0, 1]
    x_values = np.random.uniform(0, 1, num_points_N)

    # Calculate f(x_i) = sqrt(1 - x_i^2) for each x_i
    integrand_values = np.sqrt(1 - x_values**2)

    # The mean of these values approximates the integral of f(x) from 0 to 1
    mean_of_integrand = np.mean(integrand_values)

    # Estimate of pi/4 is mean_of_integrand
    # Estimate of pi is 4 * mean_of_integrand
    pi_estimate = 4 * mean_of_integrand
    return pi_estimate

# --- Part 2: Function to run multiple samples and get stats ---
def run_mc_integration_series(K_num_runs, N_points_per_run):
    """
    Runs a series of Pi estimations using MC integration and calculates statistics.

    Args:
        K_num_runs (int): The number of individual Pi estimation runs (e.g., 20).
                           This is K in the error analysis.
        N_points_per_run (int): The number of random points (N) used in each run.

    Returns:
        tuple: (mean_pi_estimate, std_dev_of_pi_estimates)
               The mean of the K Pi estimates, and the standard deviation of these K Pi estimates.
    """
    pi_estimates_from_runs = []
    for _ in range(K_num_runs):
        estimate = estimate_pi_mc_integration_single_run(N_points_per_run)
        pi_estimates_from_runs.append(estimate)

    pi_estimates_array = np.array(pi_estimates_from_runs)
    
    # Filter out NaN values that might result from N_points_per_run <= 0, if not handled earlier
    pi_estimates_array = pi_estimates_array[~np.isnan(pi_estimates_array)]
    if len(pi_estimates_array) == 0:
        return np.nan, np.nan

    mean_pi = np.mean(pi_estimates_array)
    # Std dev of the K individual Pi estimates (each derived from N points)
    std_dev_pi = np.std(pi_estimates_array) # ddof=0 by default

    return mean_pi, std_dev_pi

# --- Main execution block ---
if __name__ == '__main__':
    K_samples_for_stats = 100  # Number of samples (runs) to average for each N
    N_values = [100, 1000, 10000, 100000, 1000000, 10000000] # Points per Pi estimate

    exact_pi = np.pi
    
    # Store results for plotting
    mean_pi_estimates_for_each_N = []
    # std_dev_of_the_K_pi_estimates_for_each_N = [] # Standard deviation of the K estimates
    absolute_errors_of_mean = []

    print("--- Monte Carlo Integration for Pi: Error Convergence Analysis ---")
    for N_val in N_values:
        # For each N, get the mean of K Pi estimates and the std dev of those K estimates
        mean_pi_of_K_runs, std_dev_of_K_runs = run_mc_integration_series(K_samples_for_stats, N_val)

        if np.isnan(mean_pi_of_K_runs):
            print(f"\nN = {N_val:10,d} points per estimate: Skipped due to invalid input or all NaN results.")
            absolute_errors_of_mean.append(np.nan) # Placeholder for plotting
            mean_pi_estimates_for_each_N.append(np.nan)
            continue

        mean_pi_estimates_for_each_N.append(mean_pi_of_K_runs)
        # std_dev_of_the_K_pi_estimates_for_each_N.append(std_dev_of_K_runs)

        # Absolute error of the mean of K estimates
        error = np.abs(exact_pi - mean_pi_of_K_runs)
        absolute_errors_of_mean.append(error)

        print(f"\nN = {N_val:10,d} points per estimate:")
        print(f"  Number of samples (K) for stats: {K_samples_for_stats}")
        print(f"  Mean of {K_samples_for_stats} Pi Estimates: {mean_pi_of_K_runs:.8f}")
        print(f"  Std Dev of these {K_samples_for_stats} Pi Estimates: {std_dev_of_K_runs:.8f}")
        print(f"  Absolute Error (|Exact Pi - Mean of {K_samples_for_stats} Estimates|): {error:.3e}")

    # --- Part 4: Plotting the error convergence ---
    plt.figure(figsize=(10, 7))

    # Plot actual absolute error vs N
    plt.plot(N_values, absolute_errors_of_mean, marker='o', linestyle='-', color='dodgerblue',
             label=f'Actual Error: |$\pi_{{exact}}$ - Mean of {K_samples_for_stats} Estimates|')

    # Theoretical error scaling:
    # The standard deviation of a single Pi estimate from N points (sigma_P_N) is (4 * sigma_f) / sqrt(N)
    # where sigma_f = sqrt(Var(sqrt(1-x^2))) = sqrt( (2/3) - (pi/4)^2 ).
    sigma_f_val = np.sqrt((2/3) - (np.pi/4)**2) # Approx 0.2231
    # The error of the mean of K_samples_for_stats estimates is expected to be of the order of
    # Standard Error of Mean (SEM) = sigma_P_N / sqrt(K_samples_for_stats)
    # SEM = (4 * sigma_f_val) / (sqrt(N) * sqrt(K_samples_for_stats))
    # So, Error ~ C_theory / sqrt(N), where C_theory = (4 * sigma_f_val) / sqrt(K_samples_for_stats)
    
    C_theory = (4 * sigma_f_val) / np.sqrt(K_samples_for_stats) # Approx 0.1995 for K=20
    
    # Generate theoretical error line
    theoretical_error_line = [C_theory / np.sqrt(n_val) for n_val in N_values if n_val > 0]
    valid_N_values_for_theory = [n_val for n_val in N_values if n_val > 0]

    if theoretical_error_line: # Only plot if we have valid N values
        plt.plot(valid_N_values_for_theory, theoretical_error_line, marker='.', linestyle='--', color='red',
                 label=f'Theoretical Error Trend ($C/\\sqrt{{N}}$)\n$C = (4 \sigma_f) / \sqrt{{{K_samples_for_stats}}} \\approx {C_theory:.4f}$')

    plt.xscale('log')
    plt.yscale('log') # Vertical log scale as requested
    plt.xlabel('Number of Points per $\pi$ Estimate (N)')
    plt.ylabel('Absolute Error (log scale)')
    plt.title(f'Monte Carlo Integration for $\pi$: Error Convergence')
    plt.legend(loc='best')
    plt.grid(True, which="both", linestyle=':', alpha=0.7) # Grid for major and minor ticks
    plt.tight_layout()
    plt.show()

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import pickle
import os # For checking file sizes

def generate_sierpinski_carpet(n_levels):
    """
    Generates a Sierpinski carpet as a NumPy array.
    0 represents the carpet (black), 1 represents holes (white).

    Args:
        n_levels (int): The number of recursive levels (N in 3^N x 3^N).

    Returns:
        numpy.ndarray: A 2D NumPy array representing the Sierpinski carpet.
    """
    if not isinstance(n_levels, int) or n_levels < 0:
        raise ValueError("n_levels (N) must be a non-negative integer.")

    image_side_length = 3**n_levels
    
    # Initialize the carpet: 0 for black (carpet material), 1 for white (holes)
    # Start with a fully black carpet.
    carpet_array = np.zeros((image_side_length, image_side_length), dtype=np.uint8)

    # Recursive helper function to "punch holes"
    def _create_holes_recursive(arr, x_start, y_start, current_side, current_level, max_level):
        # Base case: if we've reached the desired number of iterations, stop
        if current_level >= max_level:
            return

        # Size of the sub-squares (integer division)
        sub_square_side = current_side // 3

        # If sub_square_side is 0, it means current_side was < 3, so no more holes can be punched.
        if sub_square_side == 0:
            return

        # Calculate coordinates for the central hole
        hole_y_start_idx = y_start + sub_square_side
        hole_y_end_idx = y_start + 2 * sub_square_side
        hole_x_start_idx = x_start + sub_square_side
        hole_x_end_idx = x_start + 2 * sub_square_side

        # "Punch" the central hole by setting its pixels to 1 (white)
        arr[hole_y_start_idx:hole_y_end_idx, hole_x_start_idx:hole_x_end_idx] = 1

        # Recursively call for the 8 surrounding sub-squares
        next_level = current_level + 1
        for i in range(3):
            for j in range(3):
                # Skip the central square (which is now a hole)
                if i == 1 and j == 1:
                    continue
                
                new_x_start = x_start + i * sub_square_side
                new_y_start = y_start + j * sub_square_side
                _create_holes_recursive(arr, new_x_start, new_y_start, sub_square_side, next_level, max_level)
    
    # Initial call to the recursive function
    # Starts with level 0, up to n_levels (exclusive for current_level in recursive calls)
    _create_holes_recursive(carpet_array, 0, 0, image_side_length, 0, n_levels)
    
    return carpet_array

# --- Main execution ---
if __name__ == "__main__":
    N = 5  # Number of levels for the Sierpinski carpet (size will be 3^N x 3^N)
    # For N=7, size is 2187x2187 pixels.
    # For N=6, size is 729x729 pixels.
    # For N=5, size is 243x243 pixels.

    image_dimension = 3**N
    base_filename = f"sierpinski_carpet_N{N}_{image_dimension}x{image_dimension}"

    print(f"Generating Sierpinski carpet for N={N} (size: {image_dimension}x{image_dimension})...")
    try:
        sierpinski_array = generate_sierpinski_carpet(N)
        print("Sierpinski carpet array generated successfully.")

        # Option 1: Save as a 1-bit PNG file using Pillow
        png_filepath = base_filename + ".png"
        try:
            pil_image = Image.fromarray(255*sierpinski_array, mode='L')
            pil_image.save(png_filepath, optimize=False) # optimize=True can reduce file size
        except ImportError:
            print("Pillow library (PIL) not found. Cannot save as PNG.")
            print("Please install it using: pip install Pillow")
        except Exception as e:
            print(f"An error occurred while saving PNG: {e}")

    except ValueError as ve:
        print(f"Error in parameter: {ve}")
    except Exception as general_e:
        print(f"An unexpected error occurred: {general_e}")

    print("Beispiel-Plot für N=5:")
    plt.matshow(sierpinski_array, cmap='gray')

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import random
import time
from scipy.stats import linregress # For fitting the line

# --- 1. Load Image and Prepare Data ---
def load_and_prepare_image(image_path):
    """
    Loads a PNG image, converts it to a binary NumPy array,
    and extracts coordinates of fractal pixels.
    Assumes input image has 0 for black (fractal) and 255 for white (background).
    Output binary_image has 1 for fractal, 0 for background.
    """
    try:
        print(f"Loading image from: {image_path}")
        img_pil = Image.open(image_path).convert('L') # Ensure grayscale
        img_array_loaded = np.array(img_pil)
        
        # Convert to binary: 1 for fractal (original black), 0 for background
        # Assuming black (carpet) pixels are 0 and white (holes) are 255 in the input image.
        binary_image = (img_array_loaded == 0).astype(np.uint8) 
        
        fractal_pixel_coords = np.argwhere(binary_image == 1) # Get (row, col) or (y, x)
        
        if fractal_pixel_coords.size == 0:
            raise ValueError("No fractal pixels (value 1 after conversion) found. Check image loading and thresholding.")
            
        print(f"Image loaded successfully: {binary_image.shape[0]}x{binary_image.shape[1]} pixels.")
        print(f"Found {len(fractal_pixel_coords)} fractal pixels.")
        return binary_image, fractal_pixel_coords
    except FileNotFoundError:
        print(f"Error: Image file not found at '{image_path}'. Please ensure the path is correct.")
        return None, None
    except Exception as e:
        print(f"Error loading or processing image: {e}")
        return None, None

# --- 2. Implement Mass Measurement Function ---
def measure_mass_in_square(binary_image_array, center_y, center_x, half_side_n):
    """
    Measures the mass (number of fractal pixels, i.e., sum of 1s) 
    in a square of side length (2*half_side_n + 1) centered at (center_y, center_x).
    Handles image boundaries by summing only within valid image coordinates.
    """
    img_height, img_width = binary_image_array.shape
    
    # Define the bounds of the square neighborhood
    y_min = max(0, center_y - half_side_n)
    y_max = min(img_height, center_y + half_side_n + 1) # Slice end is exclusive
    x_min = max(0, center_x - half_side_n)
    x_max = min(img_width, center_x + half_side_n + 1) # Slice end is exclusive
    
    # Extract the neighborhood and sum the fractal pixels (where value is 1)
    neighborhood = binary_image_array[y_min:y_max, x_min:x_max]
    mass = np.sum(neighborhood)
    return mass

# --- Main Analysis Function ---
def run_mass_radius_analysis(image_path, list_of_n_values, 
                             num_centers_per_exp, num_experiments):
    """
    Performs the mass-radius analysis by sampling points from the fractal.
    Returns lists of L values, mean masses, and std devs of mean masses.
    """
    binary_image, all_fractal_coords = load_and_prepare_image(image_path)
    if binary_image is None or all_fractal_coords.size == 0:
        print("Failed to load or process image data. Aborting analysis.")
        return [], [], []

    img_height, img_width = binary_image.shape
    
    results_L_values = []
    results_mean_of_means_mass = []
    results_std_dev_of_means_mass = []

    total_fractal_pixels = len(all_fractal_coords)

    for n_val in list_of_n_values:
        current_L_value = 2 * n_val + 1
        print(f"\nProcessing scale L = {current_L_value} (n = {n_val})...")
        
        # For this method, any fractal pixel can be a center.
        # The measure_mass_in_square function handles boundaries.
        # We don't need to filter centers unless L is very large relative to image size
        # and we want to avoid too many truncated boxes, but averaging helps.

        if total_fractal_pixels == 0:
            print("  No fractal pixels to sample centers from. Skipping scale.")
            continue
        
        experiment_mean_masses_for_this_L = [] # Store the mean mass from each experiment

        for exp_idx in range(num_experiments):
            batch_masses_for_current_experiment = []
            
            # Sample 'num_centers_per_exp' indices from all_fractal_coords
            # Using np.random.choice for potentially large arrays.
            # With replacement is fine and simpler if num_centers_per_exp can be > total_fractal_pixels (unlikely here).
            # If sampling without replacement, ensure num_centers_per_exp <= total_fractal_pixels.
            
            num_to_sample = min(num_centers_per_exp, total_fractal_pixels)
            if num_to_sample == 0 : continue # Should not happen if total_fractal_pixels > 0

            # replace=True allows sampling more than available if needed, but usually we want unique centers if possible.
            # Let's stick to sampling 'num_to_sample' without replacement if possible, or with if needed.
            use_replacement = num_centers_per_exp > total_fractal_pixels
            
            sampled_indices = np.random.choice(total_fractal_pixels, 
                                               size=num_centers_per_exp, 
                                               replace=use_replacement)
            
            centers_for_this_batch = all_fractal_coords[sampled_indices]

            for center_y, center_x in centers_for_this_batch:
                mass = measure_mass_in_square(binary_image, center_y, center_x, n_val)
                batch_masses_for_current_experiment.append(mass)
            
            if batch_masses_for_current_experiment:
                mean_mass_this_experiment = np.mean(batch_masses_for_current_experiment)
                experiment_mean_masses_for_this_L.append(mean_mass_this_experiment)
        
        if not experiment_mean_masses_for_this_L:
            print(f"  No mass data collected for L={current_L_value}. Possible issue with sampling or small N. Skipping.")
            continue

        # Calculate the mean of these experiment means, and their standard deviation
        final_mean_M_L = np.mean(experiment_mean_masses_for_this_L)
        final_std_dev_M_L = np.std(experiment_mean_masses_for_this_L) # Std dev of the means
        
        results_L_values.append(current_L_value)
        results_mean_of_means_mass.append(final_mean_M_L)
        results_std_dev_of_means_mass.append(final_std_dev_M_L)
        
        print(f"  L={current_L_value}: Mean of ({num_experiments}) Mean Masses = {final_mean_M_L:.2f}, "
              f"StdDev of these Means = {final_std_dev_M_L:.3f}")

    return results_L_values, results_mean_of_means_mass, results_std_dev_of_means_mass

# --- 6. Visualize Results and Fit for Power Law ---
def visualize_and_fit_results(L_values, mean_mass_values, std_dev_of_means):
    if not L_values or not mean_mass_values:
        print("No data available for visualization.")
        return

    # Filter out non-positive values before taking log, if any (e.g. if mass was 0)
    valid_indices = [i for i, m in enumerate(mean_mass_values) if m > 0]
    if not valid_indices:
        print("No positive mean mass values to plot on log scale.")
        return
        
    L_values_log = np.log([L_values[i] for i in valid_indices])
    mean_mass_log = np.log([mean_mass_values[i] for i in valid_indices])
    # Std dev of means can be used for error bars, but interpretation on log plot needs care.
    # For now, we'll plot points and fit.

    plt.figure(figsize=(10, 7))
    plt.plot(L_values_log, mean_mass_log, 'o', markersize=8, alpha=0.7, label='Log(Mean Mass) vs Log(L)')
    
    # Attempt to fit a line to the log-log data to find the slope (Fractal Dimension)
    # Students should be encouraged to inspect the plot and select an appropriate range for fitting.
    # Here, we'll try to fit a central portion or all if few points.
    if len(L_values_log) >= 2: # Need at least 2 points for a line
        # Heuristic for fitting range (e.g., middle 50% of points if enough points)
        num_pts_log = len(L_values_log)
        fit_start_idx = num_pts_log // 4
        fit_end_idx = (3 * num_pts_log // 4) 
        if fit_end_idx <= fit_start_idx : # if too few points, use all available for fit.
             fit_start_idx = 0
             fit_end_idx = num_pts_log
        if fit_end_idx - fit_start_idx < 2 and num_pts_log >=2: # ensure at least 2 points
            fit_start_idx = 0
            fit_end_idx = num_pts_log


        if fit_end_idx - fit_start_idx >= 2:
            try:
                slope, intercept, r_value, p_value, std_err = linregress(
                    L_values_log[fit_start_idx:fit_end_idx], 
                    mean_mass_log[fit_start_idx:fit_end_idx]
                )
                print(f"\n--- Linear Regression Fit (Log-Log Scale) ---")
                print(f"Fit Range (log_L indices): {fit_start_idx} to {fit_end_idx-1}")
                print(f"  Estimated Fractal Dimension (Slope D): {slope:.4f}")
                print(f"  Intercept: {intercept:.4f}")
                print(f"  R-squared: {r_value**2:.4f}")
                print(f"  Standard error of the slope: {std_err:.4f}")
                
                # Plot the fitted line
                L_fit_plot = np.array(L_values_log) # Use all L values for plotting the line extent
                mass_fit_plot = slope * L_fit_plot + intercept
                plt.plot(L_fit_plot, mass_fit_plot, 'r--', linewidth=2,
                         label=f'Fit: D = {slope:.3f} (R²={r_value**2:.3f})')
            except Exception as e:
                print(f"Could not perform linear regression: {e}")
        else:
            print("Not enough data points in selected range for a meaningful linear fit.")
    else:
        print("Not enough data points for linear regression.")

    plt.xlabel('Log(L) (L = side length of square)')
    plt.ylabel('Log(Mean Mass)')
    plt.title(f'Mass-Radius Method (Squares) - Fractal Dimension Estimation')
    plt.legend()
    plt.grid(True, which="both", linestyle=':', alpha=0.7)
    plt.show()

# --- Main Execution ---
if __name__ == '__main__':
    # IMPORTANT: Students should replace this with the actual path to their generated image
    # For N=8 (6561x6561), use: "sierpinski_carpet_N8_6561x6561_grayscale.png"
    # For N=7 (2187x2187), use: "sierpinski_carpet_N7_2187x2187_grayscale.png"
    # Ensure this file exists in the same directory as the script, or provide a full path.
    IMAGE_FILENAME = "sierpinski_carpet_N8_6561x6561.png" 
    
    # Define scales 'n' (half_side_n for square of side 2n+1).
    # L = 2n+1. Example: n=1 -> L=3; n=2 -> L=5; n=4 -> L=9, etc.
    # These values should be chosen based on image size. Max L should be < image_size/2.
    # For 2187x2187 (N=7), max n could be ~500 (L~1001).
    # For 6561x6561 (N=8), max n could be ~1500 (L~3001).
    # Using a geometric-like progression for n:
    n_values_to_test = [1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233] # Fibonacci-like for L spread
    # Max n=233 => L = 467. This should be fine for N=7 and N=8 images.
    # Filter n_values if L gets too large for smaller images if necessary, but this range is likely okay.

    # Monte Carlo parameters
    num_centers_to_sample_per_experiment = 100 # Number of random centers to average mass over
    num_experiments = 20                       # Number of times to repeat the above to get stats on the mean mass

    print("Starting Mass-Radius Analysis for Fractal Dimension...")
    start_time = time.time()
    
    L_results, mean_mass_results, std_dev_results = run_mass_radius_analysis(
        IMAGE_FILENAME, 
        n_values_to_test, 
        num_centers_to_sample_per_experiment, 
        num_experiments
    )
    
    end_time = time.time()
    total_time = end_time - start_time
    print(f"\n--- Analysis Complete ---")
    print(f"Total execution time: {total_time:.2f} seconds ({total_time/60:.2f} minutes).")

    if L_results and mean_mass_results:
        visualize_and_fit_results(L_results, mean_mass_results, std_dev_results)
    else:
        print("Analysis did not yield results to visualize. Please check for errors above.")