# -*- 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 trim140 it to once an hour. Make three separate datasets 

SESonceanhour_1 = SESdata.iloc[1::3, :]


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

SEStrim140SESonceanhour_1["atn"].plot()
plt.title('Hour_1')

#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 =SEStrim140SESonceanhour_1
        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_1.csv")


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

#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)
#140
#make a new time series of this rolling mean rolling mean 140 hours
SEStrim140SESonceanhour_1['rolling_atn_140'] = SEStrim140SESonceanhour_1['atn'].rolling(rollingmeanUSE0_100, min_periods=23).mean()

sns.lineplot(data=SEStrim140SESonceanhour_1, x='Date_time', y='atn')
sns.lineplot(data=SEStrim140SESonceanhour_1, x='Date_time', y='rolling_atn_140')
plt.xticks(rotation=45);
plt.title('Hour_1')


#Now check the slope of this time series
#subset columns
SESMARS140atn = SEStrim140SESonceanhour_1[["atn", "rolling_atn_140"]]

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


#subset to remove first 140 rows (to eliminate artifacts of using rolling mean)
SESslope140trim140 = SESMARS140atn.iloc[125:, ]

#find max slope after this
trigger140slope = max(SESslope140trim140['diff'])
print(trigger140slope)
#0.01438800002857138

#make that new max the trigger
SESslope140trim140["trigger140slope"] = np.where(SESslope140trim140['diff'] >= trigger140slope, '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
trigger80pc140slope = 0.8 * trigger140slope
print(trigger80pc140slope)
#0.011510400022857104

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

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

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

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

#turn into data frame
daily_SESatnatndf = pd.DataFrame([daily_SESatn]).T


testtrim1401H_1 = SESslope140trim140.loc['2023-05-08':'2023-08-25']  
testtrim1401D = daily_SESatnatndf.loc['2023-05-08':'2023-08-25']  
trigger140H = testtrim1401H_1[testtrim1401H_1["trigger80pc140slope"].str.contains("Trigger")]
trigger140Hfull = testtrim1401H_1[testtrim1401H_1["trigger140slope"].str.contains("Trigger")]

fig, ax = plt.subplots()
sns.lineplot(data=SESdatatrim140, x="Date_time", y="atn", ax=ax, color = "gray")
sns.lineplot(data=testtrim1401H_1, x="Date_time", y="atn", ax=ax, color = "blue")
sns.lineplot(data=testtrim1401D, x="Date_time", y="atn", ax=ax, color = "orange")
sns.lineplot(data=daily_SESatntrim140, x="Date_time", y="atn", ax=ax, color = "black")
sns.scatterplot(data=trigger140H, x="Date_time", y="atn", hue="trigger80pc140slope")
sns.scatterplot(data=trigger140Hfull, x="Date_time", y="atn", color = "orange")
plt.title('Hour_1 rm140')
plt.xticks(rotation=45);