# -*- coding: utf-8 -*-
"""
Created on Tue Apr 23 09:15:57 2024

@author: chuffard
"""
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 16 15:47:46 2024

@author: chuffard
"""

# -*- coding: utf-8 -*-
"""
Created on Wed Dec 13 09:34:42 2023

@author: chuffard
"""

import numpy as np

import pandas as pd

#import spectrum #doesn't want to install...tried everything I can find

import matplotlib.pyplot as plt

#from spectrum import Periodogram, data_cosine

from datetime import datetime

import scipy as SP

import scipy.signal as signal

import seaborn as sns
import statistics
#import endaq as endaq

import statsmodels.tsa.stattools as sts

from statsmodels.tsa.stattools import acf

import pyflakes as pyf

#import lzip as lz

from pandas import read_csv

from statsmodels.graphics.tsaplots import plot_acf

from statsmodels.graphics.tsaplots import plot_pacf


import ephem

import pygam

import re

import numpy as np
from scipy import stats

from matplotlib.backends.backend_pdf import PdfPages

# read in data; use first line as header
#Main_Hourly_SES_MARS = pd.read_csv('G:/Analyses/MainSES_fluoro_atn_moon_sun_current.csv', parse_dates=True, index_col="datetime")

# read in data; use first line as header
SESdata = pd.read_csv('G:/Analyses/SES_out_atn_Mars2023Deployment.csv', parse_dates=True, index_col="Date_time")

#Now trim it to once an hour. Make three separate datasets 

SESonceanhour_2 = SESdata.iloc[2::3, :]


#remove infinite values
SESonceanhour_2= SESonceanhour_2[SESonceanhour_2['atn']!= np.inf]
SESonceanhour_2.dtypes
SEStrimSESonceanhour_2 = SESonceanhour_2.drop(['Date', 'Time'], axis=1)

SEStrimSESonceanhour_2["atn"].plot()
plt.title('Hour_2')

#fix weird time zone stuff
# for col in Main_Hourly_SES_MARS.select_dtypes(include=['datetime64[ns, UTC]']).columns:
#     Main_Hourly_SES_MARS[col] = Main_Hourly_SES_MARS[col].apply(lambda x: x.tz_localize(None))

#run a loop that uses running means from 24 to 336 hours, and see how many triggers it comes up with

count=0
R_atncount=[]
R_atnmeancount=[]
n_count=[]
dfcount=[]
SErcount=[]
CIhighr=[]
trigger=[]
trigger_count = []
wait_count = []
triggerlist=[]

#4 days to 18 days
for x in range(23, 336):
        data =SEStrimSESonceanhour_2
        variablename = str(x)
        data[variablename] = data['atn'].rolling(x, min_periods=23).mean()
        R_atncount = data['atn'].rolling(x, min_periods=23).mean()
        R_atnmeancount = R_atncount.mean()
        n_count = len(R_atncount)
        dfcount = n_count-1
        SErcount = stats.sem(R_atncount.dropna())
        CIlevel = 0.9999
        CI99atnrcount = stats.t.interval(CIlevel, dfcount, R_atnmeancount, scale = SErcount)
        CIhighr = float(CI99atnrcount[1])
        triggercol = str(x)
        data[triggercol]=np.where(data[variablename] >= CIhighr, 'Trigger', 'Wait')
        triggercount = data[triggercol].value_counts()['Trigger']
        waitcount = data[triggercol].value_counts()['Wait']
        totalcount = triggercount+waitcount
        triggerlist.append(triggercount)
    
triggerlistdf = pd.DataFrame(triggerlist)

triggerlistdf= triggerlistdf.rename(columns={0: "number_of_triggers"})
triggerlistdf = triggerlistdf.reset_index()
triggerlistdf['rollingmean_count'] = triggerlistdf['index']+23
#triggerlistdf= triggerlistdf.rename(columns={"index": "rollingmean_count"})
triggerlistdf.to_csv("triggerlistdf_2.csv")


sns.lineplot(data=triggerlistdf, x='rollingmean_count', y='number_of_triggers')
plt.title('Hour_2')

#find the rolling mean with the fewest triggers
#There's a local min at around 130 hours. 
#But first the chances of trigger are very low, so need to exclude this
triggerlistdf0_100 = triggerlistdf.iloc[100:]
rowwithlottrigger0_100 = triggerlistdf0_100[triggerlistdf0_100.number_of_triggers==triggerlistdf0_100.number_of_triggers.min()]
rollingmeanUSE0_100 = rowwithlottrigger0_100.rollingmean_count.min()
print(rollingmeanUSE0_100)
#132
#make a new time series of this rolling mean rolling mean 132 hours
SEStrimSESonceanhour_2['rolling_atn_132'] = SEStrimSESonceanhour_2['atn'].rolling(rollingmeanUSE0_100, min_periods=23).mean()

sns.lineplot(data=SEStrimSESonceanhour_2, x='Date_time', y='atn')
sns.lineplot(data=SEStrimSESonceanhour_2, x='Date_time', y='rolling_atn_132')
plt.xticks(rotation=45);
plt.title('Hour_2 rolling 132')


#Now check the slope of this time series
#subset columns
SESMARS132atn = SEStrimSESonceanhour_2[["atn", "rolling_atn_132"]]

#find proxy for slope (percent change from one sample to next)
SESMARS132atn['diff'] = SESMARS132atn['rolling_atn_132'].diff()


#subset to remove first 132 rows (to eliminate artifacts of using rolling mean)
SESslope132TRIM = SESMARS132atn.iloc[125:, ]

#find max slope after this
trigger132slope = max(SESslope132TRIM['diff'])
print(trigger132slope)
#0.013459882401515133

#make that new max the trigger
SESslope132TRIM["trigger132slope"] = np.where(SESslope132TRIM['diff'] >= trigger132slope, 'TriggerSlope', 'WaitSlope')

#this is a little high and picks up the peak. maybe we can get a little before the peak
#trigger 80% of max slope
trigger80pc132slope = 0.8 * trigger132slope
print(trigger80pc132slope)
#0.010767905921212108

#list which samples are trigger vs wait
SESslope132TRIM["trigger80pc132slope"] = np.where(SESslope132TRIM['diff'] >= trigger80pc132slope, 'TriggerSlope', 'WaitSlope')

#plot to verify.
sns.scatterplot(data=SESslope132TRIM, x="Date_time", y="atn", hue="trigger80pc132slope")
plt.xticks(rotation=45);
plt.show()

sns.scatterplot(data=SESslope132TRIM, x="Date_time", y="atn", hue="trigger132slope")
plt.xticks(rotation=45);
plt.show()

#find daily average of SES time series
daily_SES132atn = SESslope132TRIM['atn'].resample('1D').mean()

#turn into data frame
daily_SES132atnatndf = pd.DataFrame([daily_SES132atn]).T


testtrim1H_2 = SESslope132TRIM.loc['2023-05-08':'2023-08-25']  
testtrim1D132 = daily_SES132atnatndf.loc['2023-05-08':'2023-08-25']  
triggerH132 = testtrim1H_2[testtrim1H_2["trigger80pc132slope"].str.contains("Trigger")]
triggerHfull132 = testtrim1H_2[testtrim1H_2["trigger132slope"].str.contains("Trigger")]

fig, ax = plt.subplots()
sns.lineplot(data=SESdatatrim, x="Date_time", y="atn", ax=ax, color = "gray")
sns.lineplot(data=testtrim1H_2, x="Date_time", y="atn", ax=ax, color = "blue")
sns.lineplot(data=testtrim1D132, x="Date_time", y="atn", ax=ax, color = "orange")
sns.lineplot(data=daily_SESatntrim, x="Date_time", y="atn", ax=ax, color = "black")
sns.scatterplot(data=triggerH132, x="Date_time", y="atn", hue="trigger80pc132slope")
sns.scatterplot(data=triggerHfull132, x="Date_time", y="atn", color = "orange")
plt.title('Hour_2 rm132')
plt.xticks(rotation=45);