My Practical Insights of Using SciPy Library

import numpy as np import scipy from scipy import optimize, linalg, integrate, interpolate from scipy import stats, signal, sparse, spatial import matplotlib.pyplot as plt # Check SciPy version and available modules print(f"SciPy version: {scipy.__version__}") print(f"Available subpackages: {[pkg for pkg in dir(scipy) if not pkg.startswith('_')][:10]}") # Quick demonstration of different modules print("\n=== SCIPY MODULES DEMONSTRATION ===") # 1. Linear Algebra (enhanced NumPy capabilities) A = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]]) b = np.array([2, 4, -1]) x = linalg.solve(A, b) print(f"Linear system solution: {x}") # 2. Optimization - find minimum of a function def objective_function(x): return x**2 + 10*np.sin(x) result = optimize.minimize_scalar(objective_function) print(f"Minimum found at x = {result.x:.4f}, f(x) = {result.fun:.4f}") # 3. Integration - compute definite integral def integrand(x): return np.exp(-x**2) integral_result, error = integrate.quad(integrand, 0, 1) print(f"Integral result: {integral_result:.6f} ± {error:.2e}") # 4. Statistics - generate random samples and fit distribution data = stats.norm.rvs(loc=0, scale=1, size=1000) mu, sigma = stats.norm.fit(data) print(f"Fitted normal distribution: μ = {mu:.4f}, σ = {sigma:.4f}") # 5. Signal processing - create and filter a signal t = np.linspace(0, 1, 500) noisy_signal = np.sin(2*np.pi*5*t) + 0.3*np.random.randn(len(t)) b, a = signal.butter(4, 0.1, 'low') filtered_signal = signal.filtfilt(b, a, noisy_signal) print(f"Signal filtered: Original std = {noisy_signal.std():.3f}, Filtered std = {filtered_signal.std():.3f}")

# Advanced optimization techniques and practical applications # 1. Multi-dimensional optimization def rosenbrock(x): """The Rosenbrock function - a classic optimization test case""" return 100 * (x[1] - x[0]**2)**2 + (1 - x[0])**2 # Different optimization methods comparison methods = ['BFGS', 'L-BFGS-B', 'Powell', 'Nelder-Mead'] initial_guess = [0, 0] print("=== OPTIMIZATION METHODS COMPARISON ===") for method in methods: result = optimize.minimize(rosenbrock, initial_guess, method=method) print(f"{method:12}: x* = [{result.x[0]:.4f}, {result.x[1]:.4f}], " f"f(x*) = {result.fun:.6f}, iterations = {result.nit}") # 2. Constrained optimization def objective(x): return -(x[0] * x[1] * x[2]) # Maximize volume (minimize negative) def constraint1(x): return x[0] + 2*x[1] + 2*x[2] - 10 # Surface area constraint def constraint2(x): return x[0] - x[1] # Design constraint constraints = [{'type': 'ineq', 'fun': constraint1}, {'type': 'ineq', 'fun': constraint2}] bounds = [(0, None), (0, None), (0, None)] # All variables positive constrained_result = optimize.minimize(objective, [1, 1, 1], method='SLSQP', bounds=bounds, constraints=constraints) print(f"\n=== CONSTRAINED OPTIMIZATION ===") print(f"Optimal dimensions: {constrained_result.x}") print(f"Maximum volume: {-constrained_result.fun:.4f}") # 3. Root finding methods def equation_system(vars): x, y = vars return [x**2 + y**2 - 4, # Circle equation x - y - 1] # Line equation # Find intersection points roots = optimize.fsolve(equation_system, [1, 1]) print(f"\n=== ROOT FINDING ===") print(f"Intersection point: x = {roots[0]:.4f}, y = {roots[1]:.4f}") print(f"Verification: x² + y² = {roots[0]**2 + roots[1]**2:.4f} (should be 4)") # 4. Curve fitting and parameter estimation def exponential_decay(t, A, k, C): """Exponential decay model: A * exp(-k*t) + C""" return A * np.exp(-k * t) + C # Generate synthetic data with noise t_data = np.linspace(0, 5, 50) true_params = [10, 0.8, 2] y_true = exponential_decay(t_data, *true_params) noise = np.random.normal(0, 0.5, len(t_data)) y_data = y_true + noise # Fit the curve popt, pcov = optimize.curve_fit(exponential_decay, t_data, y_data) print(f"\n=== CURVE FITTING ===") print(f"True parameters: A={true_params[0]}, k={true_params[1]}, C={true_params[2]}") print(f"Fitted parameters: A={popt[0]:.2f}, k={popt[1]:.2f}, C={popt[2]:.2f}") print(f"Parameter uncertainties: {np.sqrt(np.diag(pcov))}") # 5. Global optimization for avoiding local minima def multi_modal_func(x): """Function with multiple local minima""" return np.sin(x[0]) * np.cos(x[1]) + 0.1 * (x[0]**2 + x[1]**2) # Local optimization (might get stuck) local_result = optimize.minimize(multi_modal_func, [1, 1], method='BFGS') # Global optimization bounds_global = [(-5, 5), (-5, 5)] global_result = optimize.differential_evolution(multi_modal_func, bounds_global) print(f"\n=== GLOBAL vs LOCAL OPTIMIZATION ===") print(f"Local minimum: x = [{local_result.x[0]:.4f}, {local_result.x[1]:.4f}], " f"f(x) = {local_result.fun:.6f}") print(f"Global minimum: x = [{global_result.x[0]:.4f}, {global_result.x[1]:.4f}], " f"f(x) = {global_result.fun:.6f}")

# Comprehensive statistical analysis and hypothesis testing # 1. Working with probability distributions print("=== PROBABILITY DISTRIBUTIONS ===") # Generate samples from different distributions normal_data = stats.norm.rvs(loc=0, scale=1, size=1000) exponential_data = stats.expon.rvs(scale=2, size=1000) uniform_data = stats.uniform.rvs(loc=-1, scale=2, size=1000) # Distribution fitting and goodness-of-fit tests distributions = [ (normal_data, 'Normal', stats.norm), (exponential_data, 'Exponential', stats.expon), (uniform_data, 'Uniform', stats.uniform) ] for data, name, dist in distributions: # Fit distribution parameters params = dist.fit(data) # Kolmogorov-Smirnov test ks_stat, ks_p = stats.kstest(data, lambda x: dist.cdf(x, *params)) print(f"{name:12}: KS test p-value = {ks_p:.4f} {'✓' if ks_p > 0.05 else '✗'}") # 2. Hypothesis testing print(f"\n=== HYPOTHESIS TESTING ===") # Generate two groups for comparison group1 = stats.norm.rvs(loc=100, scale=15, size=50) group2 = stats.norm.rvs(loc=105, scale=15, size=50) # Independent t-test t_stat, t_p = stats.ttest_ind(group1, group2) print(f"T-test: t = {t_stat:.4f}, p = {t_p:.4f}") # Mann-Whitney U test (non-parametric) u_stat, u_p = stats.mannwhitneyu(group1, group2) print(f"Mann-Whitney U test: U = {u_stat:.0f}, p = {u_p:.4f}") # Chi-square test for independence observed = np.array([[10, 10, 20], [20, 20, 40]]) chi2_stat, chi2_p, dof, expected = stats.chi2_contingency(observed) print(f"Chi-square test: χ² = {chi2_stat:.4f}, p = {chi2_p:.4f}") # 3. Correlation analysis x = np.random.randn(100) y = 2*x + np.random.randn(100) # Pearson correlation pearson_r, pearson_p = stats.pearsonr(x, y) print(f"\n=== CORRELATION ANALYSIS ===") print(f"Pearson correlation: r = {pearson_r:.4f}, p = {pearson_p:.4f}") # Spearman rank correlation spearman_r, spearman_p = stats.spearmanr(x, y) print(f"Spearman correlation: ρ = {spearman_r:.4f}, p = {spearman_p:.4f}") # 4. Bootstrap confidence intervals def bootstrap_mean(data, n_bootstrap=1000, confidence=0.95): """Calculate bootstrap confidence interval for the mean""" bootstrap_means = [] n = len(data) for _ in range(n_bootstrap): bootstrap_sample = np.random.choice(data, size=n, replace=True) bootstrap_means.append(np.mean(bootstrap_sample)) alpha = 1 - confidence lower = np.percentile(bootstrap_means, 100 * alpha/2) upper = np.percentile(bootstrap_means, 100 * (1 - alpha/2)) return np.mean(bootstrap_means), lower, upper sample_data = stats.norm.rvs(loc=50, scale=10, size=30) boot_mean, boot_lower, boot_upper = bootstrap_mean(sample_data) print(f"\n=== BOOTSTRAP CONFIDENCE INTERVALS ===") print(f"Sample mean: {np.mean(sample_data):.2f}") print(f"Bootstrap 95% CI: [{boot_lower:.2f}, {boot_upper:.2f}]") # 5. Multiple testing correction np.random.seed(42) n_tests = 20 p_values = [] for i in range(n_tests): # Generate random data (most should show no effect) if i < 3: # First 3 tests have real effects group_a = stats.norm.rvs(0, 1, 30) group_b = stats.norm.rvs(0.8, 1, 30) # Effect size = 0.8 else: group_a = stats.norm.rvs(0, 1, 30) group_b = stats.norm.rvs(0, 1, 30) # No effect _, p = stats.ttest_ind(group_a, group_b) p_values.append(p) # Benjamini-Hochberg correction rejected, corrected_p, _, _ = stats.multipletests(p_values, method='fdr_bh') print(f"\n=== MULTIPLE TESTING CORRECTION ===") print(f"Significant tests (uncorrected): {sum(np.array(p_values) < 0.05)}") print(f"Significant tests (FDR corrected): {sum(rejected)}") print(f"True positives detected: {sum(rejected[:3])}/3")

# Advanced linear algebra operations and matrix decompositions # 1. Matrix decompositions A = np.random.randn(5, 5) A = A @ A.T # Make positive definite print("=== MATRIX DECOMPOSITIONS ===") # Cholesky decomposition (for positive definite matrices) try: L = linalg.cholesky(A, lower=True) print(f"Cholesky decomposition successful: A = L @ L.T") print(f"Reconstruction error: {np.linalg.norm(A - L @ L.T):.2e}") except linalg.LinAlgError: print("Matrix is not positive definite") # QR decomposition B = np.random.randn(6, 4) # Overdetermined system Q, R = linalg.qr(B) print(f"QR decomposition: B({B.shape}) = Q({Q.shape}) @ R({R.shape})") print(f"Q orthogonality check: ||Q.T @ Q - I|| = {np.linalg.norm(Q.T @ Q - np.eye(Q.shape[1])):.2e}") # Singular Value Decomposition (SVD) U, s, Vt = linalg.svd(B) print(f"SVD: B = U({U.shape}) @ diag(s) @ V.T({Vt.shape})") print(f"Condition number: {s[0]/s[-1]:.2f}") # 2. Solving different types of linear systems print(f"\n=== SOLVING LINEAR SYSTEMS ===") # Standard system Ax = b A_sys = np.random.randn(4, 4) b_sys = np.random.randn(4) x_solve = linalg.solve(A_sys, b_sys) print(f"Standard solve residual: {np.linalg.norm(A_sys @ x_solve - b_sys):.2e}") # Triangular system (more efficient) T = np.triu(A_sys) # Upper triangular x_triangular = linalg.solve_triangular(T, b_sys) print(f"Triangular solve residual: {np.linalg.norm(T @ x_triangular - b_sys):.2e}") # Least squares for overdetermined systems A_over = np.random.randn(10, 4) b_over = np.random.randn(10) x_lstsq, residuals, rank, s = linalg.lstsq(A_over, b_over) print(f"Least squares: residual norm = {residuals[0]:.3f}, rank = {rank}") # 3. Eigenvalue problems print(f"\n=== EIGENVALUE PROBLEMS ===") # Standard eigenvalue problem eigenvals, eigenvecs = linalg.eigh(A) # For symmetric matrices print(f"Eigenvalues: {eigenvals}") print(f"Largest eigenvalue: {eigenvals[-1]:.4f}") # Generalized eigenvalue problem Ax = λBx B_gen = np.random.randn(5, 5) B_gen = B_gen @ B_gen.T + np.eye(5) # Make positive definite gen_eigenvals, gen_eigenvecs = linalg.eigh(A, B_gen) print(f"Generalized eigenvalues: {gen_eigenvals}") # 4. Matrix functions and special operations print(f"\n=== MATRIX FUNCTIONS ===") # Matrix exponential A_small = np.array([[0, 1], [-1, 0]]) # Rotation matrix generator exp_A = linalg.expm(A_small) print(f"Matrix exponential computed") print(f"exp(A) properties: det = {linalg.det(exp_A):.4f} (should be ~1)") # Matrix square root sqrt_A = linalg.sqrtm(A) print(f"Matrix square root error: {np.linalg.norm(sqrt_A @ sqrt_A - A):.2e}") # Matrix logarithm log_A = linalg.logm(exp_A) print(f"Matrix log(exp(A)) error: {np.linalg.norm(log_A - A_small):.2e}") # 5. Sparse matrix operations print(f"\n=== SPARSE MATRIX OPERATIONS ===") # Create a sparse matrix from scipy.sparse import csr_matrix, linalg as sparse_linalg # Random sparse matrix (10% non-zero elements) n = 1000 density = 0.1 A_sparse_dense = np.random.randn(n, n) mask = np.random.random((n, n)) > density A_sparse_dense[mask] = 0 A_sparse = csr_matrix(A_sparse_dense) print(f"Sparse matrix: {A_sparse.shape}, {A_sparse.nnz} non-zeros ({A_sparse.nnz/n**2*100:.1f}%)") # Sparse eigenvalues (only compute a few) eigenvals_sparse, eigenvecs_sparse = sparse_linalg.eigs(A_sparse, k=6, which='LM') print(f"Largest magnitude eigenvalues: {np.abs(eigenvals_sparse)}") # Memory comparison dense_memory = A_sparse_dense.nbytes / 1024**2 sparse_memory = (A_sparse.data.nbytes + A_sparse.indices.nbytes + A_sparse.indptr.nbytes) / 1024**2 print(f"Memory usage: Dense = {dense_memory:.1f}MB, Sparse = {sparse_memory:.1f}MB") print(f"Memory savings: {(1 - sparse_memory/dense_memory)*100:.1f}%")

# Numerical integration and solving differential equations # 1. Numerical integration of functions print("=== NUMERICAL INTEGRATION ===") # Simple integration def simple_func(x): return np.exp(-x**2) * np.cos(x) result, error = integrate.quad(simple_func, 0, 2) print(f"Single integral: ∫₀² exp(-x²)cos(x) dx = {result:.6f} ± {error:.2e}") # Double integration def double_integrand(y, x): return np.exp(-(x**2 + y**2)) double_result, double_error = integrate.dblquad(double_integrand, 0, 1, 0, 1) print(f"Double integral: ∫∫ exp(-(x²+y²)) dxdy = {double_result:.6f} ± {double_error:.2e}") # Integration with parameters def parametric_func(x, a, b): return a * x**2 + b * np.sin(x) param_result, _ = integrate.quad(parametric_func, 0, np.pi, args=(2, 3)) print(f"Parametric integral result: {param_result:.6f}") # 2. Solving Ordinary Differential Equations (ODEs) print(f"\n=== ORDINARY DIFFERENTIAL EQUATIONS ===") # Simple ODE: dy/dt = -k*y (exponential decay) def exponential_decay(t, y, k): return -k * y # Initial conditions t_span = (0, 5) y0 = [10.0] k = 0.5 # Solve using different methods methods = ['RK45', 'RK23', 'BDF'] for method in methods: sol = integrate.solve_ivp(exponential_decay, t_span, y0, args=(k,), method=method, dense_output=True) # Evaluate at specific points t_eval = np.linspace(0, 5, 11) y_eval = sol.sol(t_eval) # Compare with analytical solution y_analytical = y0[0] * np.exp(-k * t_eval) error = np.mean(np.abs(y_eval[0] - y_analytical)) print(f"{method:5}: Final value = {y_eval[0][-1]:.4f}, " f"Mean error = {error:.2e}, Steps = {len(sol.t)}") # 3. System of ODEs: Predator-Prey model (Lotka-Volterra) def lotka_volterra(t, z, alpha, beta, gamma, delta): """ Lotka-Volterra predator-prey model z[0] = prey population, z[1] = predator population """ x, y = z dxdt = alpha * x - beta * x * y dydt = delta * x * y - gamma * y return [dxdt, dydt] # Parameters and initial conditions alpha, beta, gamma, delta = 1.0, 0.5, 0.75, 0.25 initial_populations = [10, 5] # [prey, predator] t_span_lv = (0, 15) # Solve the system lv_solution = integrate.solve_ivp(lotka_volterra, t_span_lv, initial_populations, args=(alpha, beta, gamma, delta), dense_output=True, rtol=1e-8) t_plot = np.linspace(0, 15, 1000) populations = lv_solution.sol(t_plot) print(f"\n=== LOTKA-VOLTERRA SYSTEM ===") print(f"Simulation time span: {t_span_lv}") print(f"Final populations: Prey = {populations[0][-1]:.2f}, Predator = {populations[1][-1]:.2f}") print(f"Population oscillation period: ~{2*np.pi/np.sqrt(alpha*gamma):.2f} time units") # 4. Boundary Value Problems (BVP) def boundary_value_problem(x, y): """ Second-order BVP: y'' + y = 0 Converted to system: y₁ = y, y₂ = y' """ return np.vstack((y[1], -y[0])) def boundary_conditions(ya, yb): """ Boundary conditions: y(0) = 0, y(π) = 1 """ return np.array([ya[0], yb[0] - 1]) # Initial guess for the solution x_bvp = np.linspace(0, np.pi, 11) y_guess = np.zeros((2, x_bvp.size)) y_guess[0] = np.sin(x_bvp) # Guess for y y_guess[1] = np.cos(x_bvp) # Guess for y' # Solve BVP bvp_solution = integrate.solve_bvp(boundary_value_problem, boundary_conditions, x_bvp, y_guess) print(f"\n=== BOUNDARY VALUE PROBLEM ===") print(f"BVP solved successfully: {bvp_solution.success}") print(f"Solution at x=π/2: y = {bvp_solution.sol(np.pi/2)[0]:.6f}") print(f"Analytical value: {np.sin(np.pi/2):.6f}") # 5. Advanced integration techniques print(f"\n=== ADVANCED INTEGRATION TECHNIQUES ===") # Oscillatory integrands def oscillatory(x): return np.sin(50*x) * np.exp(-x) # Standard quad might struggle with this osc_result, osc_error = integrate.quad(oscillatory, 0, 10, limit=100) print(f"Oscillatory integral: {osc_result:.6f} ± {osc_error:.2e}") # Infinite limits def decay_func(x): return np.exp(-x) / (1 + x**2) inf_result, inf_error = integrate.quad(decay_func, 0, np.inf) print(f"Infinite integral: {inf_result:.6f} ± {inf_error:.2e}") # Sample-based integration (Monte Carlo) def monte_carlo_integration(func, a, b, n_samples=100000): """Simple Monte Carlo integration""" x_samples = np.random.uniform(a, b, n_samples) y_samples = func(x_samples) return (b - a) * np.mean(y_samples) mc_result = monte_carlo_integration(lambda x: x**2, 0, 1, 100000) analytical_mc = 1/3 print(f"Monte Carlo ∫₀¹ x² dx: {mc_result:.6f} (analytical: {analytical_mc:.6f})") print(f"MC Error: {abs(mc_result - analytical_mc):.6f}")

# Digital signal processing and frequency analysis # 1. Signal generation and basic operations print("=== SIGNAL GENERATION ===") # Create a composite signal fs = 1000 # Sampling frequency t = np.linspace(0, 2, fs*2, endpoint=False) # Multiple frequency components + noise f1, f2, f3 = 50, 120, 250 # Hz signal = (2*np.sin(2*np.pi*f1*t) + 1.5*np.sin(2*np.pi*f2*t) + np.sin(2*np.pi*f3*t) + 0.3*np.random.randn(len(t))) print(f"Generated signal: {len(t)} samples at {fs} Hz") print(f"Signal statistics: mean = {np.mean(signal):.3f}, std = {np.std(signal):.3f}") # 2. Filter design and application print(f"\n=== DIGITAL FILTERING ===") # Design different types of filters # Low-pass filter nyquist = 0.5 * fs cutoff_low = 100 # Hz normalized_cutoff = cutoff_low / nyquist b_low, a_low = signal.butter(4, normalized_cutoff, btype='low') filtered_low = signal.filtfilt(b_low, a_low, signal) # High-pass filter cutoff_high = 80 normalized_cutoff_high = cutoff_high / nyquist b_high, a_high = signal.butter(4, normalized_cutoff_high, btype='high') filtered_high = signal.filtfilt(b_high, a_high, signal) # Band-pass filter low_freq, high_freq = 40, 200 normalized_band = [low_freq/nyquist, high_freq/nyquist] b_band, a_band = signal.butter(4, normalized_band, btype='band') filtered_band = signal.filtfilt(b_band, a_band, signal) print(f"Filter results (RMS values):") print(f"Original signal: {np.sqrt(np.mean(signal**2)):.3f}") print(f"Low-pass filtered: {np.sqrt(np.mean(filtered_low**2)):.3f}") print(f"High-pass filtered: {np.sqrt(np.mean(filtered_high**2)):.3f}") print(f"Band-pass filtered: {np.sqrt(np.mean(filtered_band**2)):.3f}") # 3. Spectral analysis print(f"\n=== SPECTRAL ANALYSIS ===") # Power Spectral Density (PSD) frequencies, psd = signal.welch(signal, fs=fs, nperseg=1024) # Find peaks in the spectrum peaks, properties = signal.find_peaks(psd, height=0.01, distance=20) peak_frequencies = frequencies[peaks] peak_powers = psd[peaks] print(f"Detected frequency peaks:") for f_peak, p_peak in zip(peak_frequencies, peak_powers): print(f" {f_peak:.1f} Hz: {10*np.log10(p_peak):.1f} dB") # Spectrogram for time-frequency analysis f_spec, t_spec, Sxx = signal.spectrogram(signal, fs=fs, nperseg=256, noverlap=128) print(f"Spectrogram shape: {Sxx.shape} (freq × time)") # 4. Cross-correlation and convolution print(f"\n=== CORRELATION AND CONVOLUTION ===") # Create a template signal (short pulse) template_t = np.linspace(0, 0.1, int(0.1*fs)) template = np.sin(2*np.pi*f1*template_t) * np.exp(-template_t*20) # Cross-correlation to find template in signal correlation = signal.correlate(signal, template, mode='valid') correlation_lags = signal.correlation_lags(len(signal), len(template), mode='valid') # Find the best match max_corr_idx = np.argmax(np.abs(correlation)) best_match_time = correlation_lags[max_corr_idx] / fs print(f"Template matching:") print(f"Template length: {len(template)} samples ({len(template)/fs:.3f} s)") print(f"Best match at: {best_match_time:.3f} s") print(f"Correlation coefficient: {np.max(np.abs(correlation)):.3f}") # 5. System identification and transfer functions print(f"\n=== SYSTEM ANALYSIS ===") # Define a system (second-order low-pass filter) # Transfer function: H(s) = ωₙ²/(s² + 2ζωₙs + ωₙ²) wn = 2*np.pi*50 # Natural frequency zeta = 0.7 # Damping ratio # Create transfer function num = [wn**2] den = [1, 2*zeta*wn, wn**2] system = signal.TransferFunction(num, den) # System response analysis w, h = signal.freqresp(system) frequencies_hz = w / (2*np.pi) # Find -3dB bandwidth magnitude_db = 20 * np.log10(np.abs(h)) idx_3db = np.where(magnitude_db <= -3)[0] if len(idx_3db) > 0: bandwidth = frequencies_hz[idx_3db[0]] print(f"System -3dB bandwidth: {bandwidth:.1f} Hz") # Step response t_step, y_step = signal.step(system) settling_time_idx = np.where(np.abs(y_step - y_step[-1]) < 0.02*y_step[-1])[0] if len(settling_time_idx) > 0: settling_time = t_step[settling_time_idx[0]] print(f"Step response settling time (2%): {settling_time:.3f} s") # Impulse response t_impulse, y_impulse = signal.impulse(system) print(f"Impulse response calculated: {len(t_impulse)} points") # 6. Advanced signal processing techniques print(f"\n=== ADVANCED TECHNIQUES ===") # Hilbert transform for analytic signal analytic_signal = signal.hilbert(filtered_band) amplitude_envelope = np.abs(analytic_signal) instantaneous_phase = np.unwrap(np.angle(analytic_signal)) instantaneous_frequency = np.diff(instantaneous_phase) / (2*np.pi) * fs print(f"Hilbert transform analysis:") print(f"Mean amplitude envelope: {np.mean(amplitude_envelope):.3f}") print(f"Mean instantaneous frequency: {np.mean(instantaneous_frequency):.1f} Hz") # Savitzky-Golay filter for smoothing smoothed_signal = signal.savgol_filter(signal, window_length=51, polyorder=3) smoothing_factor = np.var(signal - smoothed_signal) / np.var(signal) print(f"Savitzky-Golay smoothing factor: {smoothing_factor:.3f}") # Detrending detrended_signal = signal.detrend(signal) print(f"Signal detrending: removed {np.mean(signal):.3f} DC offset") # Peak detection with advanced parameters peaks_advanced, props_advanced = signal.find_peaks( np.abs(signal), height=1.0, # Minimum peak height distance=100, # Minimum peak separation (samples) width=10, # Minimum peak width prominence=0.5 # Peak prominence ) print(f"Advanced peak detection: found {len(peaks_advanced)} significant peaks") if len(peaks_advanced) > 0: peak_times = t[peaks_advanced] print(f"Peak times: {peak_times[:5]} ...") # Show first 5 peaks

# Production-ready scientific computing with SciPy import time import warnings from functools import wraps # 1. Performance monitoring and optimization def performance_monitor(func): """Decorator to monitor function performance""" @wraps(func) def wrapper(*args, **kwargs): start_time = time.time() start_memory = None try: import psutil process = psutil.Process() start_memory = process.memory_info().rss / 1024**2 except ImportError: pass result = func(*args, **kwargs) end_time = time.time() execution_time = end_time - start_time end_memory = None if start_memory is not None: end_memory = process.memory_info().rss / 1024**2 memory_delta = end_memory - start_memory print(f"Function: {func.__name__}") print(f"Execution time: {execution_time:.4f}s") if end_memory is not None: print(f"Memory delta: {memory_delta:+.1f}MB") return result return wrapper # 2. Robust optimization with error handling @performance_monitor def robust_optimization(objective_func, initial_guess, method='BFGS', **kwargs): """Robust optimization with fallback methods""" methods_to_try = [method, 'L-BFGS-B', 'Powell', 'Nelder-Mead'] for i, opt_method in enumerate(methods_to_try): try: with warnings.catch_warnings(): warnings.simplefilter("ignore") result = optimize.minimize( objective_func, initial_guess, method=opt_method, **kwargs ) if result.success: print(f"Optimization successful with {opt_method}") return result else: print(f"Method {opt_method} converged but reports failure") if i == len(methods_to_try) - 1: return result except Exception as e: print(f"Method {opt_method} failed: {e}") if i == len(methods_to_try) - 1: raise return None # Test robust optimization def challenging_function(x): """A function that might cause numerical issues""" return np.sum(x**4) + 0.1 * np.sum(x**2) + np.prod(np.abs(x)) result = robust_optimization(challenging_function, np.array([1.0, -0.5, 2.0])) print(f"Optimized solution: {result.x}") # 3. Numerical stability and conditioning @performance_monitor def analyze_matrix_conditioning(A, tolerance=1e-12): """Analyze matrix conditioning and suggest solutions""" condition_number = linalg.cond(A) print(f"\n=== MATRIX CONDITIONING ANALYSIS ===") print(f"Condition number: {condition_number:.2e}") if condition_number > 1e12: print("⚠️ Matrix is ill-conditioned") # Suggest regularization regularization_param = tolerance * np.trace(A) / A.shape[0] A_regularized = A + regularization_param * np.eye(A.shape[0]) cond_regularized = linalg.cond(A_regularized) print(f"Suggested regularization parameter: {regularization_param:.2e}") print(f"Regularized condition number: {cond_regularized:.2e}") return A_regularized, True elif condition_number > 1e8: print("⚠️ Matrix is moderately ill-conditioned") print("Consider using higher precision or iterative methods") else: print("✓ Matrix is well-conditioned") return A, False # Test with different matrix types print("Testing well-conditioned matrix:") A_good = np.random.randn(5, 5) A_good = A_good @ A_good.T + np.eye(5) analyze_matrix_conditioning(A_good) print("\nTesting ill-conditioned matrix:") A_bad = linalg.hilbert(8) # Hilbert matrix is notoriously ill-conditioned analyze_matrix_conditioning(A_bad) # 4. Efficient batch processing @performance_monitor def batch_statistical_analysis(data_batches, confidence_level=0.95): """Efficiently process multiple datasets""" results = [] n_batches = len(data_batches) print(f"\n=== BATCH STATISTICAL ANALYSIS ===") print(f"Processing {n_batches} datasets...") for i, data in enumerate(data_batches): try: # Basic statistics mean_val = np.mean(data) std_val = np.std(data, ddof=1) # Confidence interval for the mean n = len(data) se = std_val / np.sqrt(n) alpha = 1 - confidence_level t_critical = stats.t.ppf(1 - alpha/2, df=n-1) ci_lower = mean_val - t_critical * se ci_upper = mean_val + t_critical * se # Normality test _, normality_p = stats.shapiro(data[:5000] if len(data) > 5000 else data) # Outlier detection using IQR method q1, q3 = np.percentile(data, [25, 75]) iqr = q3 - q1 outlier_bounds = [q1 - 1.5*iqr, q3 + 1.5*iqr] outliers = np.sum((data < outlier_bounds[0]) | (data > outlier_bounds[1])) results.append({ 'dataset': i+1, 'n_samples': n, 'mean': mean_val, 'std': std_val, 'ci_lower': ci_lower, 'ci_upper': ci_upper, 'normality_p': normality_p, 'n_outliers': outliers, 'outlier_pct': 100 * outliers / n }) except Exception as e: print(f"Error processing dataset {i+1}: {e}") results.append({'dataset': i+1, 'error': str(e)}) return results # Test batch processing test_datasets = [ stats.norm.rvs(0, 1, 1000), stats.exponential.rvs(scale=2, size=1500), stats.uniform.rvs(-1, 2, 800), np.concatenate([stats.norm.rvs(0, 1, 900), stats.norm.rvs(5, 1, 100)]) # Bimodal ] batch_results = batch_statistical_analysis(test_datasets) print("Batch analysis results:") for result in batch_results: if 'error' not in result: print(f"Dataset {result['dataset']}: " f"mean = {result['mean']:.3f} ± {result['std']:.3f}, " f"outliers = {result['outlier_pct']:.1f}%") # 5. Memory-efficient large-scale computations def memory_efficient_correlation_matrix(data, chunk_size=1000): """Compute correlation matrix for large datasets using chunking""" n_features = data.shape[1] correlation_matrix = np.zeros((n_features, n_features)) print(f"\n=== MEMORY-EFFICIENT CORRELATION ===") print(f"Computing {n_features}×{n_features} correlation matrix") print(f"Data shape: {data.shape}, Chunk size: {chunk_size}") # Standardize the entire dataset first (needed for correlation) data_standardized = stats.zscore(data, axis=0) # Compute correlation in chunks for i in range(0, n_features, chunk_size): end_i = min(i + chunk_size, n_features) for j in range(i, n_features, chunk_size): end_j = min(j + chunk_size, n_features) # Compute correlation for this chunk chunk_corr = np.corrcoef(data_standardized[:, i:end_i].T, data_standardized[:, j:end_j].T) # Extract the relevant part and store if i == j: # Diagonal block correlation_matrix[i:end_i, j:end_j] = chunk_corr[:end_i-i, :end_j-j] else: # Off-diagonal blocks (symmetric) block = chunk_corr[:end_i-i, end_i-i:] correlation_matrix[i:end_i, j:end_j] = block correlation_matrix[j:end_j, i:end_i] = block.T return correlation_matrix # Test with moderate-sized data test_data_large = np.random.randn(5000, 50) corr_matrix = memory_efficient_correlation_matrix(test_data_large, chunk_size=20) print(f"Correlation matrix computed: {corr_matrix.shape}") print(f"Diagonal elements (should be 1): {np.diag(corr_matrix)[:5]}") # 6. Error handling and validation class ScientificComputingValidator: """Validation utilities for scientific computing""" @staticmethod def validate_optimization_result(result): """Validate optimization results""" issues = [] if not result.success: issues.append(f"Optimization failed: {result.message}") if hasattr(result, 'hess_inv') and result.hess_inv is not None: try: eigenvals = linalg.eigvals(result.hess_inv) if np.any(eigenvals <= 0): issues.append("Hessian is not positive definite (not at minimum)") except: pass if hasattr(result, 'jac') and result.jac is not None: gradient_norm = np.linalg.norm(result.jac) if gradient_norm > 1e-3: issues.append(f"Large gradient norm: {gradient_norm:.2e}") return issues @staticmethod def validate_statistical_test(statistic, p_value, alpha=0.05): """Validate statistical test results""" results = { 'statistic': statistic, 'p_value': p_value, 'significant': p_value < alpha, 'effect_size': 'unknown' } warnings = [] if p_value < 1e-16: warnings.append("P-value extremely small, check for numerical issues") if not 0 <= p_value <= 1: warnings.append(f"Invalid p-value: {p_value}") results['warnings'] = warnings return results # Test validation validator = ScientificComputingValidator() # Test optimization validation test_result = optimize.minimize(lambda x: x[0]**2 + x[1]**2, [1, 1]) opt_issues = validator.validate_optimization_result(test_result) print(f"\n=== VALIDATION RESULTS ===") print(f"Optimization issues: {opt_issues if opt_issues else 'None'}") # Test statistical validation sample1 = np.random.normal(0, 1, 100) sample2 = np.random.normal(0.5, 1, 100) t_stat, p_val = stats.ttest_ind(sample1, sample2) stat_results = validator.validate_statistical_test(t_stat, p_val) print(f"Statistical test: significant = {stat_results['significant']}, " f"warnings = {stat_results['warnings'] if stat_results['warnings'] else 'None'}") print(f"\n=== PRODUCTION TIPS SUMMARY ===") print("✓ Always monitor performance and memory usage") print("✓ Implement robust error handling with fallback methods") print("✓ Check matrix conditioning before solving linear systems") print("✓ Use batch processing for large-scale analyses") print("✓ Validate results and check for numerical issues") print("✓ Consider memory-efficient algorithms for large datasets")