"""
Complete SciPy Scientific Computing Example
This example demonstrates a comprehensive workflow combining multiple SciPy modules
Author: RSK World
Website: https://rskworld.in
Email: help@rskworld.in
Phone: +91 93305 39277
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy import integrate
from scipy.interpolate import interp1d
from scipy import stats
from scipy import signal
from scipy.fft import fft, fftfreq

def main():
    """
    Complete workflow example combining optimization, integration, 
    interpolation, statistics, and signal processing
    """
    print("=" * 70)
    print("Complete SciPy Scientific Computing Example")
    print("Author: RSK World - https://rskworld.in")
    print("=" * 70 + "\n")
    
    # ===================================================================
    # 1. OPTIMIZATION: Find optimal parameters for a model
    # ===================================================================
    print("1. OPTIMIZATION: Finding optimal model parameters")
    print("-" * 70)
    
    # Generate noisy experimental data
    np.random.seed(42)
    x_data = np.linspace(0, 10, 20)
    y_true = 3 * np.exp(-0.5 * x_data) * np.sin(2 * x_data)
    y_data = y_true + 0.2 * np.random.randn(len(x_data))
    
    # Model function: y = a * exp(-b*x) * sin(c*x)
    def model(x, params):
        a, b, c = params
        return a * np.exp(-b * x) * np.sin(c * x)
    
    # Objective: minimize sum of squared errors
    def objective(params):
        y_pred = model(x_data, params)
        return np.sum((y_data - y_pred)**2)
    
    # Initial guess
    params0 = [2, 0.3, 1.5]
    
    # Optimize
    result = minimize(objective, params0, method='BFGS')
    optimal_params = result.x
    
    print(f"Optimal parameters: a={optimal_params[0]:.4f}, "
          f"b={optimal_params[1]:.4f}, c={optimal_params[2]:.4f}")
    print(f"Minimum error: {result.fun:.6f}\n")
    
    # ===================================================================
    # 2. INTEGRATION: Calculate area under the fitted curve
    # ===================================================================
    print("2. INTEGRATION: Calculating area under the fitted curve")
    print("-" * 70)
    
    def fitted_function(x):
        return model(x, optimal_params)
    
    # Integrate from 0 to 10
    area, error = integrate.quad(fitted_function, 0, 10)
    print(f"Area under curve from 0 to 10: {area:.6f}")
    print(f"Integration error estimate: {error:.2e}\n")
    
    # ===================================================================
    # 3. INTERPOLATION: Create smooth curve from fitted model
    # ===================================================================
    print("3. INTERPOLATION: Creating smooth interpolated curve")
    print("-" * 70)
    
    # Generate fine grid
    x_fine = np.linspace(0, 10, 200)
    y_fine = model(x_fine, optimal_params)
    
    # Create interpolation function
    interp_func = interp1d(x_fine, y_fine, kind='cubic')
    
    # Evaluate at new points
    x_new = np.linspace(0, 10, 50)
    y_interp = interp_func(x_new)
    
    print(f"Interpolated {len(x_new)} points from {len(x_fine)} original points\n")
    
    # ===================================================================
    # 4. STATISTICS: Analyze residuals
    # ===================================================================
    print("4. STATISTICS: Analyzing model residuals")
    print("-" * 70)
    
    # Calculate residuals
    y_pred = model(x_data, optimal_params)
    residuals = y_data - y_pred
    
    # Statistical analysis
    mean_residual = np.mean(residuals)
    std_residual = np.std(residuals)
    
    # Test if residuals are normally distributed
    _, p_value = stats.normaltest(residuals)
    
    print(f"Mean residual: {mean_residual:.6f}")
    print(f"Std residual: {std_residual:.6f}")
    print(f"Normality test p-value: {p_value:.6f}")
    print(f"Residuals are {'normally distributed' if p_value > 0.05 else 'not normally distributed'}\n")
    
    # ===================================================================
    # 5. SIGNAL PROCESSING: Filter the fitted signal
    # ===================================================================
    print("5. SIGNAL PROCESSING: Filtering the fitted signal")
    print("-" * 70)
    
    # Add noise to fitted signal
    y_noisy = y_fine + 0.1 * np.random.randn(len(y_fine))
    
    # Apply low-pass filter
    b, a = signal.butter(4, 0.1, 'low')
    y_filtered = signal.filtfilt(b, a, y_noisy)
    
    # Frequency analysis
    fft_vals = fft(y_filtered)
    fft_freq = fftfreq(len(y_filtered), x_fine[1] - x_fine[0])
    psd = np.abs(fft_vals)**2
    
    # Find dominant frequency
    positive_freq_idx = fft_freq > 0
    dominant_freq_idx = np.argmax(psd[positive_freq_idx])
    dominant_freq = fft_freq[positive_freq_idx][dominant_freq_idx]
    
    print(f"Dominant frequency: {dominant_freq:.4f} Hz")
    print(f"Expected frequency: {optimal_params[2]/(2*np.pi):.4f} Hz\n")
    
    # ===================================================================
    # 6. VISUALIZATION: Create comprehensive plot
    # ===================================================================
    print("6. VISUALIZATION: Creating comprehensive visualization")
    print("-" * 70)
    
    fig = plt.figure(figsize=(16, 12))
    
    # Plot 1: Data and fitted model
    ax1 = plt.subplot(3, 2, 1)
    ax1.scatter(x_data, y_data, color='red', s=50, alpha=0.7, label='Experimental data')
    ax1.plot(x_data, y_true, 'k--', linewidth=2, alpha=0.5, label='True function')
    ax1.plot(x_fine, y_fine, 'b-', linewidth=2, label='Fitted model')
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('y', fontsize=12)
    ax1.set_title('1. Optimization: Model Fitting', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Plot 2: Integration area
    ax2 = plt.subplot(3, 2, 2)
    ax2.plot(x_fine, y_fine, 'b-', linewidth=2, label='Fitted function')
    ax2.fill_between(x_fine, 0, y_fine, alpha=0.3, color='blue', 
                     label=f'Area = {area:.4f}')
    ax2.set_xlabel('x', fontsize=12)
    ax2.set_ylabel('y', fontsize=12)
    ax2.set_title('2. Integration: Area Calculation', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Plot 3: Interpolation
    ax3 = plt.subplot(3, 2, 3)
    ax3.plot(x_fine, y_fine, 'b-', linewidth=1, alpha=0.5, label='Original')
    ax3.scatter(x_new, y_interp, color='red', s=30, label='Interpolated points')
    ax3.set_xlabel('x', fontsize=12)
    ax3.set_ylabel('y', fontsize=12)
    ax3.set_title('3. Interpolation: Smooth Curve', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Plot 4: Residuals
    ax4 = plt.subplot(3, 2, 4)
    ax4.scatter(x_data, residuals, color='green', s=50, alpha=0.7)
    ax4.axhline(0, color='black', linestyle='--', linewidth=1)
    ax4.axhline(mean_residual, color='red', linestyle='--', linewidth=1, 
                label=f'Mean = {mean_residual:.4f}')
    ax4.fill_between(x_data, mean_residual - std_residual, mean_residual + std_residual,
                     alpha=0.2, color='red', label=f'±1 std = {std_residual:.4f}')
    ax4.set_xlabel('x', fontsize=12)
    ax4.set_ylabel('Residual', fontsize=12)
    ax4.set_title('4. Statistics: Residual Analysis', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # Plot 5: Signal filtering
    ax5 = plt.subplot(3, 2, 5)
    ax5.plot(x_fine, y_noisy, 'r-', linewidth=1, alpha=0.5, label='Noisy signal')
    ax5.plot(x_fine, y_filtered, 'b-', linewidth=2, label='Filtered signal')
    ax5.set_xlabel('x', fontsize=12)
    ax5.set_ylabel('y', fontsize=12)
    ax5.set_title('5. Signal Processing: Filtering', fontsize=12, fontweight='bold')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # Plot 6: Frequency spectrum
    ax6 = plt.subplot(3, 2, 6)
    positive_freq = fft_freq[positive_freq_idx]
    positive_psd = psd[positive_freq_idx]
    ax6.plot(positive_freq, positive_psd, 'b-', linewidth=2)
    ax6.axvline(dominant_freq, color='red', linestyle='--', linewidth=2,
                label=f'Dominant: {dominant_freq:.4f} Hz')
    ax6.set_xlabel('Frequency (Hz)', fontsize=12)
    ax6.set_ylabel('Power Spectral Density', fontsize=12)
    ax6.set_title('6. Signal Processing: Frequency Analysis', fontsize=12, fontweight='bold')
    ax6.set_xlim(0, 2)
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('complete_example.png', dpi=300, bbox_inches='tight')
    print("Complete visualization saved as 'complete_example.png'\n")
    plt.close()
    
    print("=" * 70)
    print("Complete example finished successfully!")
    print("=" * 70)


if __name__ == "__main__":
    main()