import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime
import numpy as np

# Set style for better-looking plots
plt.style.use('default')
sns.set_palette("husl")

# Set matplotlib to use non-interactive backend to ensure plots are saved
plt.ioff()

def load_and_plot_data():
    """Load the RCarson CSV data and create comprehensive plots with standard deviation"""
    
    # Load the data
    try:
        df = pd.read_csv('data/RCarson_20250812_GFET001_A001001_noniso.csv')
        print(f"Data loaded successfully! Shape: {df.shape}")
        print(f"Columns: {list(df.columns)}")
        print(f"Data range: {df['Datetime'].min()} to {df['Datetime'].max()}")
        print(f"Number of measurements: {len(df)}")
    except FileNotFoundError:
        print("Error: Could not find the CSV file. Please check the path.")
        return
    except Exception as e:
        print(f"Error loading data: {e}")
        return
    
    # Convert Datetime to datetime object
    df['Datetime'] = pd.to_datetime(df['Datetime'])
    
    # Clean the data by removing rows where Vrse = -2.048 (default/invalid value)
    initial_count = len(df)
    df_cleaned = df[df['Vrse'] != -2.048].copy()
    removed_count = initial_count - len(df_cleaned)
    
    print(f"\nData Cleaning:")
    print(f"  Initial rows: {initial_count}")
    print(f"  Rows with Vrse = -2.048: {removed_count}")
    print(f"  Cleaned rows: {len(df_cleaned)}")
    print(f"  Data reduction: {removed_count/initial_count*100:.1f}%")
    
    # Trim the first 60 samples to remove initial sensor stabilization data
    if len(df_cleaned) > 60:
        df_trimmed = df_cleaned.iloc[60:].copy()
        trimmed_count = len(df_cleaned) - len(df_trimmed)
        print(f"\nData Trimming:")
        print(f"  Rows before trimming: {len(df_cleaned)}")
        print(f"  First 60 samples removed: {trimmed_count}")
        print(f"  Final trimmed rows: {len(df_trimmed)}")
        print(f"  Total data reduction: {(removed_count + trimmed_count)/initial_count*100:.1f}%")
    else:
        df_trimmed = df_cleaned.copy()
        trimmed_count = 0
        print(f"\nData Trimming:")
        print(f"  Warning: Only {len(df_cleaned)} rows available, cannot trim 60 samples")
        print(f"  No trimming applied")
    
    # Use trimmed data for plotting
    df = df_trimmed
    
    # Create a comprehensive figure with subplots
    fig, axes = plt.subplots(3, 2, figsize=(15, 12))
    fig.suptitle('RCarson GliderFET001 Non-ISO Sensor Data Analysis (Cleaned & Trimmed Data with Standard Deviation)', fontsize=16, fontweight='bold')
    
    # Plot 1: Vrse (Reference Electrode Voltage) over time with std
    axes[0, 0].plot(df['Datetime'], df['Vrse'], 'r-', linewidth=1, alpha=0.8, label='Vrse')
    if 'Vrse_std' in df.columns and df['Vrse_std'].sum() > 0:  # Only plot std if it has non-zero values
        axes[0, 0].fill_between(df['Datetime'], 
                                df['Vrse'] - df['Vrse_std'], 
                                df['Vrse'] + df['Vrse_std'], 
                                alpha=0.3, color='red', label='±1σ')
    axes[0, 0].set_title('Vrse (Reference Electrode Voltage) Over Time (Cleaned & Trimmed)', fontweight='bold')
    axes[0, 0].set_ylabel('Vrse (V)')
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].tick_params(axis='x', rotation=45)
    axes[0, 0].legend()
    
    # Plot 2: Vk (Working Electrode Voltage) over time with std
    axes[0, 1].plot(df['Datetime'], df['Vk'], 'm-', linewidth=1, alpha=0.8, label='Vk')
    if 'Vk_std' in df.columns and df['Vk_std'].sum() > 0:  # Only plot std if it has non-zero values
        axes[0, 1].fill_between(df['Datetime'], 
                                df['Vk'] - df['Vk_std'], 
                                df['Vk'] + df['Vk_std'], 
                                alpha=0.3, color='magenta', label='±1σ')
    axes[0, 1].set_title('Vk (Working Electrode Voltage) Over Time (Cleaned & Trimmed)', fontweight='bold')
    axes[0, 1].set_ylabel('Vk (V)')
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].tick_params(axis='x', rotation=45)
    axes[0, 1].legend()
    
    # Plot 3: VbiasPos (Bias Voltage) over time
    axes[1, 0].plot(df['Datetime'], df['VbiasPos'], 'purple', linewidth=1, alpha=0.8)
    axes[1, 0].set_title('VbiasPos (Bias Voltage) Over Time (Cleaned & Trimmed)', fontweight='bold')
    axes[1, 0].set_ylabel('VbiasPos (V)')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].tick_params(axis='x', rotation=45)
    
    # Plot 4: Current measurements (Ik and Ib) over time with dual y-axes
    ax4 = axes[1, 1]
    ax4_twin = ax4.twinx()
    
    # Plot Ik on left y-axis
    line1 = ax4.plot(df['Datetime'], df['Ik'], 'c-', linewidth=1, alpha=0.8, label='Ik (Working Current)')
    ax4.set_ylabel('Ik (Working Current) - nA', color='c')
    ax4.tick_params(axis='y', labelcolor='c')
    
    # Plot Ib on right y-axis
    line2 = ax4_twin.plot(df['Datetime'], df['Ib'], 'orange', linewidth=1, alpha=0.8, label='Ib (Bias Current)')
    ax4_twin.set_ylabel('Ib (Bias Current) - nA', color='orange')
    ax4_twin.tick_params(axis='y', labelcolor='orange')
    
    ax4.set_title('Current Measurements Over Time (Cleaned & Trimmed)', fontweight='bold')
    ax4.grid(True, alpha=0.3)
    ax4.tick_params(axis='x', rotation=45)
    
    # Create combined legend
    lines = line1 + line2
    labels = [l.get_label() for l in lines]
    ax4.legend(lines, labels, loc='upper left')
    
    # Plot 5: Standard deviation analysis - Enhanced visualization
    if 'Vrse_std' in df.columns and 'Vk_std' in df.columns:
        # Create a more prominent standard deviation plot
        axes[2, 0].plot(df['Datetime'], df['Vrse_std'], linewidth=2, alpha=0.8, label='Vrse std', color='red')
        axes[2, 0].plot(df['Datetime'], df['Vk_std'], linewidth=2, alpha=0.8, label='Vk std', color='magenta')
        
        # Add horizontal lines for mean values
        vrse_std_mean = df['Vrse_std'].mean()
        vk_std_mean = df['Vk_std'].mean()
        axes[2, 0].axhline(y=vrse_std_mean, color='red', linestyle='--', alpha=0.6, label=f'Vrse std mean: {vrse_std_mean:.6f}')
        axes[2, 0].axhline(y=vk_std_mean, color='magenta', linestyle='--', alpha=0.6, label=f'Vk std mean: {vk_std_mean:.6f}')
        
        # Add statistics text
        axes[2, 0].text(0.02, 0.98, f'Vrse std: μ={vrse_std_mean:.6f}, σ={df["Vrse_std"].std():.6f}', 
                        transform=axes[2, 0].transAxes, fontsize=9, verticalalignment='top',
                        bbox=dict(boxstyle='round', facecolor='red', alpha=0.3))
        axes[2, 0].text(0.02, 0.90, f'Vk std: μ={vk_std_mean:.6f}, σ={df["Vk_std"].std():.6f}', 
                        transform=axes[2, 0].transAxes, fontsize=9, verticalalignment='top',
                        bbox=dict(boxstyle='round', facecolor='magenta', alpha=0.3))
        
        axes[2, 0].set_title('Standard Deviation Over Time (Enhanced)', fontweight='bold', fontsize=12)
        axes[2, 0].set_ylabel('Standard Deviation (V)', fontsize=11)
        axes[2, 0].set_xlabel('Time', fontsize=11)
        axes[2, 0].grid(True, alpha=0.3)
        axes[2, 0].tick_params(axis='x', rotation=45)
        axes[2, 0].legend(fontsize=9, loc='upper right')
        
        # Set y-axis limits to focus on the std values
        y_max = max(df['Vrse_std'].max(), df['Vk_std'].max()) * 1.1
        axes[2, 0].set_ylim(0, y_max)
        
    else:
        # If no std data, show a different plot
        axes[2, 0].text(0.5, 0.5, 'No Standard Deviation Data Available', 
                        ha='center', va='center', transform=axes[2, 0].transAxes, fontsize=12)
        axes[2, 0].set_title('Standard Deviation Over Time', fontweight='bold')
    
    # Plot 6: Data quality and statistics
    axes[2, 1].text(0.1, 0.8, f'Total Measurements: {len(df)}', fontsize=10, transform=axes[2, 1].transAxes)
    axes[2, 1].text(0.1, 0.7, f'Time Range: {df["Datetime"].min().strftime("%Y-%m-%d %H:%M")} to {df["Datetime"].max().strftime("%Y-%m-%d %H:%M")}', 
                    fontsize=10, transform=axes[2, 1].transAxes)
    axes[2, 1].text(0.1, 0.6, f'Vrse Range: {df["Vrse"].min():.3f} to {df["Vrse"].max():.3f} V', 
                    fontsize=10, transform=axes[2, 1].transAxes)
    axes[2, 1].text(0.1, 0.5, f'Vk Range: {df["Vk"].min():.3f} to {df["Vk"].max():.3f} V', 
                    fontsize=10, transform=axes[2, 1].transAxes)
    axes[2, 1].text(0.1, 0.4, f'Ik Range: {df["Ik"].min():.1f} to {df["Ik"].max():.1f} nA', 
                    fontsize=10, transform=axes[2, 1].transAxes)
    axes[2, 1].text(0.1, 0.3, f'VbiasPos: {df["VbiasPos"].iloc[0]:.2f} V (constant)', 
                    fontsize=10, transform=axes[2, 1].transAxes)
    axes[2, 1].text(0.1, 0.2, f'Removed {removed_count} rows with Vrse=-2.048', 
                    fontsize=9, transform=axes[2, 1].transAxes, color='red')
    axes[2, 1].text(0.1, 0.1, f'Trimmed first {trimmed_count} samples', 
                    fontsize=9, transform=axes[2, 1].transAxes, color='blue')
    axes[2, 1].set_title('Data Summary (Cleaned & Trimmed)', fontweight='bold')
    axes[2, 1].axis('off')
    
    # Adjust layout
    plt.tight_layout()
    
    # Save the plot
    plt.savefig('rcarson_gliderfet_data_analysis_cleaned_trimmed.png', dpi=300, bbox_inches='tight')
    print("Main plot saved as 'rcarson_gliderfet_data_analysis_cleaned_trimmed.png'")
    
    # Close the figure to free memory
    plt.close(fig)
    
    # Create additional analysis plots
    create_additional_plots(df, removed_count, trimmed_count)

def create_additional_plots(df, removed_count, trimmed_count):
    """Create additional analysis plots for RCarson data"""
    
    try:
        # Create a new figure for correlation and distribution analysis
        fig2, axes2 = plt.subplots(2, 3, figsize=(18, 10))  # Changed to 2x3 layout
        fig2.suptitle('RCarson Data Additional Analysis (Cleaned & Trimmed Data)', fontsize=16, fontweight='bold')
        
        # Plot 1: Vrse vs Vk scatter plot
        axes2[0, 0].scatter(df['Vrse'], df['Vk'], alpha=0.6, s=20)
        axes2[0, 0].set_xlabel('Vrse (V)')
        axes2[0, 0].set_ylabel('Vk (V)')
        axes2[0, 0].set_title('Vrse vs Vk (Cleaned & Trimmed)', fontweight='bold')
        axes2[0, 0].grid(True, alpha=0.3)
        
        # Plot 2: Ik vs Ib scatter plot
        axes2[0, 1].scatter(df['Ik'], df['Ib'], alpha=0.6, s=20)
        axes2[0, 1].set_xlabel('Ik (Working Current) - nA')
        axes2[0, 1].set_ylabel('Ib (Bias Current) - nA')
        axes2[0, 1].set_title('Ik vs Ib (Cleaned & Trimmed)', fontweight='bold')
        axes2[0, 1].grid(True, alpha=0.3)
        
        # Plot 3: Vrse distribution histogram
        axes2[0, 2].hist(df['Vrse'], bins=30, alpha=0.7, color='red', edgecolor='black')
        axes2[0, 2].set_xlabel('Vrse (V)')
        axes2[0, 2].set_ylabel('Frequency')
        axes2[0, 2].set_title('Vrse Distribution (Cleaned & Trimmed)', fontweight='bold')
        axes2[0, 2].grid(True, alpha=0.3)
        
        # Plot 4: Ik distribution histogram
        axes2[1, 0].hist(df['Ik'], bins=50, alpha=0.7, color='cyan', edgecolor='black')
        axes2[1, 0].set_xlabel('Ik (Working Current) - nA')
        axes2[1, 0].set_ylabel('Frequency')
        axes2[1, 0].set_title('Ik Distribution (Cleaned & Trimmed)', fontweight='bold')
        axes2[1, 0].grid(True, alpha=0.3)
        
        # Plot 5: Vrse_std distribution histogram
        if 'Vrse_std' in df.columns:
            axes2[1, 1].hist(df['Vrse_std'], bins=30, alpha=0.7, color='darkred', edgecolor='black')
            axes2[1, 1].set_xlabel('Vrse Standard Deviation (V)')
            axes2[1, 1].set_ylabel('Frequency')
            axes2[1, 1].set_title('Vrse Standard Deviation Distribution', fontweight='bold')
            axes2[1, 1].grid(True, alpha=0.3)
            # Add statistics text
            vrse_std_mean = df['Vrse_std'].mean()
            vrse_std_std = df['Vrse_std'].std()
            axes2[1, 1].text(0.02, 0.98, f'μ={vrse_std_mean:.6f}\nσ={vrse_std_std:.6f}', 
                            transform=axes2[1, 1].transAxes, fontsize=10, verticalalignment='top',
                            bbox=dict(boxstyle='round', facecolor='darkred', alpha=0.3))
        
        # Plot 6: Vk_std distribution histogram
        if 'Vk_std' in df.columns:
            axes2[1, 2].hist(df['Vk_std'], bins=30, alpha=0.7, color='darkmagenta', edgecolor='black')
            axes2[1, 2].set_xlabel('Vk Standard Deviation (V)')
            axes2[1, 2].set_ylabel('Frequency')
            axes2[1, 2].set_title('Vk Standard Deviation Distribution', fontweight='bold')
            axes2[1, 2].grid(True, alpha=0.3)
            # Add statistics text
            vk_std_mean = df['Vk_std'].mean()
            vk_std_std = df['Vk_std'].std()
            axes2[1, 2].text(0.02, 0.98, f'μ={vk_std_mean:.6f}\nσ={vk_std_std:.6f}', 
                            transform=axes2[1, 2].transAxes, fontsize=10, verticalalignment='top',
                            bbox=dict(boxstyle='round', facecolor='darkmagenta', alpha=0.3))
        
        plt.tight_layout()
        
        # Save the additional plots
        plt.savefig('rcarson_additional_analysis_cleaned_trimmed.png', dpi=300, bbox_inches='tight')
        print("Additional plots saved as 'rcarson_additional_analysis_cleaned_trimmed.png'")
        
        # Close the figure to free memory
        plt.close(fig2)
        
    except Exception as e:
        print(f"Error creating additional plots: {e}")
    
    # Print summary statistics
    print_summary_statistics(df, removed_count, trimmed_count)

def print_summary_statistics(df, removed_count, trimmed_count):
    """Print summary statistics for the RCarson data"""
    
    print("\n" + "="*60)
    print("RCARSON DATA SUMMARY STATISTICS (CLEANED & TRIMMED)")
    print("="*60)
    
    # Data processing summary
    print(f"\nData Processing Summary:")
    print(f"  Rows removed (Vrse = -2.048): {removed_count}")
    print(f"  First 60 samples trimmed: {trimmed_count}")
    print(f"  Total data reduction: {(removed_count + trimmed_count)/(len(df) + removed_count + trimmed_count)*100:.1f}%")
    
    # Basic statistics for numerical columns
    numerical_cols = ['VbiasPos', 'Vrse', 'Vrse_std', 'Vk', 'Vk_std', 'Ik', 'Ib']
    
    for col in numerical_cols:
        if col in df.columns:
            print(f"\n{col}:")
            print(f"  Mean: {df[col].mean():.6f}")
            print(f"  Std:  {df[col].std():.6f}")
            print(f"  Min:  {df[col].min():.6f}")
            print(f"  Max:  {df[col].max():.6f}")
    
    # Time range information
    print(f"\nTime Range:")
    print(f"  Start: {df['Datetime'].min()}")
    print(f"  End:   {df['Datetime'].max()}")
    print(f"  Duration: {df['Datetime'].max() - df['Datetime'].min()}")
    
    # Data quality information
    print(f"\nData Quality:")
    print(f"  Total measurements (final): {len(df)}")
    print(f"  Missing values: {df.isnull().sum().sum()}")
    
    # Calculate measurement frequency
    time_diff = df['Datetime'].diff().dropna()
    avg_interval = time_diff.mean()
    print(f"  Average measurement interval: {avg_interval}")
    
    # Standard deviation analysis
    if 'Vrse_std' in df.columns and 'Vk_std' in df.columns:
        print(f"\nStandard Deviation Analysis:")
        print(f"  Vrse std - Mean: {df['Vrse_std'].mean():.6f}, Max: {df['Vrse_std'].max():.6f}")
        print(f"  Vk std - Mean: {df['Vk_std'].mean():.6f}, Max: {df['Vk_std'].max():.6f}")
        
        # Check if std values are meaningful
        if df['Vrse_std'].sum() == 0 and df['Vk_std'].sum() == 0:
            print(f"  Note: All standard deviation values are 0 - this may indicate single measurements")
        else:
            print(f"  Note: Standard deviation values vary - indicates multiple measurements per timepoint")
    
    # Vrse analysis after processing
    print(f"\nVrse Analysis (After Processing):")
    print(f"  Valid Vrse range: {df['Vrse'].min():.3f} to {df['Vrse'].max():.3f} V")
    print(f"  Vrse values above -2.0V: {(df['Vrse'] > -2.0).sum()}")
    print(f"  Vrse values below -2.0V: {(df['Vrse'] <= -2.0).sum()}")
    
    # Data processing benefits
    print(f"\nData Processing Benefits:")
    print(f"  Removed sensor initialization artifacts (Vrse = -2.048)")
    print(f"  Eliminated early sensor stabilization period (first 60 samples)")
    print(f"  Improved data quality for analysis")

if __name__ == "__main__":
    # Run the main analysis
    load_and_plot_data()