# -*- 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 numpy.ma as ma
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
Station_M_POCdata = pd.read_csv('G:/Analyses/triggering/Station_M_POCdata.csv')

Station_M_POCdata.columns

#Index(['Collect_date', 'Mass_flux', 'POC_flux'], dtype='object')


Station_M_POCdata['Date']= pd.to_datetime(Station_M_POCdata['Collect_date'], format='%m/%d/%Y')


sns.scatterplot(data=Station_M_POCdata, x="Date", y="POC_flux")
plt.xticks(rotation=45);

stdStaM = Station_M_POCdata.std()
# Mass_flux                       93.447106
# POC_flux                         7.444236
# Date         3142 days 22:48:45.324293440
# dtype: object

stdevPOC = stdStaM[1]
stdevMass = stdStaM[0]


meanStaM = Station_M_POCdata.mean()
# Mass_flux    112.859681
# POC_flux       8.351651
# dtype: float64

meanPOC = meanStaM[1]
meanMass = meanStaM[0]

PulsethresholdPOC = meanPOC+stdevPOC+stdevPOC
print(PulsethresholdPOC)

#make pulse column
Station_M_POCdata["SDtrigger"] = np.where(Station_M_POCdata['POC_flux'] >= PulsethresholdPOC, 'TriggerSD', 'WaitSlope')


sns.scatterplot(data=Station_M_POCdata, x="Date", y="POC_flux", hue="SDtrigger")
plt.xticks(rotation=45);
plt.show()

len(Station_M_POCdata)

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

count=0
R_POCcount=[]
R_POCmeancount=[]
n_count=[]
dfcount=[]
SErcount=[]
CIhighr=[]
trigger=[]
trigger_count = []
wait_count = []
triggerlist=[]

for x in range(90, 800):
        data =Station_M_POCdata
        variablename = str(x)
        data[variablename] = data['POC_flux'].rolling(x, min_periods=25).mean()
        R_POCcount = data['POC_flux'].rolling(x, min_periods=25).mean()
        R_POCmeancount = R_POCcount.mean()
        n_count = len(R_POCcount)
        dfcount = n_count-1
        SErcount = stats.sem(R_POCcount.dropna())
        CIlevel = 0.9999
        CI99POCrcount = stats.t.interval(CIlevel, dfcount, R_POCmeancount, scale = SErcount)
        CIhighr = float(CI99POCrcount[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']+90
#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()

sns.lineplot(data=triggerlistdf, x="rollingmean_count", y="number_of_triggers")
plt.xticks(rotation=45);
plt.show()

#There's a dip around 300 days. find the local min there
#301



#make a new time series of this rolling mean
Station_M_POCdata['rolling_POC_flux_301'] = Station_M_POCdata['POC_flux'].rolling(rollingmeanUSE, min_periods=30).mean()

Station_M_POCdata['rolling_POC_flux_301'].astype(float).pct_change(1)

#subset to remove first 301 rows (to eliminate artifacts of using rolling mean)
Station_M_POCdataTRIM = Station_M_POCdata.iloc[301:, ]
Station_M_POCdataTRIM = Station_M_POCdataTRIM[["Date", "rolling_POC_flux_301", "POC_flux"]]

Station_M_POCdataTRIM['diff'] = Station_M_POCdataTRIM['rolling_POC_flux_301'].diff()


#find max slope after this
trigger301slope = Station_M_POCdataTRIM['diff'].max()
print(trigger301slope)
#0.5152874753894743

#make that new max the trigger
Station_M_POCdataTRIM["trigger301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger301slope, '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
trigger80pc301slope = 0.8 * trigger301slope
print(trigger80pc301slope)
#0.4122299803115794

trigger75pc301slope = 0.75 * trigger301slope
print(trigger75pc301slope)
#0.3864656065421057


Station_M_POCdataTRIM["trigger301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger301slope, 'TriggerSlope', 'WaitSlope')

Station_M_POCdataTRIM["trigger80pc301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger80pc301slope, 'TriggerSlope', 'WaitSlope')

Station_M_POCdataTRIM["trigger75pc301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger75pc301slope, 'TriggerSlope', 'WaitSlope')

# sns.lineplot(data=Station_M_POCdataTRIM, x="Date", y="diff", hue="trigger301slope")
# plt.xticks(rotation=45);
# plt.show()


# sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="diff", hue="trigger80pc301slope")
# plt.xticks(rotation=45);
# plt.show()

# sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", hue="trigger80pc301slope")
# plt.xticks(rotation=45);
# plt.show()

# sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", hue="trigger75pc301slope")
# plt.xticks(rotation=45);
# plt.show()

# #this is hard to see so let's check from 2010 onward
# POC2010onward =(Station_M_POCdataTRIM['Date'] > '2010-01-01') 

# #filter so only rows with trigger command
# Station_M_POCdataTRIMtrigger = Station_M_POCdataTRIM.loc[Station_M_POCdataTRIM['trigger75pc301slope'] == 'TriggerSlope']

# fig, ax = plt.subplots()
# sns.lineplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", ax=ax)
# sns.scatterplot(data=Station_M_POCdataTRIMtrigger, x="Date", y="POC_flux", ax=ax, color = 'orange')
# plt.title('trigger75pc301slope')

trigger50pc301slope = 0.50 * trigger301slope
print(trigger50pc301slope)
#0.3864655065421057

Station_M_POCdataTRIM["trigger50pc301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger50pc301slope, 'TriggerSlope', 'WaitSlope')

sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="diff", hue="trigger50pc301slope")
plt.xticks(rotation=45);
plt.show()

Station_M_POCdataTRIMtrigger50 = Station_M_POCdataTRIM.loc[Station_M_POCdataTRIM['trigger50pc301slope'] == 'TriggerSlope']

fig, ax = plt.subplots()
sns.lineplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", ax=ax)
sns.scatterplot(data=Station_M_POCdataTRIMtrigger50 , x="Date", y="POC_flux", ax=ax, color = 'orange')
plt.title('trigger50pc301slope')





trigger60pc301slope = 0.60 * trigger301slope
print(trigger60pc301slope)
#0.3864656065421057

Station_M_POCdataTRIM["trigger60pc301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger60pc301slope, 'TriggerSlope', 'WaitSlope')

sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="diff", hue="trigger80pc301slope")
plt.xticks(rotation=45);
plt.show()

Station_M_POCdataTRIMtrigger60 = Station_M_POCdataTRIM.loc[Station_M_POCdataTRIM['trigger60pc301slope'] == 'TriggerSlope']

fig, ax = plt.subplots()
sns.lineplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", ax=ax)
sns.scatterplot(data=Station_M_POCdataTRIMtrigger60 , x="Date", y="POC_flux", ax=ax, color = 'orange')
plt.title('trigger60pc301slope')


trigger70pc301slope = 0.70 * trigger301slope
print(trigger70pc301slope)
#0.3864657065421057

Station_M_POCdataTRIM["trigger70pc301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger70pc301slope, 'TriggerSlope', 'WaitSlope')

sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="diff", hue="trigger80pc301slope")
plt.xticks(rotation=45);
plt.show()

Station_M_POCdataTRIMtrigger70 = Station_M_POCdataTRIM.loc[Station_M_POCdataTRIM['trigger70pc301slope'] == 'TriggerSlope']

fig, ax = plt.subplots()
sns.lineplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", ax=ax)
sns.scatterplot(data=Station_M_POCdataTRIMtrigger70 , x="Date", y="POC_flux", ax=ax, color = 'orange')
plt.title('trigger70pc301slope')



trigger80pc301slope = 0.80 * trigger301slope
print(trigger80pc301slope)
#0.3864658065421057

Station_M_POCdataTRIM["trigger80pc301slope"] = np.where(Station_M_POCdataTRIM['diff'] >= trigger80pc301slope, 'TriggerSlope', 'WaitSlope')

sns.scatterplot(data=Station_M_POCdataTRIM, x="Date", y="diff", hue="trigger80pc301slope")
plt.xticks(rotation=45);
plt.show()

Station_M_POCdataTRIMtrigger80 = Station_M_POCdataTRIM.loc[Station_M_POCdataTRIM['trigger80pc301slope'] == 'TriggerSlope']

fig, ax = plt.subplots()
sns.lineplot(data=Station_M_POCdataTRIM, x="Date", y="POC_flux", ax=ax)
sns.scatterplot(data=Station_M_POCdataTRIMtrigger80 , x="Date", y="POC_flux", ax=ax, color = 'orange')
plt.title('trigger80pc301slope')


Station_M_POCdataTRIM.to_csv("Station_M_POCdataTRIM.csv")


fig, ax = plt.subplots(2,2)
sns.lineplot(ax=ax[0, 0], data=Station_M_POCdataTRIM, x="Date", y="POC_flux")
sns.scatterplot(ax=ax[0, 0], data=Station_M_POCdataTRIMtrigger50 , x="Date", y="POC_flux", color = 'orange')
ax[0][0].set_title('trigger50pc301slope')
ax[0][0].tick_params(axis='x', rotation=45)
ax[0][0].axhline(PulsethresholdPOC, ls='--', linewidth=1, color='red')


sns.lineplot(ax=ax[0, 1], data=Station_M_POCdataTRIM, x="Date", y="POC_flux")
sns.scatterplot(ax=ax[0, 1], data=Station_M_POCdataTRIMtrigger60 , x="Date", y="POC_flux", color = 'orange')
ax[0][1].set_title('trigger60pc301slope')
ax[0][1].tick_params(axis='x', rotation=45)
ax[0][1].axhline(PulsethresholdPOC, ls='--', linewidth=1, color='red')


sns.lineplot(ax=ax[1, 0], data=Station_M_POCdataTRIM, x="Date", y="POC_flux")
sns.scatterplot(ax=ax[1, 0], data=Station_M_POCdataTRIMtrigger70 , x="Date", y="POC_flux", color = 'orange')
ax[1][0].set_title('trigger70pc301slope')
ax[1][0].tick_params(axis='x', rotation=45)
ax[1][0].axhline(PulsethresholdPOC, ls='--', linewidth=1, color='red')

sns.lineplot(ax=ax[1, 1], data=Station_M_POCdataTRIM, x="Date", y="POC_flux")
sns.scatterplot(ax=ax[1, 1], data=Station_M_POCdataTRIMtrigger80 , x="Date", y="POC_flux", color = 'orange')
ax[1][1].set_title('trigger80pc301slope')
ax[1][1].tick_params(axis='x', rotation=45)
plt.tight_layout()
ax[1][1].axhline(PulsethresholdPOC, ls='--', linewidth=1, color='red')