# -*- coding: utf-8 -*-
"""
Created on Tue Jul 11 11:41:41 2023

@author: fabio
"""
#%%
from astropy.io import fits
import sys
import numpy as np
import scipy as sci
from matplotlib import pyplot  as plt
from gtda.time_series import Resampler, SlidingWindow, takens_embedding_optimal_parameters,    TakensEmbedding, PermutationEntropy
import os
import pandas as pd
import urllib as ull
from multiprocessing import Process


def plot_histgraph(arg1, arg2, arg3, tit):
    if( arg3 == 'p'):
        plt.plot(arg1,arg2)
    elif (arg3=='h'):
        plt.hist(arg1,arg2)
    else:
        print('options not valid')
        self.join()
    plt.title(tit)
    plt.show()


#%%
############# Functions for data analysis ###########
# see: https://www.kaggle.com/code/scratchpad/notebook33b71266b0/edit
#############################################################
#import math
#def takensEmbedding (data, delay, dimension):
#    "This function returns the Takens embedding of data with delay into dimension, delay*dimension must be < len(data)"
#    if delay*dimension > len(data):
#        raise NameError('Delay times dimension exceed length of data!')    
#    embeddedData = np.array([data[0:len(data)-delay*dimension]])
#     for i in range(1, dimension):
#         embeddedData = np.append(embeddedData, [data[i*delay:len(data) - delay*(dimension - i)]], axis=0)
#     return embeddedData;

# def mutualInformation(data, delay, nBins):
#     "This function calculates the mutual information given the delay"
#     I = 0;
#     xmax = max(data);
#     xmin = min(data);
#     delayData = data[delay:len(data)];
#     shortData = data[0:len(data)-delay];
#     sizeBin = abs(xmax - xmin) / nBins;
#     #the use of dictionaries makes the process a bit faster
#     probInBin = {};
#     conditionBin = {};
#     conditionDelayBin = {};
#     for h in range(0,nBins):
#         if h not in probInBin:
#             conditionBin.update({h : (shortData >= (xmin + h*sizeBin)) & (shortData < (xmin + (h+1)*sizeBin))})
#             probInBin.update({h : len(shortData[conditionBin[h]]) / len(shortData)});
#         for k in range(0,nBins):
#             if k not in probInBin:
#                 conditionBin.update({k : (shortData >= (xmin + k*sizeBin)) & (shortData < (xmin + (k+1)*sizeBin))});
#                 probInBin.update({k : len(shortData[conditionBin[k]]) / len(shortData)});
#             if k not in conditionDelayBin:
#                 conditionDelayBin.update({k : (delayData >= (xmin + k*sizeBin)) & (delayData < (xmin + (k+1)*sizeBin))});
#             Phk = len(shortData[conditionBin[h] & conditionDelayBin[k]]) / len(shortData);
#             if Phk != 0 and probInBin[h] != 0 and probInBin[k] != 0:
#                 I -= Phk * math.log( Phk / (probInBin[h] * probInBin[k]));
#     return I;

col=4
rows=15
lag=[[0] * col for i in range(rows)]

emb=[[0]*col for i in range(rows)] 


#%%
#last_tab='grb_table_1694558565.txt'
colnames=['GRB','Time [UT]', 'Trigger Number', 
          'BAT T90 [sec]','BAT Fluence (15-150 keV) [10^-7 erg/cm^2]',
          'BAT Fluence 90% Error (15-150 keV) [10^-7 erg/cm^2]','Redshift']
#['GRB','Time [UT]','Trigger Number', 'RA','Dec','T90', 'Fluence','Mag', 'RedSh']
base_dir=os.getcwd()
#%%
if (len(sys.argv) ==2) and (sys.argv[1]=='N'): #Network
    base='https://swift.gsfc.nasa.gov/'

############### produce last table at https://swift.gsfc.nasa.gov/archive/grb_table/fullview/ 
    url_base='https://swift.gsfc.nasa.gov/archive/grb_table/fullview/'
    ff=ull.request.urlopen(url_base)
    url_content=ff.read()
    testo=url_content.decode(encoding='utf-8', errors='ignore')
    index=testo.find('href="/archive/grb_table/tmp/grb_table')
    prelen=len('href="/archive/grb_table/tmp/')
    last_tab=testo[index+prelen:index+prelen+24]
    print(last_tab)
    ff.close()

    url=base+'archive/grb_table/tmp/'+last_tab

    tab=pd.read_csv(url,sep='\t',header=0,encoding='Latin-1')

#ColNames = tab.columns
    last_table=base_dir+'/GrbTableNew.xlsx'
    tab.to_excel(last_table,columns=colnames,index=False)
else:
    last_table=base_dir+'/GrbTableNew.xlsx'
    tab=pd.read_excel(last_table,header=0)
#####insert here candidates file
evts=pd.read_csv('events_1.txt',header=0)
#evts=pd.read_csv('events_2.txt',header=0)
ColNames = tab.columns
#%%
# histogram of log10(T90)
T90vect=[]
for (i, row) in tab.iterrows():
    print(row)
    event=tab.loc[tab['GRB']==row.GRB]
    T90=event['BAT T90 [sec]'].iloc[0]
    try:
        T90nr=float(T90)
    except:
        pass
    T90vect.append(T90nr)
    title="T90 Graph"
p = Process(target=plot_histgraph, args=[np.log10(T90vect),int(np.sqrt(1800)), 'h', title])
p.start()
#plt.hist(np.log10(T90vect),int(np.sqrt(1800)))
#plt.show()
del(event, T90, T90nr, T90vect)

########## MAIN LOOP  #############
ranget_s=200
T90vect=[]
baseDB='https://heasarc.gsfc.nasa.gov/'
basesearch=baseDB+'FTP/swift/data/obs'
j=0
for (i, row) in evts.iterrows():
    GRBname=row.GRB
    print(str(i)+' '+GRBname)
    
    event=tab.loc[tab['GRB']==row.GRB]
    T90=event['BAT T90 [sec]'].iloc[0]
    T90nr=float(T90)
    T90vect.append(T90nr)
    ranget=ranget_s/64e-3
#    ranget=int(10*T90nr/64e-3) #T90 is in seconds, sampling is 64ms;this gets a nr of samples 
    rangetn=int(ranget/4)
    rangetp=int(3*ranget/4)
    Trg=event['Trigger Number'].iloc[0]
    year=GRBname[0:2]
    month=GRBname[2:4]
    day=GRBname[4:6]
    if(len(GRBname)>6):
        seq=GRBname[6:7]
    urlname= basesearch+'/20'+year+'_'+month+'/'
    filename1='00'+Trg+'000'
    f=ull.request.urlopen(urlname)
    urlstring=f.read()
    indice=urlstring.decode().find('href=\"'+filename1) 
    f.close()
    if(indice>0):
        folder=urlname+filename1+'/bat/rate/'
        fitsfilename='sw'+filename1+'brtms.lc'
        #save file locally - function soon deprecated, 
#        ull.request.urlretrieve(folder+fitsfilename+'.gz', fitsfilename+'.gz')
        file=fits.open(folder+fitsfilename+'.gz')
        info=file.info()
        trgtime=file['PRIMARY'].header['TRIGTIME']
        data=file['RATE'].data
        t=data['TIME']-trgtime
        sample_rate=1/(t[2]-t[1])
        print(f"sample_rate {sample_rate}")
        counts=data['COUNTS'] #4 counters 
        nz=np.where(t>0)
        n=nz[0][0]
        tit='GRB'+GRBname
#        p1=Process(target=plot_histgraph, args=[t[n-rangetn:n+rangetp],counts[n-rangetn:n+rangetp][:],'p', tit])
#        p1.start()
##        plt.plot(t[n-rangetn:n+rangetp],counts[n-rangetn:n+rangetp][:])
##        plt.title('GRB'+GRBname)
       # plt.legend()
       # plt.show()
        file.close()
#%%
        analyse=True
        if analyse == True:
            max_time_delay = 3
            max_embedding_dimension = 10
            stride = 1
            for k in [0,1,2,3]: #k number of detector i=evt nr,
                optimal_time_delay, optimal_embedding_dimension = takens_embedding_optimal_parameters(counts[n-rangetn:n+rangetp, k], max_time_delay, max_embedding_dimension,  stride=stride)
                lag[i][k]= optimal_time_delay
                emb[i][k]= optimal_embedding_dimension
                print(f"detector {k}")
                print(f"Optimal embedding time delay based on mutual information: {optimal_time_delay}")
                print(f"Optimal embedding dimension based on false nearest neighbors: {optimal_embedding_dimension}")
        
nemb=np.array(emb)                
plt.hist(nemb)
for i in range(12):
    print(T90vect[i],nemb[;,0])