# -*- coding: utf-8 -*-
"""
Created on Thu Feb 19 15:29:05 2026

@author: chuffard
"""
from scipy import stats
import numpy as np
import seaborn as sns
import os
import glob
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from datetime import datetime
import sys
import os.path
import warnings
# Path to the folder
from pandas.api.types import is_datetime64_any_dtype as is_datetime
path = '//atlas/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2025 MARS deployment/2025July_Nov'

# Get all matching file paths as a list
atn_files = glob.glob(path+'/*sip*')

atn = pd.concat((pd.read_csv(f) for f in atn_files), ignore_index=True)
atn.columns
atn.dtypes
#make atn numerical
atn['atn'] = pd.to_numeric(atn['atn'], errors='coerce')

atn.to_csv(path+'/Station_M_fluoro.csv')


# Get all matching file paths as a list
fluoro_files = glob.glob(path+'/*fluorom*')

fluorescence = pd.concat((pd.read_csv(f) for f in fluoro_files), ignore_index=True)
fluorescence.columns
# Index(['SES Timestamp (ISO601)', 'Dark Ref Integration Time',
#        'Dark Ref 690 Counts', 'Dark Ref 590 Counts',
#        'Measure Ref Integration Time', 'Measure Ref 690 Counts',
#        'Measure Ref 590 Counts', 'Luminescense Integration Time',
#        'Luminescense 690 Counts', 'Luminescense 590 Counts',
#        'Fluorescense Integration Time', 'Fluorescense 690 Counts',
#        'Fluorescense 590 Counts', 'Phosphorescense Integration Time',
#        'Phosphorescense 690 Counts', 'Phosphorescense 590 Counts',
#        'Constant Target Fluorescense Integration Time',
#        'Constant Target Fluorescense 690 Counts',
#        'Constant Target Fluorescense 590 Counts', 'Code'],
#       dtype='object')

fluorescence['Phosphorescense_690'] = fluorescence['Phosphorescense 690 Counts']-fluorescence['Dark Ref 690 Counts']
fluorescence['Phosphorescense_590'] = fluorescence['Phosphorescense 590 Counts']-fluorescence['Dark Ref 590 Counts']

fluorescence['Fluorescense_690'] = fluorescence['Fluorescense 690 Counts']-fluorescence['Dark Ref 690 Counts']
fluorescence['Fluorescense_590'] = fluorescence['Fluorescense 590 Counts']-fluorescence['Dark Ref 590 Counts']

fluorescence['Luminescense_690'] = fluorescence['Luminescense 690 Counts']-fluorescence['Dark Ref 690 Counts']
fluorescence['Luminescense_590'] = fluorescence['Luminescense 590 Counts']-fluorescence['Dark Ref 590 Counts']

fluorescence['Constant_690'] = fluorescence['Constant Target Fluorescense 690 Counts']-fluorescence['Dark Ref 690 Counts']
fluorescence['Constant_590'] = fluorescence['Constant Target Fluorescense 590 Counts']-fluorescence['Dark Ref 590 Counts']

fluorescence= fluorescence[['SES Timestamp (ISO601)', 'Phosphorescense_690', 'Phosphorescense_590','Fluorescense_690','Fluorescense_590',
                           'Luminescense_690','Luminescense_590', 'Constant_590', 'Constant_690']]

# Get all matching file paths as a list
red_pH_files = glob.glob(path+'/*red tapeTRIM*')

red_pH = pd.concat((pd.read_csv(f, sep='\t',encoding='Windows-1252') for f in red_pH_files), ignore_index=True)
red_pH['height above seafloor (m)_red']= 1.3
red_pH.columns
# Index(['DateTime (YYYY-MM-DD hh:mm:ss)', 'Fraction of Second (ms)', 'Comment',
#        'status', 'dphi (0.001 °)', 'umolar (0.001 umol/L)',
#        'mbar (0.001 mbar)', 'airSat (0.001 %air sat)', 'tempSample (0.001 °C)',
#        'tempCase (0.001 °C)', 'signalIntensity (0.001 mV)',
#        'ambientLight (0.001 mV)', 'pressure (0.001 mbar)',
#        'humidity (0.001 %RH)', 'resistorTemp (0.001 Ohm or 0.001 mV)',
#        'percentO2 (0.001 %O2)', 'tempOptical (0.001 °C)', 'pH (0.001 pH)',
#        'R (10e-6)', '-', '-.1', 'height above seafloor (m)'],
#       dtype='object')

red_pH=red_pH[['DateTime (YYYY-MM-DD hh:mm:ss)','tempOptical (0.001 °C)','pressure (0.001 mbar)',
              'pH (0.001 pH)' ,'height above seafloor (m)_red']]

blue_pH_files = glob.glob(path+'/*blue tapeTRIM*')

blue_pH = pd.concat((pd.read_csv(f, sep='\t',encoding='Windows-1252') for f in blue_pH_files), ignore_index=True)
blue_pH['height above seafloor (m)_blue']= 0.3
blue_pH.columns

# Index(['DateTime (YYYY-MM-DD hh:mm:ss)', 'Fraction of Second (ms)', 'Comment',
#        'status', 'dphi (0.001 °)', 'umolar (0.001 umol/L)',
#        'mbar (0.001 mbar)', 'airSat (0.001 %air sat)', 'tempSample (0.001 °C)',
#        'tempCase (0.001 °C)', 'signalIntensity (0.001 mV)',
#        'ambientLight (0.001 mV)', 'pressure (0.001 mbar)',
#        'humidity (0.001 %RH)', 'resistorTemp (0.001 Ohm or 0.001 mV)',
#        'percentO2 (0.001 %O2)', 'tempOptical (0.001 °C)', 'pH (0.001 pH)',
#        'R (10e-6)', '-', '-.1', 'height above seafloor (m)'],
#       dtype='object')
blue_pH=blue_pH[['DateTime (YYYY-MM-DD hh:mm:ss)','tempOptical (0.001 °C)','pressure (0.001 mbar)',
              'pH (0.001 pH)' ,'height above seafloor (m)_blue']]
               
# df_minute = []
# def datetime_processing_avg(df, timecol):
#     df['timestamp_colNone'] = df[timecol].dt.tz_localize(None) #native
#     df['timestamp_colUTC'] = df['timestamp_colNone'].dt.tz_localize('UTC')#UTC
#     df['Local_time_LA'] = df['timestamp_colUTC'].dt.tz_convert('America/Los_Angeles')#America/Los_Angeles
#     df['RoundedMinute'] = df['Local_time_LA'].dt.floor('min')#UTC
#     df_minute = df.groupby(['RoundedMinute']).mean()
#     new_name = f"{df}_minute"
#     return df_minute, new_name

def datetime_processing_avg(df, time_col):
#def datetime_processing_avg(df, time_col, target_tz):   
    """
    1. Converts a time column to a new timezone.
    2. Rounds timestamps to the nearest hour.
    3. Calculates the mean of all other data columns, averaged by hour.
    """
    # Ensure column is datetime and handle timezone conversion
    # If naive, localize first; if aware, use tz_convert
    df["DateTime"] = pd.to_datetime(df[time_col])
    
    # if df[time_col].dt.tz is None:
    #     df[time_col] = df[time_col].dt.tz_localize('America/Los_Angeles') # Assume UTC if naive
    
    # Change timezone
 #   df[time_col] = df[time_col].dt.tz_convert(target_tz)
    
    # Round to the nearest hour
    df["DateTime_Hour"] = df["DateTime"].dt.floor('h', ambiguous='NaT')
    
    # Group by the rounded hour and average all numeric data
    # resample() is the preferred way for time-based groupby
    hourly_averaged_df = df.resample('h', on="DateTime_Hour").mean(numeric_only=True)
    
    return hourly_averaged_df

blue_pH_hour = datetime_processing_avg(blue_pH, 'DateTime (YYYY-MM-DD hh:mm:ss)')
blue_pH_hour.columns
# Index(['tempOptical (0.001 °C)', 'pressure (0.001 mbar)', 'pH (0.001 pH)',
#        'height above seafloor (m)_blue'],
#       dtype='object')
cols_to_suffix = ['tempOptical (0.001 °C)', 'pressure (0.001 mbar)', 'pH (0.001 pH)']
# Add '_new' suffix only to columns A and B
blue_pH_hour = blue_pH_hour.rename(columns={c: c + '_blue' for c in cols_to_suffix})
red_pH_hour = datetime_processing_avg(red_pH, 'DateTime (YYYY-MM-DD hh:mm:ss)')
red_pH_hour = red_pH_hour.rename(columns={c: c + '_red' for c in cols_to_suffix})
fluoro_hour = datetime_processing_avg(fluorescence, 'SES Timestamp (ISO601)')
atn_hour = datetime_processing_avg(atn, 'Datetime')
##give blue and red their own data columns

#atenuance is already once an hour
#dfs = [blue_pH_hour, red_pH_hour, fluoro_hour, atn_hour] 
merged_df = pd.concat([blue_pH_hour, red_pH_hour, fluoro_hour, atn_hour], axis=1)
merged_df= merged_df.reset_index()
#doesnt have atn
#####atn times don't match up. start here and78iu fix tomorrow
#Trim to columns within the deployment period only
deploy_range = pd.date_range(start='2025-07-17', end='2025-11-09')

merged_df = merged_df[merged_df['DateTime_Hour'].between('2025-07-17', '2025-11-09')]

numeric_cols = merged_df.select_dtypes(include=['number']).columns
datetime_column_name = merged_df['DateTime_Hour']
is_datetime(merged_df['DateTime_Hour'])

# 3. Save plots to a multi-page PDF
with PdfPages('all_sensor_plots.pdf') as pdf:
    for col in numeric_cols:
        plt.figure(figsize=(10, 6)) # Create a new figure
        plt.plot(datetime_column_name, merged_df[col], label=col)
        plt.title(f'Line Plot of {col}')
   #     plt.xlabel('DateTime_Hour')
        plt.ylabel('Value')
        plt.legend()
        plt.xticks(rotation=45) 
        plt.grid(True)
        # Save the current figure to the PDF
        pdf.savefig()
        plt.close() # Close to free up memory

print("PDF generated successfully: all_sensor_plots.pdf")

print(os.getcwd())

merged_df.columns
# Index(['DateTime_Hour', 'tempOptical (0.001 °C)_blue',
#        'pressure (0.001 mbar)_blue', 'pH (0.001 pH)_blue',
#        'height above seafloor (m)_blue', 'tempOptical (0.001 °C)_red',
#        'pressure (0.001 mbar)_red', 'pH (0.001 pH)_red',
#        'height above seafloor (m)_red', 'Phosphorescense_690',
#        'Phosphorescense_590', 'Fluorescense_690', 'Fluorescense_590',
#        'Luminescense_690', 'Luminescense_590', 'Constant_590', 'Constant_690',
#        'atn', 'Number in background'],
#       dtype='str')


#NEED TO FILL IN ATN GAP
##########################Basic plotting
fig, ax = plt.subplots()

# Plot each series with distinct properties and labels
ax.scatter(merged_df['Fluorescense_690'], merged_df['pH (0.001 pH)_red'], label='pH red', color='red', marker='o', alpha = 0.3) # 'o' for circle marker
ax.scatter(merged_df['Fluorescense_690'], merged_df['pH (0.001 pH)_blue'], label='pH blue', color='blue', marker='s', alpha = 0.3)  # 's' for square marker

# Add labels, title, and legend
ax.set_xlabel('fluoro chl')
ax.set_ylabel('pH')
ax.legend() # Displays the labels for each series

# Display the plot
plt.show()
################
fig, ax = plt.subplots()

# Plot each series with distinct properties and labels
ax.plot(merged_df['date_time_format'], merged_df['pH (0.001 pH)_red'], label='pH red', color='red', marker='o', alpha = 0.3) # 'o' for circle marker
ax.plot(merged_df['date_time_format'], merged_df['pH (0.001 pH)_blue'], label='pH blue', color='blue', marker='s', alpha = 0.3)  # 's' for square marker

# Add labels, title, and legend
ax.set_xlabel('date')
ax.set_ylabel('pH')
ax.legend() # Displays the labels for each series

# Display the plot
plt.show()

################

ax = sns.scatterplot(data=merged_df, x='pH (0.001 pH)_red', y='pH (0.001 pH)_blue', hue='Fluorescense_690', palette='viridis', alpha = 0.3, legend=True)

# Create a mappable to sync the colorbar with the plot
norm = plt.Normalize(merged_df['Fluorescense_690'].min(), merged_df['Fluorescense_690'].max())
sm = plt.cm.ScalarMappable(cmap="viridis", norm=norm)
sm.set_array([])

output_path = "//atlas/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2025 MARS deployment/2025July_Nov/pHscatter_plot.pdf"
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.show()



################

ax = sns.scatterplot(data=merged_df, x='DateTime_Hour', y='pH (0.001 pH)_red', hue='Fluorescense_690', palette='viridis',alpha = 0.3, legend=True)

# Create a mappable to sync the colorbar with the plot
norm = plt.Normalize(merged_df['Fluorescense_690'].min(), merged_df['Fluorescense_690'].max())
sm = plt.cm.ScalarMappable(cmap="viridis", norm=norm)
sm.set_array([])

output_path = "//atlas/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2025 MARS deployment/2025July_Nov/pH_red.pdf"
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.show()



################

ax = sns.scatterplot(data=merged_df, x='date_time_format', y='pH (0.001 pH)_blue', hue='Fluorescense_690', palette='viridis',alpha = 0.3,legend=True)

# Create a mappable to sync the colorbar with the plot
norm = plt.Normalize(merged_df['Fluorescense_690'].min(), merged_df['Fluorescense_690'].max())
sm = plt.cm.ScalarMappable(cmap="viridis", norm=norm)
plt.xticks(rotation=45) 
sm.set_array([])

output_path = "//atlas/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2025 MARS deployment/2025July_Nov/pH_blue.pdf"
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.xticks(rotation=45) 
plt.show()

#%%
#### 
use_min_period = 24
atn = atn.sort_values(by='Datetime')

atn['rm150atn'] = atn['atn'].rolling(150, min_periods=23).mean()
atn['Slope_rm_150atn'] = atn['rm150atn'].diff()

##what percentil are these values
atn['Percentile_Rank'] = atn['Slope_rm_150atn'].rank(pct=True) * 100

lowtrigger = 0.007 #97.05 th percentile (1 away from 97.1th)
hightrigger = 0.011# = 99.1th percentile
atn["triggerNrecslope_atn"] = np.where(atn['Slope_rm_150atn'] >= lowtrigger, 'LowTriggeratn', 'WaitSlopeatn')

nsamples_6weeks = 32
if (len(atn) > nsamples_6weeks):
    #find the rolling rate low triggers in past ~1.5 months
    atn['low_trigger_rate'] = (atn['triggerNrecslope_atn'] == 'LowTriggeratn').rolling(nsamples_6weeks, min_periods=2).sum()

    #Line to help keep tabs on whether the switch has been flipped and we're now using the high trigger

    atn['max_rate'] = atn['low_trigger_rate'].cummax()

    try:
        # Ensure 'max_rate' column is not empty and the values are valid
        if 'max_rate' in atn.columns and not atn['max_rate'].isnull().all():
            # Apply condition element-wise to the 'max_rate' column
            atn["Fire_slope_trigger"] = np.where(
                atn['max_rate'].abs() > 3,
                np.where(atn['Slope_rm_150atn'] >= hightrigger, 'Fire', 'Wait'),
                np.where(atn['Slope_rm_150atn'] >= lowtrigger, 'Fire', 'Wait')
            )
        else:
            print("The 'max_rate' column is missing or contains all NaN values.")
    except Exception as e:
        print(f"An error occurred: {e}")

#if doesn't fire in ~21 days, fire anyway
matchlist = ['Wait'] * 30
matchlist.insert(0,'Fire')

def check_sequence(df, col_name, seq):
    """
    Checks if a sequence of values in a column matches a given list.

    Args:
        df: The DataFrame.
        col_name: The name of the column to check.
        seq: The sequence of strings to match.

    Returns:
        A new DataFrame with a 'match' column indicating if the sequence was found.
    """

    df['match'] = False
    for i in range(len(df) - len(seq) + 1):
        if df[col_name].iloc[i:i + len(seq)].tolist() == seq:
            df.loc[i + len(seq) - 1, 'match'] = True
    return df


# Example usage:
df = atn
seq_to_match = matchlist

df = check_sequence(df, 'Fire_slope_trigger', seq_to_match)
df['match'] = df['match'].astype(str)
df["Fire_anyway"] = np.where(df['match'].str.contains('T'), 'Fire', 'Wait')

if (len(df) > nsamples_6weeks):
    df['SendTrigger'] = np.where((df['Fire_anyway'] == 'Fire') |
                        (df['Fire_slope_trigger'] =='Fire'), 'SendTrigger', 'Wait')
    result = np.where((df['Fire_anyway'] == 'Fire') |
             (df['Fire_slope_trigger'] =='Fire'), 'Event detected', 'Wait')
else:
    df['SendTrigger'] = np.where((df['Fire_anyway'] == 'Fire') |
                        False, 'SendTrigger', 'Wait')
    result = np.where((df['Fire_anyway'] == 'Fire') |
             False, 'Event detected', 'Wait')

# Emit result
if (len(result) > 0):
    print(result[len(result)-1])

df.to_csv('//atlas/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2025 MARS deployment/2025July_Nov/SES_triggered.csv', index=False)

df['Slope_rm_150atn'].hist(bins=100) 
df['z_score690slopeatn'] = stats.zscore(df['Slope_rm_150atn'], nan_policy='omit')

numeric_colsfluoro = df.select_dtypes(include=['number']).columns
datetime_column_name = df['Datetime']
is_datetime(df['Datetime'])

# 3. Save plots to a multi-page PDF
with PdfPages('all_sensor_plotsatn.pdf') as pdf:
    for col in numeric_colsfluoro:
        plt.figure(figsize=(10, 6)) # Create a new figure
        plt.plot(datetime_column_name, df[col], label=col)
        plt.title(f'Line Plot of {col}')
   #     plt.xlabel('DateTime_Hour')
        plt.ylabel('Value')
        plt.legend()
        plt.xticks(rotation=45) 
        plt.grid(True)
        # Save the current figure to the PDF
        pdf.savefig()
        plt.close() # Close to free up memory

print("PDF generated successfully: all_sensor_plots.pdf")

#%%
#! /usr/bin/env python

use_min_period = 24
fluoro = fluorescence.sort_values(by='DateTime')
fluoro.columns
fluoro['rm150fluoro'] = fluoro['Fluorescense_690'].rolling(150, min_periods=23).mean()
fluoro['Slope_rm_150fluoro'] = fluoro['rm150fluoro'].diff()
fluoro['z_score690slopefluoro'] = stats.zscore(df['Slope_rm_150fluoro'], nan_policy='omit')


threshold5 = fluoro['Slope_rm_150fluoro'].quantile(0.95)
threshold3 = fluoro['Slope_rm_150fluoro'].quantile(0.97)

lowtrigger =threshold5
hightrigger = threshold3
fluoro["triggerNrecslope_fluoro"] = np.where(fluoro['Slope_rm_150fluoro'] >= lowtrigger, 'LowTriggerfluoro', 'WaitSlopefluoro')

nsamples_6weeks = 32
if (len(fluoro) > nsamples_6weeks):
    #find the rolling rate low triggers in past ~1.5 months
    fluoro['low_trigger_rate'] = (fluoro['triggerNrecslope_fluoro'] == 'LowTriggerfluoro').rolling(nsamples_6weeks, min_periods=2).sum()

    #Line to help keep tabs on whether the switch has been flipped and we're now using the high trigger

    fluoro['max_rate'] = fluoro['low_trigger_rate'].cummax()

    try:
        # Ensure 'max_rate' column is not empty and the values are valid
        if 'max_rate' in fluoro.columns and not fluoro['max_rate'].isnull().all():
            # Apply condition element-wise to the 'max_rate' column
            fluoro["Fire_slope_trigger"] = np.where(
                fluoro['max_rate'].abs() > 3,
                np.where(fluoro['Slope_rm_150fluoro'] >= hightrigger, 'Fire', 'Wait'),
                np.where(fluoro['Slope_rm_150fluoro'] >= lowtrigger, 'Fire', 'Wait')
            )
        else:
            print("The 'max_rate' column is missing or contains all NaN values.")
    except Exception as e:
        print(f"An error occurred: {e}")

#if doesn't fire in ~21 days, fire anyway
matchlist = ['Wait'] * 30
matchlist.insert(0,'Fire')

def check_sequence(df, col_name, seq):
    """
    Checks if a sequence of values in a column matches a given list.

    Args:
        df: The DataFrame.
        col_name: The name of the column to check.
        seq: The sequence of strings to match.

    Returns:
        A new DataFrame with a 'match' column indicating if the sequence was found.
    """

    df['match'] = False
    for i in range(len(df) - len(seq) + 1):
        if df[col_name].iloc[i:i + len(seq)].tolist() == seq:
            df.loc[i + len(seq) - 1, 'match'] = True
    return df


# Example usage:
df = fluoro
seq_to_match = matchlist

df = check_sequence(df, 'Fire_slope_trigger', seq_to_match)
df['match'] = df['match'].astype(str)
df["Fire_anyway"] = np.where(df['match'].str.contains('T'), 'Fire', 'Wait')

if (len(df) > nsamples_6weeks):
    df['SendTrigger'] = np.where((df['Fire_anyway'] == 'Fire') |
                        (df['Fire_slope_trigger'] =='Fire'), 'SendTrigger', 'Wait')
    result = np.where((df['Fire_anyway'] == 'Fire') |
             (df['Fire_slope_trigger'] =='Fire'), 'Event detected', 'Wait')
else:
    df['SendTrigger'] = np.where((df['Fire_anyway'] == 'Fire') |
                        False, 'SendTrigger', 'Wait')
    result = np.where((df['Fire_anyway'] == 'Fire') |
             False, 'Event detected', 'Wait')

# Emit result
if (len(result) > 0):
    print(result[len(result)-1])

df.to_csv('//atlas/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2025 MARS deployment/2025July_Nov/SES_triggeredFluoro.csv', index=False)


numeric_colsfluoro = df.select_dtypes(include=['number']).columns
datetime_column_name = df['DateTime']
is_datetime(df['DateTime'])

# 3. Save plots to a multi-page PDF
with PdfPages('all_sensor_plotsfluoro.pdf') as pdf:
    for col in numeric_colsfluoro:
        plt.figure(figsize=(10, 6)) # Create a new figure
        plt.plot(datetime_column_name, df[col], label=col)
        plt.title(f'Line Plot of {col}')
   #     plt.xlabel('DateTime_Hour')
        plt.ylabel('Value')
        plt.legend()
        plt.xticks(rotation=45) 
        plt.grid(True)
        # Save the current figure to the PDF
        pdf.savefig()
        plt.close() # Close to free up memory

print("PDF generated successfully: all_sensor_plots.pdf")
