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





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

#remove infinite values
SESdata= SESdata[SESdata['atn']!= np.inf]

SESdata.dtypes
SEStrim = SESdata.drop(['Date', 'Time'], axis=1)

SEStrim["atn"].plot()

#find hourly average of SES time series
hourly_SESatn = SEStrim['atn'].resample('1H').mean()

#transpose
hourly_SESatndf = pd.DataFrame([hourly_SESatn]).T

#fix weird time zone stuff
for col in hourly_SESatndf.select_dtypes(include=['datetime64[ns, UTC]']).columns:
    hourly_SESatndf[col] = hourly_SESatndf[col].apply(lambda x: x.tz_localize(None))
    
#Trim to period w prior to crab clog
trimhourly = hourly_SESatndf.loc['2023-06-27':'2023-08-14']  
trimhourly.plot()
plt.title('hourly, trimmed to no crab period')
plt.ylabel('attenuance')

#fill na with interpolation
timehourlyinterpolate = trimhourly.interpolate(method='linear')


#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(timehourlyinterpolate['atn'])
plt.subplot(2, 1, 1)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [units]')
plt.ylabel('PSD [V**2/units]')
plt.title('hourly, trimmed to no crab period, welch')
plt.show()

#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/

plot_acf(timehourlyinterpolate['atn'], lags=60)
plt.title('hourly, trimmed to no crab period, autocorrelation')

#adjusting N doesn't make a difference
#fft foesn't make a difference
#drop na no difference

plot_acf(timehourlyinterpolate['atn'], lags=60, missing = 'drop')
plt.title('hourly, trimmed to no crab period, autocorrelation drop na')

#adjusting N doesn't make a difference

#These are the acf results for hourly averages
acfarray = acf(timehourlyinterpolate['atn'], adjusted=False, 
                   nlags=60, qstat=True, fft=False, alpha=0.05, 
                   bartlett_confint=True, missing='none')

results = pd.DataFrame()
results["acf"] = acf(timehourlyinterpolate['atn'], nlags=60)

results.to_csv("C:/Python_temp/acfhourly.csv")


plot_pacf(timehourlyinterpolate['atn'], lags=60)
plt.title('hourly, trimmed to no crab period, partial autocorrelation')



#plot by hours of the day
timehourlyinterpolate.index.dtype #dtype('<M8[ns]')
timehourlyinterpolate['datetime'] = timehourlyinterpolate.index.strftime('%Y-%m-%d %H:%M:%S')
timehourlyinterpolate['datetime']  = pd.to_datetime(timehourlyinterpolate['datetime'] , format='%Y-%m-%d %H:%M:%S')
timehourlyinterpolate.datetime.dtype #dtype('O')
timehourlyinterpolate['hour'] = timehourlyinterpolate['datetime'].dt.hour


boxplot = timehourlyinterpolate.boxplot(by='hour')
plt.savefig('bot plots atn by hour.eps',format='eps')



sns.relplot(data=timehourlyinterpolate, x="hour", y="atn", kind="line")
plt.savefig('atn by 24 hour.eps',format='eps')


##see if different hours of the day have different cycles

midnight_atn = (timehourlyinterpolate[timehourlyinterpolate.hour == 0])

resultsmidnight = pd.DataFrame()
resultsmidnight["acf"] = acf(timehourlyinterpolate['atn'])


plot_acf(midnight_atn['atn'], missing = 'drop')
plt.title('midnight, trimmed to no crab period, autocorrelation drop na')


##see if different hours of the day have different cycles
a = pd.DataFrame()
b = pd.DataFrame()
c = []
plot_title = [] 

for i in range(24):
    a = (timehourlyinterpolate[timehourlyinterpolate.hour == i])
    b = acf(a['atn'])
    c.append(b)
    plot_title = i, 'trimmed to no crab period, autocorrelation drop na'
    plot_acf(a['atn'], missing = 'drop')
    plt.title(plot_title)
    plt.savefig(str(i)+'trimmed to no crab period, autocorrelation drop na.pdf')
    
    #spectral density by hour
    
for i in range(24):
    a = (timehourlyinterpolate[timehourlyinterpolate.hour == i])
    b = acf(a['atn'])
    fs = 1
    plot_title = i, '_spectral density'
    f, Pxx_den = signal.welch(a['atn'])
    plt.semilogy(f, Pxx_den)
    plt.ylim([0.5e-3, 1])
    plt.xlabel('frequency [units]')
    plt.ylabel('PSD [V**2/units]')
    plt.title(plot_title)
    plt.savefig(str(i)+'spectral density.pdf')
    plt.close()


#plot all against attenuation time series

#moon phase illumination
timehourlyinterpolate['year'] = timehourlyinterpolate['datetime'].dt.year
timehourlyinterpolate['month'] = timehourlyinterpolate['datetime'].dt.month
timehourlyinterpolate['day'] = timehourlyinterpolate['datetime'].dt.day
timehourlyinterpolate['simpleEdatetime'] = timehourlyinterpolate['year'].map(str) +'/'+ timehourlyinterpolate['month'].map(str) + '/' + timehourlyinterpolate['day'].map(str)

#test finding moonphase
#date = '1974/11/11'
m = ephem.Moon('1974/11/11')
print(m.moon_phase)

m = ephem.Moon('2023/06/27')
print(m.moon_phase)

date = timehourlyinterpolate['datetime'].iloc[0]
m = ephem.Moon(date)
print(m.moon_phase)


#how many days in sample 8/14/2023 - 6/27/2023
n = timehourlyinterpolate.nunique(axis=0)
print(n) #simpleEdatetime

m = ephem.Moon()
m.compute('2023/06/27')
j = ephem.Date('2023/06/27')
for add in range(0, 49):
    print (ephem.date(j + add))
    print(m.moon_phase)



b = pd.DataFrame()
c = []
e=[]
l=[]
#needs padding of the date for interpolating the moon phase
j = ephem.Date('2023/06/23')
for add in range(0, 55):
    d = (ephem.date(j + add))
    m = ephem.Moon(d)
    moonphase = m.moon_phase
    c.append(moonphase)
    print(ephem.date(d + ephem.hour))
    local = ephem.localtime(d)
    e.append(d)
    l.append(local)
    
moonphasedata = pd.DataFrame(
    {'local_time': l,
     'UTCtime': e,
     'moonphase': c
    })


COPYtimehourlyinterpolate = timehourlyinterpolate.copy()

#Needs to get into same time zone. this says UTC but it's actually local
COPYtimehourlyinterpolate['local_time'] = COPYtimehourlyinterpolate['datetime']
COPYtimehourlyinterpolate['UTC'] = pd.to_datetime(COPYtimehourlyinterpolate['local_time'], utc = True)
moonphasedata['UTC'] = pd.to_datetime(moonphasedata['local_time'], utc = True)
atnmoonphase = pd.concat([COPYtimehourlyinterpolate, moonphasedata])
#this does funky things to the dates and times, so aggregate by datetime
useatnmoonphase = pd.DataFrame(atnmoonphase, columns=['local_time', 'atn','moonphase'])

useatnmoonphase = useatnmoonphase.resample('1H', on='local_time').mean()

#now turn index to column
#turn index into a column
useatnmoonphase['datetime'] = useatnmoonphase.index

#now interpolate moonphase
useatnmoonphase['mooninterp'] = useatnmoonphase['moonphase'].interpolate(method='linear')
#now remove the original moonphase because of all the NAs
useatnmoonphase= useatnmoonphase.drop(columns=['moonphase'])

#plot moon phase against interpolated moon phase to make sure it's working
ax1 = useatnmoonphase.plot(x = 'datetime', y = ["mooninterp"])
ax2 = useatnmoonphase.plot.scatter(x = 'datetime', y = ['moonphase'], marker = '^', ax = ax1)
plt.xticks(rotation=45);



#now bring back the hour of the day
useatnmoonphase['hour'] = useatnmoonphase['datetime'].dt.hour
useatnmoonphase['year'] = useatnmoonphase['datetime'].dt.year
useatnmoonphase['month'] = useatnmoonphase['datetime'].dt.month
useatnmoonphase['day'] = useatnmoonphase['datetime'].dt.day
useatnmoonphase['Year-Week'] = useatnmoonphase['datetime'].dt.strftime('%Y-%U')
useatnmoonphase['Year-Week'] = useatnmoonphase['datetime'].dt.strftime('%Y-%U')

mako = sns.color_palette("rocket", 9)
sdiverging_colors = sns.color_palette("RdBu", 24)
#rfilter to complete cases
useatnmoonphase= useatnmoonphase.loc[useatnmoonphase.notna().all(axis='columns')]

#now to get the two weeks column
len(useatnmoonphase) #1176
len(useatnmoonphase) /336
twoweeks = 336 # number of times you repeat each item

yourlist = list( range(1, 337))*4
df = pd.DataFrame(yourlist)
df2 = pd.concat(df, useatnmoonphase)



sns.scatterplot(data=useatnmoonphase, x='mooninterp', y='atn', hue='month')
sns.scatterplot(data=useatnmoonphase, x='hour', y='atn', hue='mooninterp')
sns.scatterplot(data=useatnmoonphase, x='mooninterp', y='atn', hue='hour')



plt.tricontourf(useatnmoonphase["hour"], useatnmoonphase["mooninterp"], useatnmoonphase["atn"],palette = 'mako')
plt.title('Attenuance by hour and lunar illumination')
plt.xlabel('hour')
plt.ylabel('lunar illumination')
plt.show()




# instanciate figure and axes instances
fig,ax = plt.subplots(1,2,figsize=(7, 4))
a = sns.lineplot(data=useatnmoonphase, 
             x="hour", 
             y="atn", 
             hue = 'Year-Week',
             palette = mako, ci = None, ax=ax[0])
ax[0].legend_.remove()
b = sns.scatterplot(data = useatnmoonphase,
               x = 'hour', 
               y = 'atn',
               hue = 'Year-Week',
               palette = mako,ax=ax[1])
b.set(ylabel=None)
# sns.move_legend(
#     ax[1],
#     loc="right",
#     bbox_to_anchor=(-0.1, 1.1),
#     ncol=2,
#     title=None,
#     frameon=False,
# )
sns.move_legend(ax[1], "upper left", bbox_to_anchor=(1, 1))
ax[0].autoscale()
ax[1].autoscale()




sns.lineplot(data=useatnmoonphase, 
             x="hour", 
             y="atn", 
             hue = 'Year-Week',
             palette = mako, ci = None)
plt.legend(bbox_to_anchor=(1.05, 1))






#See how hourly averages compare to daily averages of atn
# useatnmoonphase['simpleEdatetime'] = useatnmoonphase['year'].map(str) +'/'+ useatnmoonphase['month'].map(str) + '/' + useatnmoonphase['day'].map(str)
# useatnmoonphase= useatnmoonphase.drop(columns=['simpleEdatetime'])
Dailyatn = useatnmoonphase.resample('1D', on='datetime').mean()
Dailyatn= Dailyatn.drop(columns=['hour'])
Dailyatn = Dailyatn.rename(columns={'atn': 'Dailyatn'})


#now plot the Daily atn vs atn for each hour of the day over time, and then as a scatter plot
#join daily and hourly datasets by year month and day
DailyHourlyatn = Dailyatn.merge(useatnmoonphase, on=['year', 'month', 'day'])

g = sns.FacetGrid(DailyHourlyatn, col='hour', hue = 'month', col_wrap=12)
plt.xticks(rotation=45)
g.map(plt.scatter, 'Dailyatn', 'atn')


plt.savefig('daily atn vs hourly atn.pdf')



# plot each our of day atn vs daily avg atn 

def annotate(DailyHourlyatn, **kws):

    # display data; see that for each Facet, hue groups are annotated separately - uncomment the following two lines
    # print(data.sex.unique())
    # display(data) 

    # get the hue group; there will be one
    g = DailyHourlyatn.month.unique()[0]

    # get the y-position from the dict
    y = yg[g]

    r, p = SP.stats.pearsonr(DailyHourlyatn['Dailyatn'], DailyHourlyatn['atn'])
    ax = plt.gca()
    ax.text(1, y, f'{g}: r={r:.2f}, p={p:.2f}')



# define a y-position for each annotation in the hue group
yg = {'May': 8, 'June': 9, 'July':10, 'Aug':11}

# plot
g = sns.lmplot(x='Dailyatn', y='atn', col='hour', data=DailyHourlyatn, hue='month', col_wrap=4, ci = None)


###try to get details on spearman rank
from scipy.stats import spearmanr



g = sns.lmplot(data=DailyHourlyatn, x='Dailyatn', y='atn', col='hour', sharey=False, sharex=False,col_wrap=4)

for ax, feature in zip(g.axes.flat, g.col_names):
    r, pvalue = spearmanr(DailyHourlyatn['atn'], DailyHourlyatn[feature])
    ax.collections[0].set_label(f'Spearman = {r:.2f}')
    ax.legend()
plt.tight_layout()
plt.show()








# d = np.transpose(c, axes=None)
    
# df2 = pd.DataFrame(d, columns=['hr0', 'hr1', 'hr2','hr3','hr4','hr5','hr6','hr7','hr8','hr9','hr10','hr11','hr12','hr13','hr14','hr15','hr16',
#                                  'hr17','hr18','hr19','hr20','hr21','hr22','hr23'])

# plot_acf(midnight_atn['atn'], missing = 'drop')
# plt.title('midnight, trimmed to no crab period, autocorrelation drop na')



# plt.scatter(df2.index.values, df2['hr2'])


#now plot by hour over time. But make hour as hours from noon
timehourlyinterpolate['hour_from_noon'] = timehourlyinterpolate['hour']-12

atn_hour_wide = timehourlyinterpolate.pivot(index="datetime", columns="hour_from_noon", values="atn")
sdiverging_colors = sns.color_palette("RdBu", 24)

sns.relplot(data=atn_hour_wide, kind="line", palette = sdiverging_colors)
plt.xticks(rotation=45)


g = sns.FacetGrid(timehourlyinterpolate, col='hour', col_wrap=12)
plt.xticks(rotation=45)
g.map(plt.plot, 'datetime', 'atn')
plt.show()

g = sns.FacetGrid(timehourlyinterpolate, col='hour', col_wrap=12)
plt.xticks(rotation=45)
g.map(plt.plot, 'datetime', 'atn')
plt.savefig('atn by hour of day over time.eps',format='eps')


#pull out each hour separately

g = sns.FacetGrid(timehourlyinterpolate, col='hour_from_noon', hue = 'hour_from_noon', col_wrap=12, palette = sdiverging_colors)
plt.xticks(rotation=45)
g.map(plt.plot, 'datetime', 'atn')
plt.show()

g = sns.FacetGrid(timehourlyinterpolate, col='hour_from_noon', hue = 'hour_from_noon', col_wrap=12, palette = sdiverging_colors)
plt.xticks(rotation=45)
g.map(plt.plot, 'datetime', 'atn')
plt.savefig('atn by hour over time.eps',format='eps')

#add moon phase
timehourlyinterpolate['moonphase'] = ephem.Moon(timehourlyinterpolate['datetime'])




######################################





#try on regular 5 min collect times
SEStrimcopy = SEStrim.copy()


#fix weird time zone stuff
for col in SEStrimcopy.select_dtypes(include=['datetime64[ns, UTC]']).columns:
    SEStrimcopy[col] = SEStrimcopy[col].apply(lambda x: x.tz_localize(None))
    
#Trim to period w prior to crab clog
trimbase = SEStrimcopy.loc['2023-06-27':'2023-08-14']  
trimbase["atn"].plot()
plt.title('5 minute collect, trim no crab period')


#fill na with interpolation
trimbaseinterpolate = trimbase.interpolate(method='linear')

#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(trimbaseinterpolate['atn'])
plt.subplot(2, 1, 1)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [units]')
plt.ylabel('PSD [V**2/units]')
plt.title('5 min collect, trimmed to no crab period, welch')
plt.show()

#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_acf

plot_acf(trimbaseinterpolate['atn'], lags=180)
plt.title('5 minute collect, trim no crab period, autocorrelation, every 20 min')
pyplot.show()


from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(trimbaseinterpolate['atn'], lags=180)
plt.title('5 minute collect, trim no crab period, partial autocorrelation, every 20 min')
pyplot.show()


#find daily average of SES time series
TwoDaily_SESatn = SEStrim['atn'].resample('2D').mean()

#transpose
TwoDaily_SESatndf = pd.DataFrame([TwoDaily_SESatn]).T

#fix weird time zone stuff
for col in TwoDaily_SESatndf.select_dtypes(include=['datetime64[ns, UTC]']).columns:
    TwoDaily_SESatndf[col] = TwoDaily_SESatndf[col].apply(lambda x: x.tz_localize(None))
    
#Trim to period w prior to crab clog
trim2Daily = TwoDaily_SESatndf.loc['2023-06-27':'2023-08-14']  
trim2Daily.plot()

#fill na with interpolation
trim2Dailyinterpolate = trim2Daily.interpolate(method='linear')

#try spectray density welch
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html
fs = 1
f, Pxx_den = signal.welch(trim2Dailyinterpolate['atn'])
plt.subplot(2, 1, 1)
plt.semilogy(f, Pxx_den)
plt.ylim([0.5e-3, 1])
plt.xlabel('frequency [units]')
plt.ylabel('PSD [V**2/units]')
plt.title('2 day, trimmed to no crab period, welch')
plt.show()
#Try autocorrelation
#from https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
from pandas import read_csv
from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_acf

plot_acf(trim2Dailyinterpolate['atn'])
plt.title('2 day, trim no crab period, autocorrelation')
pyplot.show()


from matplotlib import pyplot
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(trim2Dailyinterpolate['atn'])
plt.title('2 day, trim no crab period,partial autocorrelation')
pyplot.show()


###delete?


#find confidence intervals for each hour

# calculate the 95% confidence intervals for the mean RT across all participantID
# and provide both the upper and lower CIs in a table
results = pd. DataFrame()
results['mean_hour_atn'] = timehourlyinterpolate.groupby('hour').mean('atn')
results['sem_hour_atn'] = timehourlyinterpolate.groupby(['hour'])['atn'].sem()
results['ci_hour_atn'] = results['sem_hour_atn'] * 1.96
results['upper_ci'] = results['mean_hour_atn'] + results['ci_hour_atn']
results['lower_ci'] = results['mean_hour_atn'] - results['ci_hour_atn']

results.to_csv('SES atn by hour.csv')

results.columns
#Index(['mean_hour_atn', 'sem_hour_atn', 'ci_hour_atn', 'upper_ci', 'lower_ci'], dtype='object')
#turn index into a column
results['hour'] = results.index

#######try radial format###########
#THIS PRODUCES SOMETHING

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)

rad = np.linspace(0, 5, 100)
azm = np.linspace(0, 2 * np.pi, 100)
r, th = np.meshgrid(rad, azm)
z = (r ** 2.0) / 4.0

plt.subplot(projection="polar")

plt.pcolormesh(th, r, z)
#plt.pcolormesh(th, z, r)

plt.plot(azm, r, color='k', ls='none') 
plt.thetagrids([theta * 15 for theta in range(360//15)])


plt.show()


format###########  
    
  
    
  #more trying