from __future__ import division
from os import path, remove
import numpy as np
import pandas as pd
import csv
from sklearn.model_selection import StratifiedShuffleSplit
from time import time
from matplotlib import pyplot as plt
import seaborn as sns
from tensorflow.contrib import rnn
from tensorflow.contrib import learn
import shutil
from tensorflow.contrib.learn.python.learn import learn_runner
from sklearn.metrics import r2_score
from fastdtw import fastdtw
from collections import OrderedDict
from scipy.spatial.distance import euclidean
from statsmodels.tsa.stattools import coint
from skopt.space.space import Integer, Real
from skopt import gp_minimize
from skopt.plots import plot_convergence
import pickle
import inspect
import dill
import sys
from pandas import read_csv
from pandas import datetime
from pandas.plotting import autocorrelation_plot
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.arima_model import ARIMA
from pandas import DataFrame
from sklearn.metrics import mean_squared_error
from price_history import PriceHistory
from arima.arima_estimator import ArimaEstimator
import warnings
from collections import OrderedDict

%matplotlib inline

"""
Explanation: http://machinelearningmastery.com/arima-for-time-series-forecasting-with-python/
http://ucanalytics.com/blogs/arima-models-manufacturing-case-study-example-part-3/
https://stats.stackexchange.com/questions/23036/estimating-same-model-over-multiple-time-series
The idea with ARIMA models is that the final residual should look like white noise otherwise, there is juice or information available in the data to extract
End of explanation
"""


data_path = '../../../../../Dropbox/data'

score_dic_filepath = data_path + "/arima/scoredic_non_stationary.npy"

ph_data_path = data_path + '/price_history'

csv_in = "../price_history_03_seq_start_suddens_trimmed.csv"

ph = PriceHistory(csv_in)

seqs = ph.extractAllSequences()

def custom_transformation(sequence, window):
    #logged_vals = np.log(sequence.values)
    logged_vals = sequence.values
    logged_seq = sequence.copy()
    logged_seq[logged_seq.index] = logged_vals
    return logged_seq - logged_seq.rolling(window=window).mean()

orig_seq = seqs[0]

window = 3

seq_tr = custom_transformation(sequence=seqs[0], window=window)

np.all(seq_tr == seq_tr)

seq = seq_tr[window-1:]

np.all(seq == seq)

plt.figure(figsize=(15,6))
seq.plot()
plt.show()

ae = ArimaEstimator(p_auto_regression_order=4, d_integration_level=1, q_moving_average=3, easy_mode=False)

target_len = 30

XX = seq[:-target_len].values.reshape((1, -1))
XX.shape

YY = seq[-target_len:].values.reshape((1, -1))
YY.shape

%%time
with warnings.catch_warnings():
    warnings.filterwarnings("ignore")
    ae.fit(inputs=XX, targets=YY)

preds = ae.preds.flatten()
yy = YY.flatten()
yy.shape, preds.shape

plt.figure(figsize=(15,6))
plt.plot(yy, label='real')
plt.plot(preds, label='preds')
plt.legend()
plt.show()

"""
Explanation: p: The number of lag observations included in the model, also called the lag order.
d: The number of times that the raw observations are differenced, also called the degree of differencing.
q: The size of the moving average window, also called the order of moving average.
End of explanation
"""


orig_seq.plot()

aa = np.zeros((orig_seq.size - window +1, orig_seq.size))
aa.shape

for ii, item in enumerate(aa):
    aa[ii, ii:ii+window] = 1/window

trans1 = aa.dot(orig_seq)
trans1.shape

from numpy.linalg import pinv

trans2 = orig_seq.rolling(window=window).mean()[window-1:].values
trans2.shape

np.allclose(trans1, trans2)

pA = pinv(aa)
pA.shape

reconstructed = pA.dot(trans1)
reconstructed.shape

plt.plot(reconstructed, label='reconstructed')
plt.plot(orig_seq.values, label='original')
plt.legend()
plt.show()

np.allclose(reconstructed, orig_seq.values, rtol=1e-3)

#np.sqrt(np.mean((reconstructed - orig_seq.values)**2)) #RMSE
np.corrcoef(orig_seq.values, reconstructed)

"""
Explanation: Reverse Moving average
https://stats.stackexchange.com/questions/67907/extract-data-points-from-moving-average
End of explanation
"""


seq.plot()

bb = 1 - aa
print bb.shape
bb

pOneMinusA = pinv(bb)
pOneMinusA.shape

custom_recons = pOneMinusA.dot(seq.values)
custom_recons.shape

np.allclose(custom_recons, orig_seq.values, rtol=1e-3)

plt.plot(custom_recons, label='reconstructed')
plt.plot(orig_seq.values, label='original')
plt.legend()
plt.show()

"""
Explanation: Conclusion
That is great we have managed to recover the moving average
Now let's see if we can recover the subtracted moving average
End of explanation
"""


np.corrcoef(orig_seq.values, custom_recons)

aa.shape, orig_seq.values.shape

trans4 = orig_seq.values[window-1:] - aa.dot(orig_seq.values)
np.allclose(trans4, seq.values)

bb = 1 - aa
bb.shape

trans3 = bb.dot(orig_seq.values)
trans3.shape

plt.plot(trans3)
#plt.plot(orig_seq.values)

"""
Explanation: WTF What is going on ??
End of explanation
"""


parameters = OrderedDict([
    ('p_auto_regression_order', range(6)), #0-5
    ('d_integration_level', range(3)), #0-2
    ('q_moving_average', range(6)), #0-5
])

def cartesian_coord(*arrays):
    grid = np.meshgrid(*arrays)        
    coord_list = [entry.ravel() for entry in grid]
    points = np.vstack(coord_list).T
    return points

cart = cartesian_coord(*parameters.values())
cart.shape

108 * 15 / 60 / 2 #in hours!

count_inds=100
random_state = np.random.RandomState(seed=16011984)

random_inds = random_state.choice(range(len(XX)), count_inds, replace=False)
xx = XX[random_inds]
yy = YY[random_inds]

filepath = 'scoredic.npy'

%%time
with warnings.catch_warnings():
    warnings.filterwarnings("ignore")
    for pp, dd, qq in cart:
        scoredic = np.load(filepath)[()] if path.isfile(filepath) else OrderedDict()
        cur_tuple = (pp, dd, qq)
        if cur_tuple in scoredic:
            continue
        else:
            ae = ArimaEstimator(p_auto_regression_order=pp, d_integration_level=dd, q_moving_average=qq,
                                easy_mode=True)
            scoredic[cur_tuple] = ae.fit(inputs=xx, targets=yy).score(xx, yy)
            np.save(filepath, scoredic)

from sklearn.model_selection import GridSearchCV

# Pick hyperparameters
# There is NO separation of training and testing set because there is no transfer of knowledge from the training 
# that could be useful for validation
# Run Arima prediction steps for all instances of training dataset and get a fastdtw score from each case
# Get the mean of the fastdtw and this is the CV score
# 
# Do this for all possible parameters

fitted = model.fit()

fitted.

model_fit = model.fit(disp=0)

model_fit.summary()

residuals = DataFrame(model_fit.resid)

residuals.plot()
plt.show()

residuals.plot(kind='kde')
plt.show()

residuals.describe()

"""
Explanation: 
End of explanation
"""


XX = series.values
XX.shape

#splitting
train_size = int(len(XX)*2/3)
train_size

train_set = XX[:train_size]
test_set = XX[train_size:]
train_set.shape, test_set.shape

train_set[0]

history = list(train_set)

history[0] = 3333

history = list(train_set)
len(history)

predictions = []

for tt in range(len(test_set)):
    output = ARIMA(history, order=(5,1,0)).fit(disp=0).forecast()
    y_hat = output[0] #0th is the index of the prediction
    predictions.append(y_hat)
    
    observation = test_set[tt]
    #history.append(observation)
    history.append(y_hat)
    
    print "predicted: {}, expected: {}".format(y_hat, observation)

error = mean_squared_error(predictions, test_set)
error

plt.figure(figsize=(15,6))
plt.plot(predictions, label='preds')
plt.plot(test_set, label='test set')
plt.legend()
plt.show()

"""
Explanation: Train - Test
End of explanation
"""


