#Advanced Features

import numpy as np

# Set random seed for reproducible results
np.random.seed(42)

# Basic random number generation
print("Basic random number generation:")
print(f"Single random number: {np.random.random()}")
print(f"Random integer: {np.random.randint(1, 10)}")
print(f"Random array: {np.random.random(5)}")
print(f"Random integer array: {np.random.randint(1, 10, 5)}")
print()

# Random numbers from different distributions
print("Random numbers from different distributions:")
print(f"Normal distribution: {np.random.normal(0, 1, 5)}")
print(f"Uniform distribution: {np.random.uniform(-1, 1, 5)}")
print(f"Exponential distribution: {np.random.exponential(2, 5)}")
print(f"Poisson distribution: {np.random.poisson(3, 5)}")
print()

# Multidimensional random arrays
print("Multidimensional random arrays:")
random_2d = np.random.random((3, 4))
print(random_2d)
print()

# Random selection from array
array = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
print("Random selection from array:")
print(f"Select 1: {np.random.choice(array)}")
print(f"Select 3: {np.random.choice(array, 3)}")
print(f"Select 3 without replacement: {np.random.choice(array, 3, replace=False)}")
print(f"Weighted selection: {np.random.choice(array, 3, p=[0.1]*10)}")

import numpy as np

# Use new random number generator
rng = np.random.default_rng(42)

print("Modern random number generator:")
print(f"Random number: {rng.random()}")
print(f"Normal distribution: {rng.normal(0, 1, 5)}")
print(f"Integers: {rng.integers(1, 10, 5)}")
print()

# Random permutation
array = np.arange(10)
print(f"Original array: {array}")
rng.shuffle(array)  # Shuffle in place
print(f"After shuffle: {array}")

# Generate permutation
original = np.arange(10)
permutation = rng.permutation(original)
print(f"Original array: {original}")
print(f"New permutation: {permutation}")
print()

# Random sampling
data = np.arange(100)
sample = rng.choice(data, size=10, replace=False)
print(f"Random sample: {sample}")

import numpy as np

# Monte Carlo method to estimate π
def estimate_pi(n_samples):
    rng = np.random.default_rng(42)
    # Generate random points in unit square
    x = rng.uniform(-1, 1, n_samples)
    y = rng.uniform(-1, 1, n_samples)
    
    # Count points inside unit circle
    inside_circle = (x**2 + y**2) <= 1
    pi_estimate = 4 * np.sum(inside_circle) / n_samples
    
    return pi_estimate

print("Monte Carlo estimation of π:")
for n in [1000, 10000, 100000, 1000000]:
    pi_est = estimate_pi(n)
    error = abs(pi_est - np.pi)
    print(f"Samples: {n:7d}, π estimate: {pi_est:.6f}, Error: {error:.6f}")
print(f"Actual π value: {np.pi:.6f}")
print()

# 2D Random Walk
def random_walk_2d(steps):
    rng = np.random.default_rng(42)
    # Random directions: up, down, left, right
    directions = np.array([[0, 1], [0, -1], [-1, 0], [1, 0]])
    
    # Randomly choose directions
    chosen_directions = rng.choice(4, steps)
    moves = directions[chosen_directions]
    
    # Calculate positions
    positions = np.cumsum(moves, axis=0)
    
    return positions

print("2D Random Walk:")
walk = random_walk_2d(20)
print(f"Starting position: (0, 0)")
print(f"Final position: ({walk[-1, 0]}, {walk[-1, 1]})")
print(f"Maximum distance: {np.max(np.sqrt(np.sum(walk**2, axis=1))):.2f}")
print(f"Average distance: {np.mean(np.sqrt(np.sum(walk**2, axis=1))):.2f}")

import numpy as np

# Create polynomial p(x) = 2x³ + 3x² - x + 5
# Coefficients in descending order
coeffs = [2, 3, -1, 5]
print(f"Polynomial coefficients: {coeffs}")
print("Polynomial: 2x³ + 3x² - x + 5")
print()

# Evaluate polynomial
x_values = np.array([0, 1, 2, -1])
results = np.polyval(coeffs, x_values)
print(f"x values: {x_values}")
print(f"p(x): {results}")
print()

# Verify calculations
for x, result in zip(x_values, results):
    manual_calc = 2*x**3 + 3*x**2 - x + 5
    print(f"p({x}) = {result}, Manual calculation: {manual_calc}")
print()

# Polynomial roots
roots = np.roots(coeffs)
print(f"Polynomial roots: {roots}")
print("Verify roots:")
for i, root in enumerate(roots):
    value = np.polyval(coeffs, root)
    print(f"  Root {i+1}: {root:.4f}, p(root) = {value:.2e}")

import numpy as np

# Define two polynomials
# p1(x) = x² + 2x + 1 = (x + 1)²
# p2(x) = x - 1
p1 = [1, 2, 1]
p2 = [1, -1]

print("Polynomial operations:")
print(f"p1(x) = x² + 2x + 1")
print(f"p2(x) = x - 1")
print()

# Polynomial addition
sum_poly = np.polyadd(p1, p2)
print(f"p1 + p2 = {sum_poly}")
print(f"i.e.: x² + 3x")
print()

# Polynomial subtraction
diff_poly = np.polysub(p1, p2)
print(f"p1 - p2 = {diff_poly}")
print(f"i.e.: x² + x + 2")
print()

# Polynomial multiplication
mul_poly = np.polymul(p1, p2)
print(f"p1 * p2 = {mul_poly}")
print(f"i.e.: x³ + x² - x - 1")
print()

# Polynomial division
quotient, remainder = np.polydiv(p1, p2)
print(f"p1 / p2:")
print(f"  Quotient: {quotient}")
print(f"  Remainder: {remainder}")
print(f"  i.e.: (x² + 2x + 1) = (x - 1) * (x + 3) + 4")
print()

# Polynomial derivative
derivative = np.polyder(p1)
print(f"Derivative of p1: {derivative}")
print(f"i.e.: 2x + 2")
print()

# Polynomial integral
integral = np.polyint(p1)
print(f"Integral of p1: {integral}")
print(f"i.e.: x³/3 + x² + x + C")

import numpy as np

# Generate noisy data
np.random.seed(42)
x_data = np.linspace(0, 10, 20)
# True function: y = 2x² - 3x + 1
y_true = 2 * x_data**2 - 3 * x_data + 1
y_noisy = y_true + np.random.normal(0, 5, len(x_data))

print("Polynomial fitting:")
print(f"Number of data points: {len(x_data)}")
print(f"True function: y = 2x² - 3x + 1")
print()

# Fit with different degrees
for degree in [1, 2, 3, 5]:
    coeffs = np.polyfit(x_data, y_noisy, degree)
    
    # Calculate goodness of fit
    y_fit = np.polyval(coeffs, x_data)
    r_squared = 1 - np.sum((y_noisy - y_fit)**2) / np.sum((y_noisy - np.mean(y_noisy))**2)
    
    print(f"Degree {degree} polynomial fit:")
    print(f"  Coefficients: {coeffs}")
    print(f"  R²: {r_squared:.4f}")
    
    if degree == 2:
        print(f"  Fit function: {coeffs[0]:.2f}x² + {coeffs[1]:.2f}x + {coeffs[2]:.2f}")
        print(f"  Compare with true: 2x² - 3x + 1")
    print()

# Predict new data points
x_new = np.array([2.5, 7.5])
best_fit_coeffs = np.polyfit(x_data, y_noisy, 2)  # Use quadratic
y_pred = np.polyval(best_fit_coeffs, x_new)
y_true_new = 2 * x_new**2 - 3 * x_new + 1

print("Prediction results:")
for x, pred, true in zip(x_new, y_pred, y_true_new):
    print(f"x = {x}: Predicted = {pred:.2f}, True = {true:.2f}, Error = {abs(pred-true):.2f}")

import numpy as np

# Create signal
t = np.linspace(0, 1, 1000, endpoint=False)
freq1, freq2 = 5, 10  # 5Hz and 10Hz
signal = np.sin(2 * np.pi * freq1 * t) + 0.5 * np.sin(2 * np.pi * freq2 * t)

print("Fourier Transform Analysis:")
print(f"Signal length: {len(signal)}")
print(f"Sampling frequency: {1/(t[1]-t[0]):.0f} Hz")
print(f"Signal contains frequencies: {freq1} Hz and {freq2} Hz")
print()

# Calculate FFT
fft_result = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(signal), t[1] - t[0])

# Calculate amplitude spectrum
amplitude = np.abs(fft_result)
phase = np.angle(fft_result)

# Consider only positive frequencies
positive_freq_idx = freqs > 0
positive_freqs = freqs[positive_freq_idx]
positive_amplitude = amplitude[positive_freq_idx]

# Find main frequency components
peak_indices = np.where(positive_amplitude > 0.1 * np.max(positive_amplitude))[0]
peak_freqs = positive_freqs[peak_indices]
peak_amplitudes = positive_amplitude[peak_indices]

print("Detected frequency components:")
for freq, amp in zip(peak_freqs, peak_amplitudes):
    if freq < 50:  # Only show low frequency components
        print(f"  Frequency: {freq:.1f} Hz, Amplitude: {amp:.1f}")
print()

# Inverse Fourier Transform
reconstructed = np.fft.ifft(fft_result)
reconstruction_error = np.mean(np.abs(signal - np.real(reconstructed)))
print(f"Reconstruction error: {reconstruction_error:.2e}")
print()

# Frequency domain filtering example
# Create low-pass filter (keep frequencies below 15Hz)
cutoff_freq = 15
filtered_fft = fft_result.copy()
filtered_fft[np.abs(freqs) > cutoff_freq] = 0

# Inverse transform to get filtered signal
filtered_signal = np.real(np.fft.ifft(filtered_fft))

print(f"Low-pass filter (cutoff: {cutoff_freq} Hz):")
print(f"Original signal energy: {np.sum(signal**2):.2f}")
print(f"Filtered energy: {np.sum(filtered_signal**2):.2f}")
print(f"Energy retention: {np.sum(filtered_signal**2)/np.sum(signal**2)*100:.1f}%")

import numpy as np

# Create 2D signal (simulated image)
x = np.linspace(0, 1, 64)
y = np.linspace(0, 1, 64)
X, Y = np.meshgrid(x, y)

# Create image with different frequency components
image = (np.sin(2 * np.pi * 5 * X) * np.cos(2 * np.pi * 3 * Y) + 
         0.5 * np.sin(2 * np.pi * 10 * X) * np.sin(2 * np.pi * 8 * Y))

print("2D Fourier Transform:")
print(f"Image size: {image.shape}")
print(f"Image statistics: mean={np.mean(image):.3f}, std={np.std(image):.3f}")
print()

# 2D FFT
fft_2d = np.fft.fft2(image)
fft_shifted = np.fft.fftshift(fft_2d)  # Shift zero frequency to center

# Calculate amplitude spectrum
amplitude_spectrum = np.abs(fft_shifted)
log_amplitude = np.log(1 + amplitude_spectrum)

print("Frequency domain analysis:")
print(f"Maximum amplitude: {np.max(amplitude_spectrum):.1f}")
print(f"Spectrum center position: ({image.shape[0]//2}, {image.shape[1]//2})")
print()

# Inverse transform verification
reconstructed = np.real(np.fft.ifft2(fft_2d))
reconstruction_error = np.mean(np.abs(image - reconstructed))
print(f"Reconstruction error: {reconstruction_error:.2e}")
print()

# Frequency domain filtering
# Create low-pass filter
center_x, center_y = image.shape[0]//2, image.shape[1]//2
y_idx, x_idx = np.ogrid[:image.shape[0], :image.shape[1]]
distance = np.sqrt((x_idx - center_x)**2 + (y_idx - center_y)**2)

# Gaussian low-pass filter
sigma = 10
gaussian_filter = np.exp(-(distance**2) / (2 * sigma**2))

# Apply filter
filtered_fft = fft_shifted * gaussian_filter
filtered_image = np.real(np.fft.ifft2(np.fft.ifftshift(filtered_fft)))

print("Gaussian low-pass filter:")
print(f"Filter sigma: {sigma}")
print(f"Original image energy: {np.sum(image**2):.2f}")
print(f"Filtered energy: {np.sum(filtered_image**2):.2f}")
print(f"Energy retention: {np.sum(filtered_image**2)/np.sum(image**2)*100:.1f}%")

import numpy as np
import tempfile
import os

# Create sample data
data = np.random.random((5, 3))
print("Original data:")
print(data)
print()

# Save to text file
with tempfile.NamedTemporaryFile(mode='w', suffix='.txt', delete=False) as f:
    txt_filename = f.name
    
np.savetxt(txt_filename, data, fmt='%.6f', delimiter=',')
print(f"Data saved to: {txt_filename}")

# Load from text file
loaded_data = np.loadtxt(txt_filename, delimiter=',')
print("Data loaded from file:")
print(loaded_data)
print(f"Data matches: {np.allclose(data, loaded_data)}")
print()

# Save multiple arrays
array1 = np.array([1, 2, 3, 4, 5])
array2 = np.array([10, 20, 30, 40, 50])
header_info = "Array1,Array2"

combined_data = np.column_stack([array1, array2])
with tempfile.NamedTemporaryFile(mode='w', suffix='.csv', delete=False) as f:
    csv_filename = f.name

np.savetxt(csv_filename, combined_data, fmt='%d', delimiter=',', header=header_info)
print(f"CSV file saved to: {csv_filename}")

# Load CSV file
loaded_csv = np.loadtxt(csv_filename, delimiter=',', skiprows=1)
print("Data loaded from CSV:")
print(loaded_csv)
print()

# Cleanup temporary files
os.unlink(txt_filename)
os.unlink(csv_filename)

import numpy as np
import tempfile
import os

# Create test data
array_int = np.array([1, 2, 3, 4, 5], dtype=np.int32)
array_float = np.array([1.1, 2.2, 3.3, 4.4, 5.5], dtype=np.float64)
array_2d = np.random.random((3, 4))

print("Binary file operations:")
print(f"Integer array: {array_int}")
print(f"Float array: {array_float}")
print(f"2D array shape: {array_2d.shape}")
print()

# Save single array to binary file
with tempfile.NamedTemporaryFile(suffix='.npy', delete=False) as f:
    npy_filename = f.name

np.save(npy_filename, array_2d)
print(f"2D array saved to: {npy_filename}")

# Load binary file
loaded_2d = np.load(npy_filename)
print(f"Loaded array shape: {loaded_2d.shape}")
print(f"Data matches: {np.allclose(array_2d, loaded_2d)}")
print()

# Save multiple arrays to one file
with tempfile.NamedTemporaryFile(suffix='.npz', delete=False) as f:
    npz_filename = f.name

np.savez(npz_filename, 
         integers=array_int, 
         floats=array_float, 
         matrix=array_2d)
print(f"Multiple arrays saved to: {npz_filename}")

# Load multiple arrays
loaded_arrays = np.load(npz_filename)
print("Loaded arrays:")
for key in loaded_arrays.files:
    print(f"  {key}: shape={loaded_arrays[key].shape}, dtype={loaded_arrays[key].dtype}")
print()

# Compressed save
with tempfile.NamedTemporaryFile(suffix='.npz', delete=False) as f:
    compressed_filename = f.name

np.savez_compressed(compressed_filename,
                   integers=array_int,
                   floats=array_float,
                   matrix=array_2d)

# Compare file sizes
original_size = os.path.getsize(npz_filename)
compressed_size = os.path.getsize(compressed_filename)
compression_ratio = original_size / compressed_size

print(f"File size comparison:")
print(f"  Original file: {original_size} bytes")
print(f"  Compressed file: {compressed_size} bytes")
print(f"  Compression ratio: {compression_ratio:.2f}")
print()

# Cleanup temporary files
os.unlink(npy_filename)
os.unlink(npz_filename)
os.unlink(compressed_filename)

import numpy as np
import tempfile
import os

# Create large array (simulating big data)
large_array = np.random.random((1000, 1000))

print("Memory-mapped files:")
print(f"Large array shape: {large_array.shape}")
print(f"Array size: {large_array.nbytes / 1024 / 1024:.2f} MB")
print()

# Save as memory-mapped file
with tempfile.NamedTemporaryFile(suffix='.dat', delete=False) as f:
    memmap_filename = f.name

# Create memory-mapped file
memmap_array = np.memmap(memmap_filename, dtype='float64', mode='w+', shape=(1000, 1000))
memmap_array[:] = large_array[:]
memmap_array.flush()  # Ensure data is written to disk

print(f"Memory-mapped file created: {memmap_filename}")
print(f"File size: {os.path.getsize(memmap_filename) / 1024 / 1024:.2f} MB")
print()

# Open memory-mapped file in read-only mode
readonly_memmap = np.memmap(memmap_filename, dtype='float64', mode='r', shape=(1000, 1000))

# Verify data
print("Data verification:")
print(f"Original array mean: {np.mean(large_array):.6f}")
print(f"Memory-mapped mean: {np.mean(readonly_memmap):.6f}")
print(f"Data matches: {np.allclose(large_array, readonly_memmap)}")
print()

# Partial read (advantage of memory mapping)
print("Partial data reading:")
subset = readonly_memmap[100:200, 100:200]  # Only read small portion
print(f"Subset shape: {subset.shape}")
print(f"Subset mean: {np.mean(subset):.6f}")
print()

# Cleanup
del memmap_array, readonly_memmap  # Close memory maps
os.unlink(memmap_filename)

import numpy as np
import sys

print("Memory usage optimization:")

# Impact of data types on memory
size = 1000000

array_float64 = np.ones(size, dtype=np.float64)
array_float32 = np.ones(size, dtype=np.float32)
array_int32 = np.ones(size, dtype=np.int32)
array_int16 = np.ones(size, dtype=np.int16)
array_int8 = np.ones(size, dtype=np.int8)

print(f"Array size: {size:,} elements")
print(f"float64: {array_float64.nbytes / 1024 / 1024:.2f} MB")
print(f"float32: {array_float32.nbytes / 1024 / 1024:.2f} MB")
print(f"int32:   {array_int32.nbytes / 1024 / 1024:.2f} MB")
print(f"int16:   {array_int16.nbytes / 1024 / 1024:.2f} MB")
print(f"int8:    {array_int8.nbytes / 1024 / 1024:.2f} MB")
print()

# View vs Copy
original = np.arange(1000000)
print(f"Original array memory: {original.nbytes / 1024 / 1024:.2f} MB")

# Create view (no extra memory)
view = original[::2]  # Every other element
print(f"View shares memory: {np.shares_memory(original, view)}")
print(f"View base object: {view.base is original}")

# Create copy (extra memory)
copy = original[::2].copy()
print(f"Copy shares memory: {np.shares_memory(original, copy)}")
print(f"Copy base object: {copy.base is None}")
print()

# In-place operations vs creating new arrays
test_array = np.random.random(100000)
original_id = id(test_array)

# In-place operation
test_array += 1
print(f"Same ID after in-place operation: {id(test_array) == original_id}")

# Creating new array operation
test_array = test_array + 1
print(f"Same ID after new array operation: {id(test_array) == original_id}")

import numpy as np
import time

def time_function(func, *args, **kwargs):
    """Measure function execution time"""
    start_time = time.time()
    result = func(*args, **kwargs)
    end_time = time.time()
    return result, end_time - start_time

print("Performance testing and optimization:")

# Create test data
size = 1000000
array1 = np.random.random(size)
array2 = np.random.random(size)

print(f"Test array size: {size:,}")
print()

# Test different calculation methods
print("Calculate sqrt(x² + y²) using different methods:")

# Method 1: Using np.sqrt
result1, time1 = time_function(lambda: np.sqrt(array1**2 + array2**2))
print(f"Method 1 (np.sqrt): {time1:.4f} seconds")

# Method 2: Using np.hypot
result2, time2 = time_function(lambda: np.hypot(array1, array2))
print(f"Method 2 (np.hypot): {time2:.4f} seconds")

# Method 3: Using np.linalg.norm
stacked = np.stack([array1, array2], axis=0)
result3, time3 = time_function(lambda: np.linalg.norm(stacked, axis=0))
print(f"Method 3 (np.linalg.norm): {time3:.4f} seconds")

print(f"Results equal: {np.allclose(result1, result2) and np.allclose(result2, result3)}")
print(f"Fastest method: Method {np.argmin([time1, time2, time3]) + 1}")
print()

# Memory layout impact on performance
print("Memory layout impact on performance:")

# C-style (row-major) vs Fortran-style (column-major)
array_c = np.random.random((1000, 1000))
array_f = np.asfortranarray(array_c)

print(f"C-style array: {array_c.flags['C_CONTIGUOUS']}")
print(f"Fortran-style array: {array_f.flags['F_CONTIGUOUS']}")

# Row sum
_, time_c_row = time_function(lambda: np.sum(array_c, axis=1))
_, time_f_row = time_function(lambda: np.sum(array_f, axis=1))

# Column sum
_, time_c_col = time_function(lambda: np.sum(array_c, axis=0))
_, time_f_col = time_function(lambda: np.sum(array_f, axis=0))

print(f"C-style row sum: {time_c_row:.4f} seconds")
print(f"C-style column sum: {time_c_col:.4f} seconds")
print(f"F-style row sum: {time_f_row:.4f} seconds")
print(f"F-style column sum: {time_f_col:.4f} seconds")

import numpy as np

def design_fir_filter(cutoff_freq, sampling_freq, num_taps):
    """Design FIR low-pass filter"""
    # Normalized cutoff frequency
    normalized_cutoff = cutoff_freq / (sampling_freq / 2)
    
    # Generate ideal low-pass filter impulse response
    n = np.arange(num_taps)
    h = np.sinc(2 * normalized_cutoff * (n - (num_taps - 1) / 2))
    
    # Apply Hamming window
    window = np.hamming(num_taps)
    h = h * window
    
    # Normalize
    h = h / np.sum(h)
    
    return h

def apply_filter(signal, filter_coeffs):
    """Apply FIR filter"""
    return np.convolve(signal, filter_coeffs, mode='same')

# Generate test signal
fs = 1000  # Sampling frequency
t = np.linspace(0, 1, fs, endpoint=False)

# Combined signal: 50Hz + 150Hz + noise
signal = (np.sin(2 * np.pi * 50 * t) + 
         0.5 * np.sin(2 * np.pi * 150 * t) + 
         0.2 * np.random.normal(0, 1, len(t)))

print("Digital Signal Processing Example:")
print(f"Sampling frequency: {fs} Hz")
print(f"Signal length: {len(signal)} samples")
print(f"Signal contains: 50Hz + 150Hz + noise")
print()

# Design low-pass filter (cutoff 100Hz)
cutoff = 100
num_taps = 101
filter_coeffs = design_fir_filter(cutoff, fs, num_taps)

print(f"FIR filter design:")
print(f"  Cutoff frequency: {cutoff} Hz")
print(f"  Filter order: {num_taps}")
print(f"  Filter coefficient range: [{np.min(filter_coeffs):.6f}, {np.max(filter_coeffs):.6f}]")
print()

# Apply filter
filtered_signal = apply_filter(signal, filter_coeffs)

# Analyze filtering effect
original_power = np.mean(signal**2)
filtered_power = np.mean(filtered_signal**2)
print(f"Filtering effect:")
print(f"  Original signal power: {original_power:.4f}")
print(f"  Filtered power: {filtered_power:.4f}")
print(f"  Power retention: {filtered_power/original_power*100:.1f}%")
print()

# Frequency domain analysis
original_fft = np.fft.fft(signal)
filtered_fft = np.fft.fft(filtered_signal)
freqs = np.fft.fftfreq(len(signal), 1/fs)

# Calculate attenuation at 50Hz and 150Hz
freq_50_idx = np.argmin(np.abs(freqs - 50))
freq_150_idx = np.argmin(np.abs(freqs - 150))

attenuation_50 = 20 * np.log10(np.abs(filtered_fft[freq_50_idx]) / np.abs(original_fft[freq_50_idx]))
attenuation_150 = 20 * np.log10(np.abs(filtered_fft[freq_150_idx]) / np.abs(original_fft[freq_150_idx]))

print(f"Frequency attenuation:")
print(f"  Attenuation at 50Hz: {attenuation_50:.1f} dB")
print(f"  Attenuation at 150Hz: {attenuation_150:.1f} dB")

import numpy as np

def create_test_image(size=128):
    """Create test image"""
    x = np.linspace(-2, 2, size)
    y = np.linspace(-2, 2, size)
    X, Y = np.meshgrid(x, y)
    
    # Create complex pattern
    image = (np.sin(3 * X) * np.cos(3 * Y) + 
            0.5 * np.sin(8 * X) * np.sin(8 * Y) + 
            0.3 * np.random.random((size, size)))
    
    # Normalize to [0, 255]
    image = ((image - np.min(image)) / (np.max(image) - np.min(image)) * 255).astype(np.uint8)
    return image

def gaussian_kernel(size, sigma):
    """Generate Gaussian kernel"""
    kernel = np.zeros((size, size))
    center = size // 2
    
    for i in range(size):
        for j in range(size):
            x, y = i - center, j - center
            kernel[i, j] = np.exp(-(x**2 + y**2) / (2 * sigma**2))
    
    return kernel / np.sum(kernel)

def convolve_2d(image, kernel):
    """2D convolution"""
    pad_h, pad_w = kernel.shape[0] // 2, kernel.shape[1] // 2
    padded = np.pad(image, ((pad_h, pad_h), (pad_w, pad_w)), mode='reflect')
    
    result = np.zeros_like(image, dtype=np.float32)
    
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            result[i, j] = np.sum(padded[i:i+kernel.shape[0], j:j+kernel.shape[1]] * kernel)
    
    return result.astype(np.uint8)

def edge_detection(image):
    """Edge detection"""
    # Sobel operators
    sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    
    # Calculate gradients
    grad_x = convolve_2d(image, sobel_x)
    grad_y = convolve_2d(image, sobel_y)
    
    # Calculate gradient magnitude
    magnitude = np.sqrt(grad_x.astype(np.float32)**2 + grad_y.astype(np.float32)**2)
    magnitude = (magnitude / np.max(magnitude) * 255).astype(np.uint8)
    
    return magnitude

# Image processing pipeline
print("Advanced Image Processing Example:")

# Create test image
original_image = create_test_image(64)
print(f"Original image size: {original_image.shape}")
print(f"Pixel value range: [{np.min(original_image)}, {np.max(original_image)}]")
print()

# Gaussian blur
gaussian_5x5 = gaussian_kernel(5, 1.0)
blurred_image = convolve_2d(original_image, gaussian_5x5)

print(f"Gaussian blur:")
print(f"  Kernel size: 5x5")
print(f"  Sigma: 1.0")
print(f"  Blurred pixel range: [{np.min(blurred_image)}, {np.max(blurred_image)}]")
print()

# Edge detection
edges = edge_detection(original_image)
print(f"Edge detection:")
print(f"  Edge pixel range: [{np.min(edges)}, {np.max(edges)}]")
print(f"  Edge pixel count: {np.sum(edges > 50)}")
print()

# Image statistics
print(f"Image statistics comparison:")
print(f"  Original image mean: {np.mean(original_image):.1f}")
print(f"  Blurred image mean: {np.mean(blurred_image):.1f}")
print(f"  Edge image mean: {np.mean(edges):.1f}")
print()
print(f"  Original image std: {np.std(original_image):.1f}")
print(f"  Blurred image std: {np.std(blurred_image):.1f}")
print(f"  Edge image std: {np.std(edges):.1f}")