<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Deep LSTM RNNs for N-D time series
Step2:  Dataset
Step3:  Check the data real quick
Step4:  Preparing the data for training
Step5:  Long short-term memory (LSTM) RNNs
An LSTM block has mechanisms to enable "memorizing" information for an extended number of time steps. We use the LSTM block with the following transformations that map inputs to outputs across blocks at consecutive layers and consecutive time steps
Step6:  Attach the gradients
Step7:  Softmax Activation
Step8:  Cross-entropy loss function
Step9:  Averaging the loss over the sequence
Step10:  Optimizer
Step12:  Define the model
Step13:  Test and visualize predictions
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function
import mxnet as mx
from mxnet import nd, autograd
import numpy as np
from collections import defaultdict
mx.random.seed(1)
# ctx = mx.gpu(0)
ctx = mx.cpu(0)
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import scipy.fftpack
from pandas.tools import plotting
from pandas.tools.plotting import autocorrelation_plot
from datetime import datetime
sns.set_style('whitegrid')
#sns.set_context('notebook')
sns.set_context('poster')
# Make inline plots vector graphics instead of raster graphics
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'png')
Explanation: Deep LSTM RNNs for N-D time series
End of explanation
def get_list_unique_block_indices(len_data=100, seq_length=5, n_samples=10):
     returns a list of unique random int that serve as index of the first element of a block of data
        args:
            len_data (int): length of the data set
            seq_length (int): length of the blocks to extract
            n_blocks (int): # of blocks to extract
    
    set1 = set(np.random.randint(len_data // seq_length, size=n_samples)*seq_length)
    full_set = set1
    while len(full_set) < n_samples:
        set2 = set(np.random.randint(len_data // seq_length, size=n_samples)*seq_length)
        full_set = full_set | set2
    returned_list = list(full_set)[0:n_samples]
    assert(len(returned_list) == n_samples)
    return returned_list
def extract_random_sequence(data, seq_length=5, block_start_index=None):
    columns_subset = ['car.count', 'day_of_week_int', 'cloudy_or_not_cloudy', 'weather', 'current_month']
    if block_start_index is None:
        block_start_index = np.random.randint(len(data)-seq_length)
    data_subset = data.reset_index().loc[block_start_index:block_start_index+seq_length-1, columns_subset]
    assert(len(data_subset) == (seq_length))
    out_data = [list(i) for i in data_subset.values]
    return out_data
def create_batch_ND_time_series(full_data, seq_length=10, num_samples=4):
    out_data = []
    # get a list of non-overlapping random sequence start indices
    all_samples_start_indices = get_list_unique_block_indices(len(full_data), seq_length, num_samples)
    assert(len(all_samples_start_indices) == num_samples)
    for one_random_start_index in all_samples_start_indices:
        out_data.append(extract_random_sequence(full_data, seq_length, one_random_start_index))
        assert(len(out_data[-1]) == (seq_length))
    return out_data
# OI data
original_data = pd.read_csv("../data/timeseries/data.csv", index_col=0)
dict_days_to_int = {'Monday': 1, 'Tuesday': 2, 'Wednesday': 3, 'Thursday': 4, 'Friday': 5, 'Saturday': 6, 'Sunday': 7}
original_data['date_']=original_data.index
original_data['current_month'] = original_data['date_'].apply(lambda x: pd.Timestamp(x).month)
original_data['day_of_week_int'] = original_data['day.of.week'].apply(lambda x: dict_days_to_int[x])
original_data['cloudy_or_not_cloudy'] = original_data['cloud.indicator'].apply(lambda x: 0 if x=='clear' else 1)
full_data = pd.DataFrame()
#############################
# JUST A RDM SAMPLE FOR NOW   /!\ /!\ /!\ /!\ /!\ /!\ /!\ careful as it shuffles the time order!
#############################
full_data = original_data#.sample(1500)
# filter out cloudy data!
full_data = full_data[full_data['cloudy_or_not_cloudy']==0]
SEQ_LENGTH = 2
NUM_FEATURES = 5
# let's divide data in train (75%), dev (15%), test (10%)
# in sequences of 5 days (SEQ_LENGTH = 5)
full_data_length = len(full_data)
# the actual length of extracted sequence is SEQ_LENGTH + 1 so that we can do the shift of +1 for labels
total_num_of_sequences = full_data_length // (SEQ_LENGTH+1) - 1
# the length of extracted sequence is SEQ_LENGTH so that we can do the shift of +1 for labels
all_random_sequences = create_batch_ND_time_series(full_data, seq_length=SEQ_LENGTH+1, num_samples=total_num_of_sequences)
n_seq_train = int(total_num_of_sequences*0.75)
n_seq_dev = int(total_num_of_sequences*0.9) - int(total_num_of_sequences*0.75)
n_seq_test = len(all_random_sequences) - int(total_num_of_sequences*0.9)
data_train = np.array(all_random_sequences[0:n_seq_train])
data_dev = np.array(all_random_sequences[n_seq_train:n_seq_train+n_seq_dev])
data_test = np.array(all_random_sequences[n_seq_train+n_seq_dev:])
print('SHAPES of ALL, TRAIN, DEV, TEST:') 
print(np.array(all_random_sequences).shape)
print(np.array(data_train).shape)
print(np.array(data_dev).shape)
print(np.array(data_test).shape)
assert(data_train.shape == (n_seq_train, SEQ_LENGTH+1, NUM_FEATURES))
assert(data_dev.shape == (n_seq_dev, SEQ_LENGTH+1, NUM_FEATURES))
assert(data_test.shape == (n_seq_test, SEQ_LENGTH+1, NUM_FEATURES))
Explanation: Dataset: "Some time-series"
End of explanation
# num_sampling_points = min(SEQ_LENGTH, 400)
# (data_train.sample(4).transpose().iloc[range(0, SEQ_LENGTH, SEQ_LENGTH//num_sampling_points)]).plot()
# print (data_train)
# print(data_train[:, :-1, :]) # inputs
# batch_size = 5
# num_batches_train = data_train.shape[0] // batch_size
# print(num_batches_train)
# print ( nd.array(data_train[:, :-1, :]).reshape((num_batches_train, 5, SEQ_LENGTH-1, NUM_FEATURES)) )
Explanation: Check the data real quick
End of explanation
batch_size = 32
batch_size_test = 1
seq_length = SEQ_LENGTH
num_batches_train = data_train.shape[0] // batch_size
num_batches_test = data_test.shape[0] // batch_size_test
num_features = NUM_FEATURES  #  we do 1D time series for now, this is like vocab_size = 1 for characters
# inputs are from t0 to t_seq_length - 1. because the last point is kept for the 
# output ("label") of the penultimate point 
data_train_inputs = data_train[:, :-1, :]
data_train_labels = data_train[:, 1:, :]
data_test_inputs = data_test[:, :-1, :]
data_test_labels = data_test[:, 1:, :]
train_data_inputs = nd.array(data_train_inputs).reshape((num_batches_train, batch_size, seq_length, num_features))
train_data_labels = nd.array(data_train_labels).reshape((num_batches_train, batch_size, seq_length, num_features))
test_data_inputs = nd.array(data_test_inputs).reshape((num_batches_test, batch_size_test, seq_length, num_features))
test_data_labels = nd.array(data_test_labels).reshape((num_batches_test, batch_size_test, seq_length, num_features))
train_data_inputs = nd.swapaxes(train_data_inputs, 1, 2)
train_data_labels = nd.swapaxes(train_data_labels, 1, 2)
test_data_inputs = nd.swapaxes(test_data_inputs, 1, 2)
test_data_labels = nd.swapaxes(test_data_labels, 1, 2)
print('num_mini-batches_train={0} | seq_length={2} | mini-batch_size={1} | num_features={3}'.format(num_batches_train, batch_size, seq_length, num_features))
print('train_data_inputs shape: ', train_data_inputs.shape)
print('train_data_labels shape: ', train_data_labels.shape)
# print(data_train_inputs.values)
# print(train_data_inputs[0]) # see what one batch looks like
Explanation: Preparing the data for training
End of explanation
num_inputs = NUM_FEATURES  # for a 1D time series, this is just a scalar equal to 1.0
num_outputs = NUM_FEATURES  # same comment
num_hidden_units = [64]  # num of hidden units in each hidden LSTM layer
num_hidden_layers = len(num_hidden_units)  # num of hidden LSTM layers
num_units_layers = [num_features] + num_hidden_units
########################
#  Weights connecting the inputs to the hidden layer
########################
Wxg, Wxi, Wxf, Wxo, Whg, Whi, Whf, Who, bg, bi, bf, bo = {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {} 
for i_layer in range(1, num_hidden_layers+1):
    num_inputs = num_units_layers[i_layer-1]
    num_hidden_units = num_units_layers[i_layer]
    Wxg[i_layer] = nd.random_normal(shape=(num_inputs,num_hidden_units), ctx=ctx) * .01
    Wxi[i_layer] = nd.random_normal(shape=(num_inputs,num_hidden_units), ctx=ctx) * .01
    Wxf[i_layer] = nd.random_normal(shape=(num_inputs,num_hidden_units), ctx=ctx) * .01
    Wxo[i_layer] = nd.random_normal(shape=(num_inputs,num_hidden_units), ctx=ctx) * .01
    ########################
    #  Recurrent weights connecting the hidden layer across time steps
    ########################
    Whg[i_layer] = nd.random_normal(shape=(num_hidden_units, num_hidden_units), ctx=ctx) * .01
    Whi[i_layer] = nd.random_normal(shape=(num_hidden_units, num_hidden_units), ctx=ctx) * .01
    Whf[i_layer] = nd.random_normal(shape=(num_hidden_units, num_hidden_units), ctx=ctx) * .01
    Who[i_layer] = nd.random_normal(shape=(num_hidden_units, num_hidden_units), ctx=ctx) * .01
    ########################
    #  Bias vector for hidden layer
    ########################
    bg[i_layer] = nd.random_normal(shape=num_hidden_units, ctx=ctx) * .01
    bi[i_layer] = nd.random_normal(shape=num_hidden_units, ctx=ctx) * .01
    bf[i_layer] = nd.random_normal(shape=num_hidden_units, ctx=ctx) * .01
    bo[i_layer] = nd.random_normal(shape=num_hidden_units, ctx=ctx) * .01
########################
# Weights to the output nodes
########################
Why = nd.random_normal(shape=(num_units_layers[-1], num_outputs), ctx=ctx) * .01
by = nd.random_normal(shape=num_outputs, ctx=ctx) * .01
Explanation: Long short-term memory (LSTM) RNNs
An LSTM block has mechanisms to enable "memorizing" information for an extended number of time steps. We use the LSTM block with the following transformations that map inputs to outputs across blocks at consecutive layers and consecutive time steps: $\newcommand{\xb}{\mathbf{x}} \newcommand{\RR}{\mathbb{R}}$
$$g_t = \text{tanh}(X_t W_{xg} + h_{t-1} W_{hg} + b_g),$$
$$i_t = \sigma(X_t W_{xi} + h_{t-1} W_{hi} + b_i),$$
$$f_t = \sigma(X_t W_{xf} + h_{t-1} W_{hf} + b_f),$$
$$o_t = \sigma(X_t W_{xo} + h_{t-1} W_{ho} + b_o),$$
$$c_t = f_t \odot c_{t-1} + i_t \odot g_t,$$
$$h_t = o_t \odot \text{tanh}(c_t),$$
where $\odot$ is an element-wise multiplication operator, and
for all $\xb = [x_1, x_2, \ldots, x_k]^\top \in \RR^k$ the two activation functions:
$$\sigma(\xb) = \left[\frac{1}{1+\exp(-x_1)}, \ldots, \frac{1}{1+\exp(-x_k)}]\right]^\top,$$
$$\text{tanh}(\xb) = \left[\frac{1-\exp(-2x_1)}{1+\exp(-2x_1)},  \ldots, \frac{1-\exp(-2x_k)}{1+\exp(-2x_k)}\right]^\top.$$
In the transformations above, the memory cell $c_t$ stores the "long-term" memory in the vector form.
In other words, the information accumulatively captured and  encoded  until time step $t$ is stored in $c_t$ and is only passed along the same layer over different time steps.
Given the inputs $c_t$ and $h_t$, the input gate $i_t$ and forget gate $f_t$ will help the memory cell to decide how to overwrite or keep the memory information. The output gate $o_t$ further lets the LSTM block decide how to retrieve the memory information to generate the current state $h_t$ that is passed to both the next layer of the current time step and the next time step of the current layer. Such decisions are made using the hidden-layer parameters $W$ and $b$ with different subscripts: these parameters will be inferred during the training phase by gluon.
Allocate parameters
End of explanation
params = []
for i_layer in range(1, num_hidden_layers+1):
    params += [Wxg[i_layer], Wxi[i_layer], Wxf[i_layer], Wxo[i_layer], Whg[i_layer], Whi[i_layer], Whf[i_layer], Who[i_layer], bg[i_layer], bi[i_layer], bf[i_layer], bo[i_layer]]
params += [Why, by]  # add the output layer
for param in params:
    param.attach_grad()
Explanation: Attach the gradients
End of explanation
def softmax(y_linear, temperature=1.0):
    lin = (y_linear-nd.max(y_linear)) / temperature
    exp = nd.exp(lin)
    partition = nd.sum(exp, axis=0, exclude=True).reshape((-1,1))
    return exp / partition
Explanation: Softmax Activation
End of explanation
def cross_entropy(yhat, y):
    return - nd.mean(nd.sum(y * nd.log(yhat), axis=0, exclude=True))
def rmse(yhat, y):
    return nd.mean(nd.sqrt(nd.sum(nd.power(y - yhat, 2), axis=0, exclude=True)))
Explanation: Cross-entropy loss function
End of explanation
def average_ce_loss(outputs, labels):
    assert(len(outputs) == len(labels))
    total_loss = 0.
    for (output, label) in zip(outputs,labels):
        total_loss = total_loss + cross_entropy(output, label)
    return total_loss / len(outputs)
def average_rmse_loss(outputs, labels):
    assert(len(outputs) == len(labels))
    total_loss = 0.
    for (output, label) in zip(outputs,labels):
        total_loss = total_loss + rmse(output, label)
    return total_loss / len(outputs)
Explanation: Averaging the loss over the sequence
End of explanation
from exceptions import ValueError
def SGD(params, learning_rate):
    for param in params:
#         print('grrrrr: ', param.grad)
        param[:] = param - learning_rate * param.grad
def adam(params, learning_rate, M , R, index_adam_call, beta1, beta2, eps):
    k = -1
    for param in params:
        k += 1
        M[k] = beta1 * M[k] + (1. - beta1) * param.grad
        R[k] = beta2 * R[k] + (1. - beta2) * (param.grad)**2
        # bias correction since we initilized M & R to zeros, they're biased toward zero on the first few iterations
        m_k_hat = M[k] / (1. - beta1**(index_adam_call))
        r_k_hat = R[k] / (1. - beta2**(index_adam_call))
        if((np.isnan(M[k].asnumpy())).any() or (np.isnan(R[k].asnumpy())).any()):
#             print('GRRRRRR ', M, K)
            raise(ValueError('Nans!!'))
#         print('grrrrr: ', param.grad)
        param[:] = param - learning_rate * m_k_hat / (nd.sqrt(r_k_hat) + eps)
#     print('m_k_hat r_k_hat', m_k_hat, r_k_hat)
    return params, M, R
Explanation: Optimizer
End of explanation
def single_lstm_unit_calcs(X, c, Wxg, h, Whg, bg, Wxi, Whi, bi, Wxf, Whf, bf, Wxo, Who, bo):
    g = nd.tanh(nd.dot(X, Wxg) + nd.dot(h, Whg) + bg)
    i = nd.sigmoid(nd.dot(X, Wxi) + nd.dot(h, Whi) + bi)
    f = nd.sigmoid(nd.dot(X, Wxf) + nd.dot(h, Whf) + bf)
    o = nd.sigmoid(nd.dot(X, Wxo) + nd.dot(h, Who) + bo)
    #######################
    c = f * c + i * g
    h = o * nd.tanh(c)
    return c, h
def deep_lstm_rnn(inputs, h, c, temperature=1.0):
    
        h: dict of nd.arrays, each key is the index of a hidden layer (from 1 to whatever). 
        Index 0, if any, is the input layer
    
    outputs = []
    # inputs is one BATCH of sequences so its shape is number_of_seq, seq_length, features_dim 
    # (latter is 1 for a time series, vocab_size for a character, n for a n different times series)
    for X in inputs:
        # X is batch of one time stamp. E.g. if each batch has 37 sequences, then the first value of X will be a set of the 37 first values of each of the 37 sequences 
        # that means each iteration on X corresponds to one time stamp, but it is done in batches of different sequences
        h[0] = X # the first hidden layer takes the input X as input 
        for i_layer in range(1, num_hidden_layers+1):
            # lstm units now have the 2 following inputs: 
            # i) h_t from the previous layer (equivalent to the input X for a non-deep lstm net), 
            # ii) h_t-1 from the current layer (same as for non-deep lstm nets)
            c[i_layer], h[i_layer] = single_lstm_unit_calcs(h[i_layer-1], c[i_layer], Wxg[i_layer], h[i_layer], Whg[i_layer], bg[i_layer], Wxi[i_layer], Whi[i_layer], bi[i_layer], Wxf[i_layer], Whf[i_layer], bf[i_layer], Wxo[i_layer], Who[i_layer], bo[i_layer])
        yhat_linear = nd.dot(h[num_hidden_layers], Why) + by
        # yhat is a batch of several values of the same time stamp
        # this is basically the prediction of the sequence, which overlaps most of the input sequence, plus one point (character or value)
#         yhat = softmax(yhat_linear, temperature=temperature)
#         yhat = nd.sigmoid(yhat_linear)
#         yhat = nd.tanh(yhat_linear)
        yhat = yhat_linear # we cant use a 1.0-bounded activation function since amplitudes can be greater than 1.0
        outputs.append(yhat) # outputs has same shape as inputs, i.e. a list of batches of data points.
#     print('some shapes... yhat outputs', yhat.shape, len(outputs) )
    return (outputs, h, c)
Explanation: Define the model
End of explanation
INDEX_TARGET_VALUE = 0
def test_prediction(one_input_seq, one_label_seq, temperature=1.0):
      # WE ASSUME the first value in input vector is the variable of interest
    #####################################
    # Set the initial state of the hidden representation ($h_0$) to the zero vector
    #####################################  # some better initialization needed??
    h, c = {}, {}
    for i_layer in range(1, num_hidden_layers+1):
        h[i_layer] = nd.zeros(shape=(batch_size_test, num_units_layers[i_layer]), ctx=ctx)
        c[i_layer] = nd.zeros(shape=(batch_size_test, num_units_layers[i_layer]), ctx=ctx)
    
    outputs, h, c = deep_lstm_rnn(one_input_seq, h, c, temperature=temperature)
    return outputs[-1][0].asnumpy()[INDEX_TARGET_VALUE], one_label_seq.asnumpy()[-1].flatten()[INDEX_TARGET_VALUE], outputs, one_label_seq
def check_prediction(index):
    if index >= len(test_data_inputs):
        index = np.random.randint(len(test_data_inputs))
        
    o, label, outputs, labels = test_prediction(test_data_inputs[index], test_data_labels[index], temperature=1.0)    
    prediction = round(o, 3)
    true_label = round(label, 3)
    outputs = [float(i.asnumpy().flatten()[INDEX_TARGET_VALUE]) for i in outputs]  # if batch_size_test=1 then this float() will work, otherwise, nope.
    true_labels = list(test_data_labels[index].asnumpy()[:,:,INDEX_TARGET_VALUE].flatten())
    
    df = pd.DataFrame([outputs, true_labels]).transpose()
    df.columns = ['predicted', 'true']
    if true_label != 0:
        rel_error = round(100. * (prediction / (true_label+1e-5) - 1.0), 2)
    else:
        rel_error = 100.
#     print('\nprediction = {0} | actual_value = {1} | rel_error = {2}'.format(prediction, true_label, rel_error))
    return df
epochs = 10000  # one epoch is one pass over the entire training set
moving_loss = 0.
learning_rate = 0.001  # 0.1 works for a [8, 8] after about 70 epochs of 32-sized batches
# Adam Optimizer stuff
beta1 = .9
beta2 = .999
index_adam_call = 0
# M & R arrays to keep track of momenta in adam optimizer. params is a list that contains all ndarrays of parameters
M = {k: nd.zeros_like(v) for k, v in enumerate(params)}
R = {k: nd.zeros_like(v) for k, v in enumerate(params)}
df_moving_loss = pd.DataFrame(columns=['Loss', 'Error'])
df_moving_loss.index.name = 'Epoch'
# needed to update plots on the fly
%matplotlib notebook
fig, axes_fig1 = plt.subplots(1,1, figsize=(6,3))
fig2, axes_fig2 = plt.subplots(1,1, figsize=(6,3))
for e in range(epochs):
    ############################
    # Attenuate the learning rate by a factor of 2 every 100 epochs
    ############################
    if ((e+1) % 1000 == 0):
        learning_rate = learning_rate / 2.0  # TODO check if its ok to adjust learning_rate when using Adam Optimizer
    h, c = {}, {}
    for i_layer in range(1, num_hidden_layers+1):
        h[i_layer] = nd.zeros(shape=(batch_size, num_units_layers[i_layer]), ctx=ctx)
        c[i_layer] = nd.zeros(shape=(batch_size, num_units_layers[i_layer]), ctx=ctx)
    for i in range(num_batches_train):
        data_one_hot = train_data_inputs[i]
        label_one_hot = train_data_labels[i]
        with autograd.record():
            outputs, h, c = deep_lstm_rnn(data_one_hot, h, c)
            loss = average_rmse_loss(outputs, label_one_hot)
            loss.backward()
#         SGD(params, learning_rate)
        index_adam_call += 1  # needed for bias correction in Adam optimizer
        params, M, R = adam(params, learning_rate, M, R, index_adam_call, beta1, beta2, 1e-8)
        
        ##########################
        #  Keep a moving average of the losses
        ##########################
        if (i == 0) and (e == 0):
            moving_loss = nd.mean(loss).asscalar()
        else:
            moving_loss = .99 * moving_loss + .01 * nd.mean(loss).asscalar()
        df_moving_loss.loc[e] = round(moving_loss, 4)
    ############################
    #  Predictions and plots
    ############################
    data_prediction_df = check_prediction(index=e)
    
    if not (e%50):
        axes_fig1.clear()
        data_prediction_df.plot(ax=axes_fig1)
        fig.canvas.draw()
    prediction = round(data_prediction_df.tail(1)['predicted'].values.flatten()[-1], 3)
    true_label = round(data_prediction_df.tail(1)['true'].values.flatten()[-1], 3)
    if true_label != 0:
        rel_error = round(100. * np.abs(prediction / (true_label+1e-5) - 1.0), 2)
    else:
        rel_error = moving_rel_error
    if not (e%50):
        print("Epoch = {0} | Loss = {1} | Prediction = {2} True = {3} Error = {4}".format(e, moving_loss, prediction, true_label, rel_error ))
    if not (e%50):
        axes_fig2.clear()
    if e == 0:
        moving_rel_error = rel_error
    else:
        moving_rel_error = .99 * moving_rel_error + .01 * rel_error
    df_moving_loss.loc[e, ['Error']] = moving_rel_error
    if not (e%50):
        axes_loss_plot = df_moving_loss.plot(ax=axes_fig2, secondary_y='Loss', color=['r','b'])
        axes_loss_plot.right_ax.grid(False)
#     axes_loss_plot.right_ax.set_yscale('log')
        fig2.canvas.draw()
    
%matplotlib inline
# print(outputs[0].asnumpy()[0].flatten())
# print(test_data_labels[0].asnumpy()[:,:,0].flatten())
# [float(i.asnumpy().flatten()) for i in outputs]
# print([i.asnumpy() for i in outputs])
# one_label_seq = test_data_labels[0]
# print(outputs[-1][0])
# print(one_label_seq)
# print(rmse(outputs[-1][0], one_label_seq))
# print(test_data_inputs[0].asnumpy()[-1].flatten()[0])
# print([i.asnumpy().flatten() for i in outputs])
# print(float(outputs[0].asnumpy()[:, 0].flatten()))
Explanation: Test and visualize predictions
End of explanation
<END_TASK>