# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Introduction
#
# Lorem ipsum dolor sit amet, consectetur adipiscing elit. Sed placerat pellentesque tortor at luctus. Cras varius dui odio, sit amet sodales ipsum ornare non. Mauris imperdiet interdum fermentum. Suspendisse ac nisl in dui feugiat pellentesque. In ac condimentum ligula. Nam nec arcu vel eros eleifend ultricies ut eu arcu. Phasellus dictum mauris a nunc tempor pellentesque vitae eget orci. Vestibulum gravida gravida ligula, eget rutrum dui pulvinar iaculis. Curabitur fermentum elementum purus, ac vulputate magna consectetur eu. Phasellus sodales facilisis tortor, nec iaculis ex aliquam a. Phasellus euismod justo a convallis tempus. Curabitur dignissim mi mauris.
#
# Maecenas congue ut lacus ac dapibus. Maecenas mollis, sem eget egestas pulvinar, eros augue aliquam neque, id porta neque lacus a augue. Vestibulum at pharetra velit, in facilisis mauris. Aenean sollicitudin elementum mi, eget pharetra nibh vestibulum sodales. Mauris in malesuada ipsum, vitae varius metus. Vestibulum non iaculis nibh. Pellentesque habitant morbi tristique senectus et netus et malesuada fames ac turpis egestas. Phasellus semper sodales metus id commodo. Quisque tincidunt, turpis quis imperdiet sollicitudin, ante dolor imperdiet nibh, nec iaculis risus massa non libero. Donec magna risus, dignissim eu semper ac, vestibulum quis tellus. Interdum et malesuada fames ac ante ipsum primis in faucibus. Integer eu justo non justo ullamcorper cursus eget vulputate erat. Nunc auctor quam posuere, varius dui in, accumsan mi. Donec aliquet lacus vitae orci ultricies feugiat. Proin viverra, felis vel euismod rutrum, ligula risus viverra orci, id maximus nisl urna vel neque. Integer sodales velit urna, in mattis leo ornare eu.
#
# Donec molestie eget lectus nec viverra. Nulla sed semper mauris, vitae suscipit mi. Vestibulum vel sodales magna. Vivamus laoreet vestibulum nibh, sed ornare lacus luctus id. Quisque fringilla lacus ac interdum iaculis. Aliquam accumsan nisl et libero dignissim eleifend. Sed magna enim, dictum sodales odio nec, interdum interdum tortor. Morbi sodales sem libero, in interdum diam cursus et. Quisque malesuada imperdiet sem, blandit luctus purus pharetra vel.
#
# # Data and Methods

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pywt
import atddm
from constants import AIRPORTS, COLORS, TZONES, CODES, BEGDT, ENDDT
import seaborn as sns
from math import log2, floor
from matplotlib import colors
from scipy import signal
from statsmodels.robust import mad

BEGDT = pd.Timestamp(BEGDT)
ENDDT = pd.Timestamp(ENDDT)
INTERVAL = 10          # in minutes
TIMESTEP = 60*INTERVAL # in seconds
NLEVEL = 8
WTYPE = 'db5'
# CMAP = colors.ListedColormap(['red', 'darkred',
#                               'coral', 'orangered',
#                               'goldenrod', 'darkgoldenrod',
#                               'limegreen', 'darkgreen',
#                               'lightseagreen', 'seagreen',
#                               'steelblue', 'cadetblue',
#                               'blue', 'navy',
#                               'darkviolet', 'purple'
#                                ])
CMAP = 'tab20'
times = np.array([1800, 3600, 10800, 21600, 43200, 86400, 604800]) # in seconds
XTICKS = 1./times
CUTFRQ = XTICKS[2]
XTKLBS = ['.5 h', '1 h', '3 h', '6 h', '12 h', '24 h', '1 w']

# +
dd = atddm.load(subset=CODES)
m3_bin = {}
# m1_bin = {}
CODES.sort()

for code in CODES:
    indx = pd.date_range(start=BEGDT, end=ENDDT,
                         freq=str(INTERVAL)+'min', tz=TZONES[code])
    m3_bin[code] = atddm.binarrivals(dd[code].M3_FL240,
                                     interval=INTERVAL,
                                     tz=TZONES[code])[indx].fillna(0)
#    m1_bin[code] = atddm.binarrivals(dd[code].M1_FL240,
#                                     interval=INTERVAL,
#                                     tz=TZONES[code])[indx].fillna(0)
# -

# ## Wavelet analysis

# +
def findm(length, n=NLEVEL):
    return floor(length/2**n)


def trimmedindex(serie, nlev=NLEVEL):
    m = findm(len(serie), nlev)
    lenmax = m * 2**nlev
    return serie.index[:lenmax]


def wvlt_analysis(serie, wtype=WTYPE, nlev=NLEVEL):
    df = pd.DataFrame(index=trimmedindex(serie, nlev))
    # df['signal'] = serie.iloc[:len(df)]
    x = serie.iloc[:len(df)]
    for j in range(nlev):
        level = j+1
        ca, cd = pywt.dwt(x, wtype, mode='per')
        x = np.copy(ca)
        for i in range(level):
            apx = pywt.idwt(ca, None, wtype, mode= 'per')
            det = pywt.idwt(None, cd, wtype, mode= 'per')
            ca = apx
            cd = det
        for lbl, vec in zip(['approx', 'detail'], [apx, det]):
            label = 'level_{:d}_{:s}'.format(level, lbl)
            df[label] = vec
    colnames = []
    for j in range(nlev):
        level = j+1
        for lbl in ['approx', 'detail']:
            label = (level, lbl)
            colnames.append(label)
    df.columns = pd.MultiIndex.from_tuples(colnames, names=['level','type'])
    df[(0, 'signal')] = serie.iloc[:len(df)]
    return df.sort_index(axis=1)


def power_spectrum(data):
    x = data - data.mean()
    ham = signal.hamming(len(data))
    x = x*ham
    return np.abs(np.fft.fft(x))**2


# +
# m1_wvlt = {}
m3_wvlt = {}
m3_ffts = {}
levels = 10
fsize = (25,35)
for code in CODES:
    # m1_wvlt[code] = wvlt_analysis(m1_bin[code], nlev=levels)
    tmp = wvlt_analysis(m3_bin[code], nlev=levels)
    m3_wvlt[code] = tmp.copy(deep=True)
    m3_ffts[code] = tmp.apply(power_spectrum)
    freqs = np.fft.fftfreq(len(tmp), TIMESTEP)
    freqs[freqs <= 0] = np.nan
    m3_ffts[code]['freqs'] = freqs
    m3_ffts[code] = m3_ffts[code].dropna().set_index('freqs')
    
    
titles = [('Level {:d} :: approximation'.format(i), 
           'Level {:d} :: detail'.format(i)) for i in range(1, levels+1)]
titles = [item for sublist in titles for item in sublist]


# +
def plot_in_time(icao):
    f, axes = plt.subplots(levels, 2, figsize=fsize)
    tmp = m3_wvlt[icao].loc[:, (slice(1,levels), slice(None))]
    tmp.plot(ax=axes, subplots=True, colormap=CMAP,
             legend=False, title=titles)
    return (f, axes)


def plot_in_freq(icao):
    f, axes = plt.subplots(levels, 2, figsize=fsize)
    tmp = m3_ffts[icao].loc[:, (slice(1,levels), slice(None))]
    tmp.plot(ax=axes, subplots=True, colormap=CMAP,
             legend=False, title=titles)
    for ax in axes:
        ax[0].set_xticks(XTICKS)
        ax[0].set_xticklabels(XTKLBS)
        ax[0].set_xlim(right=CUTFRQ)
        ax[1].set_xticks(XTICKS)
        ax[1].set_xticklabels(XTKLBS)
        ax[1].set_xlim(left=CUTFRQ)
    return (f, axes)


# -

# # Denoising

code='EGLL'
deno = m3_bin[code].copy(deep=True)
noisy_coefs = pywt.wavedec(deno, 'db5', mode='per')

sigma = mad(noisy_coefs[-1])
uthresh = sigma*np.sqrt(2*np.log(len(deno)))
denoised = noisy_coefs[:]
denoised[1:] = (pywt.threshold(i, uthresh, 'soft') for i in denoised[1:])

m3_denoised = pywt.waverec(denoised, WTYPE, mode='per').flatten()
deno = deno.to_frame(name='original')
deno['denoised'] = pd.Series(m3_denoised.flatten(), index=m3.index)

f, ax = plt.subplots()
deno.loc['2016-08-01':'2016-08-07', 'original'].plot(ax=ax, color='cadetblue')
deno.loc['2016-08-01':'2016-08-07', 'denoised'].plot(ax=ax, color='navy')
plt.show()

# # Time-Frequency Analysis

code='EDDF'
m3 = m3_bin[code].copy(deep=True)
foo = m3.loc['2016-08-01':'2016-08-07']
xticks = [144*i for i in range(7)]
xticklabels = ['2016-08-01', '2016-08-02', '2016-08-03', '2016-08-04', '2016-08-05', '2016-08-06', '2016-08-07']

f, ax = plt.subplots()
widths = np.arange(1, 31)
cwtmatr, freqs = pywt.cwt(foo, widths, 'mexh')
im = ax.imshow(cwtmatr, extent=[0, len(foo), 1, 31], cmap='PRGn', aspect='auto',
               vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels, rotation=30)
plt.show()

# # Decomposition
#
# ## Frankfurt

code = 'EDDF'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## London Heathrow

code = 'EGLL'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## London Gatwick

code = 'EGKK'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## <NAME>

code = 'EHAM'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## <NAME>

code = 'LFPG'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## Madrid-Barajas

code = 'LEMD'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## <NAME>

code = 'LIRF'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# ## Athens International

code = 'LGAV'

f, axes = plot_in_time(code)
f.tight_layout()
plt.show()

f, axes = plot_in_freq(code)
f.tight_layout()
plt.show()

# # Disruptions

f, ax = plt.subplots(2, 2)
df = m3_wvlt['LIRF']
df.loc['2016-06-27':'2016-07-10',(7, 'approx')].plot(ax=ax[0,0],
                                                     title='Rome, Alitalia pilots strike, Jul 5')
df = m3_wvlt['EDDF']
df.loc['2016-07-14':'2016-07-27',(7, 'approx')].plot(ax=ax[0,1],
                                                     title='Frankfurt, Unknown event Jul 23')
df.loc['2016-08-23':'2016-09-05',(7, 'approx')].plot(ax=ax[1,0],
                                                     title='Frankfurt, woman evades security check Aug 31')
df = m3_wvlt['LFPG']
df.loc['2016-07-22':'2016-08-04',(7, 'approx')].plot(ax=ax[1,1],
                                                     title='Paris, Air France pilots strike, Jul 27 - Aug 2')
f.tight_layout()
plt.show()



StarcoderdataJupyterScriptsDedupFiltered