# -*- 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")

#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 275 to 1500 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=[]

for x in range(274, 1500):
        data =Main_Hourly_SES_MARS
        variablename = str(x)
        data[variablename] = data['atn'].rolling(x, min_periods=50).mean()
        R_atncount = data['atn'].rolling(x, min_periods=50).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']+274
#triggerlistdf= triggerlistdf.rename(columns={"index": "rollingmean_count"})
triggerlistdf.to_csv("triggerlistdf.csv")


#find the rolling mean with the fewest triggers
rowwithlottrigger = triggerlistdf[triggerlistdf.number_of_triggers==triggerlistdf.number_of_triggers.min()]
rollingmeanUSE = rowwithlottrigger.rollingmean_count.min()
print(rollingmeanUSE)

#make a new time series of this rolling mean
Main_Hourly_SES_MARS['rolling_atn_1153'] = Main_Hourly_SES_MARS['atn'].rolling(rollingmeanUSE, min_periods=200).mean()

#Now check the slope of this time series
#subset columns
SESMARS1153atn = Main_Hourly_SES_MARS[["atn", "rolling_atn_1153"]]

#find proxy for slope (percent change from one sample to next)
SESMARS1153atn['Percentchange'] = SESMARS1153atn['rolling_atn_1153'].astype(float).pct_change(1)

#subset to remove first 1153 rows (to eliminate artifacts of using rolling mean)
SESslope1153TRIM = SESMARS1153atn.iloc[1153:, ]

#find max slope after this
trigger1153slope = max(SESslope1153TRIM['Percentchange'])
print(trigger1153slope)
#0.006808192293436877

#make that new max the trigger
SESslope1153TRIM["trigger1153slope"] = np.where(SESslope1153TRIM['Percentchange'] >= trigger1153slope, '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
trigger80pc1153slope = 0.8 * trigger1153slope
print(trigger80pc1153slope)
#0.005446553834749502

#list which samples are trigger vs wait
SESslope1153TRIM["trigger80pc1153slope"] = np.where(SESslope1153TRIM['Percentchange'] >= trigger80pc1153slope, 'TriggerSlope', 'WaitSlope')

#plot to verify.
sns.scatterplot(data=SESslope1153TRIM, x="datetime", y="atn", hue="trigger80pc1153slope")
plt.xticks(rotation=45);
plt.show()