# -*- 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 itertools

import ephem

import pygam

import re

import numpy as np
from scipy import stats

import d3heatmap

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")
#trim to time frame minus crab
SESdata = SESdata.loc['2023-05-08':'2023-08-25']  
#remove infinite values
SESdata= SESdata[SESdata['atn']!= np.inf]

daily_SESatn = SESdata['atn'].resample('1D').mean()

#24 hour rolling mean is 24 * 3 per hour = rm 72 samples
SESdata['rm_24H'] = SESdata['atn'].rolling(72, min_periods=24,center=True).mean()

#3D rolling mean is 72*3 per hour = rm 216 samples
SESdata['rm_3D'] = SESdata['atn'].rolling(216, min_periods=72,center=True).mean()

#5D rolling mean
SESdata['rm_5D'] = SESdata['atn'].rolling(360, min_periods=100,center=True).mean()

#7D rolling mean
SESdata['rm_7D'] = SESdata['atn'].rolling(504, min_periods=150,center=True).mean()

#14D rolling mean
SESdata['rm_14D'] = SESdata['atn'].rolling(1008, min_periods=272,center=True).mean()

rawmean = np.nanmean(SESdata['atn'])
print(rawmean)
rawSD = np.nanstd(SESdata['atn'])
print(rawSD)
meanplus2SD = rawmean+rawSD+rawSD
print(meanplus2SD)

rm_24Hmean = np.nanmean(SESdata['rm_24H'])
rm_24HSD = np.nanstd(SESdata['rm_24H'])
rm_24meanplus2SD = rm_24Hmean+rm_24HSD+rm_24HSD

rm_3Dmean = np.nanmean(SESdata['rm_3D'])
rm_3DSD = np.nanstd(SESdata['rm_3D'])
rm_3Dmeanplus2SD = rm_3Dmean+rm_3DSD+rm_3DSD

rm_5Dmean = np.nanmean(SESdata['rm_5D'])
rm_5DSD = np.nanstd(SESdata['rm_5D'])
rm_5Dmeanplus2SD = rm_5Dmean+rm_5DSD+rm_5DSD

rm_7Dmean = np.nanmean(SESdata['rm_7D'])
rm_7DSD = np.nanstd(SESdata['rm_7D'])
rm_7Dmeanplus2SD = rm_7Dmean+rm_7DSD+rm_7DSD

rm_14Dmean = np.nanmean(SESdata['rm_14D'])
rm_14DSD = np.nanstd(SESdata['rm_14D'])
rm_14Dmeanplus2SD = rm_14Dmean+rm_14DSD+rm_14DSD


##NOTE- there are gaps in the time series above. this doesn't account for that since it's a quick look
SESdata['atn'].plot(label='raw', color="lightgray")
#daily_SESatn['atn'].plot(label='1D') broken and can't figure out why
SESdata['rm_24H'].plot(label='rm24H')
SESdata['rm_3D'].plot(label='rm3D')
SESdata['rm_5D'].plot(label='rm5D') 
SESdata['rm_7D'].plot(label='rm7D') 
SESdata['rm_14D'].plot(label='rm14D')
plt.title('various attenuance rolling means')
plt.hlines(y=meanplus2SD, xmin='2023-05-08', xmax='2023-08-25', color='darkgray', label='raw mean + 2SD')
plt.hlines(y=rm_24meanplus2SD, xmin='2023-05-08', xmax='2023-08-25', color='cornflowerblue', label='24rm mean + 2SD')
plt.hlines(y=rm_3Dmeanplus2SD, xmin='2023-05-08', xmax='2023-08-25', color='orange', label='3Drm mean + 2SD')
plt.hlines(y=rm_5Dmeanplus2SD, xmin='2023-05-08', xmax='2023-08-25', color='mediumseagreen', label='5Drm mean + 2SD')
plt.hlines(y=rm_7Dmeanplus2SD, xmin='2023-05-08', xmax='2023-08-25', color='lightcoral', label='7Drm mean + 2SD')
plt.hlines(y=rm_14Dmeanplus2SD, xmin='2023-05-08', xmax='2023-08-25', color='thistle', label='14Drm mean + 2SD')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))


SESdata['trigger raw']=np.where(SESdata['atn'] >= meanplus2SD, 'Pulse', 'non_p')
SESdata['trigger 24rm']=np.where(SESdata['rm_24H'] >= rm_24meanplus2SD, 'Pulse', 'non_p')
SESdata['trigger rm_3D']=np.where(SESdata['rm_3D'] >= rm_3Dmeanplus2SD, 'Pulse', 'non_p')
SESdata['trigger rm_5D']=np.where(SESdata['rm_5D'] >= rm_5Dmeanplus2SD, 'Pulse', 'non_p')
SESdata['trigger rm_7D']=np.where(SESdata['rm_7D'] >= rm_7Dmeanplus2SD, 'Pulse', 'non_p')
SESdata['trigger rm_14D']=np.where(SESdata['rm_14D'] >= rm_14Dmeanplus2SD, 'Pulse', 'non_p')

#Not identify these as respective pulse periods
SESdata.to_csv("SESdata.csv")

#how many methods picked up pulse for each particular sample
#SESdataPulses= SESdata[SESdata.apply(lambda row: row.astype(str).str.contains('Pulse', case=False).any(), axis=1)]
SESdata['sum'] = SESdata.apply(lambda x: x.str.contains('Pulse'), axis=1).sum(axis=1)
SESdata['sum'].plot()

#Anything with three or 4 hits takes place in August. Let's zoom in on that period
SESdataAUG = SESdata.loc['2023-08-01':'2023-08-25']  
#SESdataAUG['sum'].plot()
#SESdataAUG['atn'].plot(label='raw', color="lightgray", alpha=0.8)
SESdataAUG['rm_24H'].plot(label='rm24H', alpha=0.7, color='cornflowerblue')
SESdataAUG['rm_3D'].plot(label='rm3D', alpha=0.7, color='orange')
SESdataAUG['rm_5D'].plot(label='rm5D', alpha=0.7, color='mediumseagreen') 
SESdataAUG['rm_7D'].plot(label='rm7D', alpha=0.7, color='lightcoral') 
SESdataAUG['rm_14D'].plot(label='rm14D', alpha=0.7, color='thistle')
#plt.hlines(y=meanplus2SD, xmin='2023-08-01', xmax='2023-08-25', color='darkgray', label='raw mean + 2SD', alpha=0.5)
plt.hlines(y=rm_24meanplus2SD, xmin='2023-08-01', xmax='2023-08-25', color='cornflowerblue', label='24rm mean + 2SD', alpha=0.5)
plt.hlines(y=rm_3Dmeanplus2SD, xmin='2023-08-01', xmax='2023-08-25', color='orange', label='3Drm mean + 2SD', alpha=0.5)
plt.hlines(y=rm_5Dmeanplus2SD, xmin='2023-08-01', xmax='2023-08-25', color='mediumseagreen', label='5Drm mean + 2SD', alpha=0.5)
plt.hlines(y=rm_7Dmeanplus2SD, xmin='2023-08-01', xmax='2023-08-25', color='lightcoral', label='7Drm mean + 2SD', alpha=0.5)
plt.hlines(y=rm_14Dmeanplus2SD, xmin='2023-08-01', xmax='2023-08-25', color='thistle', label='14Drm mean + 2SD', alpha=0.5)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#Now see how well we can pick up on these with a rolling mean that can only look backward, using a slope-based trigger
SESdata["rm150"] = SESdata['atn'].rolling(150, min_periods=50).mean()
SESdata['diff_rm150'] = SESdata['rm150'].diff()

#find max slope after this
SESslope150TRIM = SESdata.iloc[151:, ]
trigger150slope = max(SESslope150TRIM['diff_rm150'])
print(trigger150slope)
#0.013205730433333396
SESslope150TRIM['diff_rm150'].plot(label='diff_rm150', alpha=0.7, color='cornflowerblue')

#top five slopes are all above 0.011: 


daily_SESatn = pd.DataFrame(daily_SESatn)

#rfilter to complete cases
daily_SESatn= daily_SESatn.loc[daily_SESatn.notna().all(axis='columns')]

SESdataAUG['atn'].plot()
daily_SESatn['atn'].plot()
plt.title('all and daily atn')


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

#trim to time frame minus crab
SESdata = SESdata.loc['2023-05-08':'2023-08-25']  
daily_SESatn = daily_SESatn.loc['2023-05-08':'2023-08-25']  
#Now trim it to once an hour. Make three separate datasets 
SESonceanhour_Adf = SESdata.iloc[0::3, :] #295 rolling mean
SESonceanhour_Bdf = SESdata.iloc[1::3, :] #326 rolling mean
SESonceanhour_Cdf = SESdata.iloc[2::3, :] #247 rolling mean

#The plot below shows the raw data (full MARS time series),
 #with moving averages of 24 hours, 3 days, 5 days, 7 days, and 14 days. 
 #Horizontal lines identify the threshold of each respective time series’ mean + 2 SD.