# -*- coding: utf-8 -*-

"""

Copy of GAN with Generator: LSTM, Discriminator: Convolutional NN with ECG Data

 Introduction

 ------------

    The aim of this script is to use a convolutional neural network with 

    a max pooling layer in the discrimiantor. 

    This was found to work well with the Physionet ECG data in a paper. 

    They used two convolutional NN so we will compare the difference between the 

    images generated using a single layer of CNN in the discriminator and 2 CNN layers 

    to see if this improves the quality of series generated.

"""

"""

Bringing in required dependencies as defined in the GitHub repo: 

    https://github.com/josipd/torch-two-sample/blob/master/torch_two_sample/permutation_test.pyx"""

from __future__ import division

import torch

from tqdm import tqdm

import numpy as np

from matplotlib import pyplot as plt

import seaborn as sns

from torchvision import transforms

from torch.autograd.variable import Variable

sns.set(rc={'figure.figsize':(11, 4)})

import datetime 

from datetime import date

today = date.today()

import random

import json as js

import pickle

import os

from data import ECGData, PD_to_Tensor

from Model import Generator, Discriminator 

device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')

if device == 'cuda:0':

    print('Using GPU : ')

    print(torch.cuda.get_device_name(device))

else :

    print('Using CPU')

"""#MMD Evaluation Metric Definition

Using MMD to determine the similarity between distributions

PDIST code comes from torch-two-sample utils code: 

    https://github.com/josipd/torch-two-sample/blob/master/torch_two_sample/util.py

"""

def pdist(sample_1, sample_2, norm=2, eps=1e-5):

    r"""Compute the matrix of all squared pairwise distances.

    Arguments

    ---------

    sample_1 : torch.Tensor or Variable

        The first sample, should be of shape ``(n_1, d)``.

    sample_2 : torch.Tensor or Variable

        The second sample, should be of shape ``(n_2, d)``.

    norm : float

        The l_p norm to be used.

    Returns

    -------

    torch.Tensor or Variable

        Matrix of shape (n_1, n_2). The [i, j]-th entry is equal to

        ``|| sample_1[i, :] - sample_2[j, :] ||_p``."""

    n_1, n_2 = sample_1.size(0), sample_2.size(0)

    norm = float(norm)

    if norm == 2.:

        norms_1 = torch.sum(sample_1**2, dim=1, keepdim=True)

        norms_2 = torch.sum(sample_2**2, dim=1, keepdim=True)

        norms = (norms_1.expand(n_1, n_2) +

                 norms_2.transpose(0, 1).expand(n_1, n_2))

        distances_squared = norms - 2 * sample_1.mm(sample_2.t())

        return torch.sqrt(eps + torch.abs(distances_squared))

    else:

        dim = sample_1.size(1)

        expanded_1 = sample_1.unsqueeze(1).expand(n_1, n_2, dim)

        expanded_2 = sample_2.unsqueeze(0).expand(n_1, n_2, dim)

        differences = torch.abs(expanded_1 - expanded_2) ** norm

        inner = torch.sum(differences, dim=2, keepdim=False)

        return (eps + inner) ** (1. / norm)

def permutation_test_mat(matrix,

                         n_1,  n_2,  n_permutations,

                          a00=1,  a11=1,  a01=0):

    """Compute the p-value of the following statistic (rejects when high)

        \sum_{i,j} a_{\pi(i), \pi(j)} matrix[i, j].

    """

    n = n_1 + n_2

    pi = np.zeros(n, dtype=np.int8)

    pi[n_1:] = 1

    larger = 0.

    count = 0

    for sample_n in range(1 + n_permutations):

        count = 0.

        for i in range(n):

            for j in range(i, n):

                mij = matrix[i, j] + matrix[j, i]

                if pi[i] == pi[j] == 0:

                    count += a00 * mij

                elif pi[i] == pi[j] == 1:

                    count += a11 * mij

                else:

                    count += a01 * mij

        if sample_n == 0:

            statistic = count

        elif statistic <= count:

            larger += 1

        np.random.shuffle(pi)

    return larger / n_permutations

"""Code from Torch-Two-Samples at https://torch-two-sample.readthedocs.io/en/latest/#"""

class MMDStatistic:

    r"""The *unbiased* MMD test of :cite:`gretton2012kernel`.

    The kernel used is equal to:

    .. math ::

        k(x, x') = \sum_{j=1}^k e^{-\alpha_j\|x - x'\|^2},

    for the :math:`\alpha_j` proved in :py:meth:`~.MMDStatistic.__call__`.

    Arguments

    ---------

    n_1: int

        The number of points in the first sample.

    n_2: int

        The number of points in the second sample."""

    def __init__(self, n_1, n_2):

        self.n_1 = n_1

        self.n_2 = n_2

        # The three constants used in the test.

        self.a00 = 1. / (n_1 * (n_1 - 1))

        self.a11 = 1. / (n_2 * (n_2 - 1))

        self.a01 = - 1. / (n_1 * n_2)

    def __call__(self, sample_1, sample_2, alphas, ret_matrix=False):

        r"""Evaluate the statistic.

        The kernel used is

        .. math::

            k(x, x') = \sum_{j=1}^k e^{-\alpha_j \|x - x'\|^2},

        for the provided ``alphas``.

        Arguments

        ---------

        sample_1: :class:`torch:torch.autograd.Variable`

            The first sample, of size ``(n_1, d)``.

        sample_2: variable of shape (n_2, d)

            The second sample, of size ``(n_2, d)``.

        alphas : list of :class:`float`

            The kernel parameters.

        ret_matrix: bool

            If set, the call with also return a second variable.

            This variable can be then used to compute a p-value using

            :py:meth:`~.MMDStatistic.pval`.

        Returns

        -------

        :class:`float`

            The test statistic.

        :class:`torch:torch.autograd.Variable`

            Returned only if ``ret_matrix`` was set to true."""

        sample_12 = torch.cat((sample_1, sample_2), 0)

        distances = pdist(sample_12, sample_12, norm=2)

        kernels = None

        for alpha in alphas:

            kernels_a = torch.exp(- alpha * distances ** 2)

            if kernels is None:

                kernels = kernels_a

            else:

                kernels = kernels + kernels_a

        k_1 = kernels[:self.n_1, :self.n_1]

        k_2 = kernels[self.n_1:, self.n_1:]

        k_12 = kernels[:self.n_1, self.n_1:]

        mmd = (2 * self.a01 * k_12.sum() +

               self.a00 * (k_1.sum() - torch.trace(k_1)) +

               self.a11 * (k_2.sum() - torch.trace(k_2)))

        if ret_matrix:

            return mmd, kernels

        else:

            return mmd

    def pval(self, distances, n_permutations=1000):

        r"""Compute a p-value using a permutation test.

        Arguments

        ---------

        matrix: :class:`torch:torch.autograd.Variable`

            The matrix computed using :py:meth:`~.MMDStatistic.__call__`.

        n_permutations: int

            The number of random draws from the permutation null.

        Returns

        -------

        float

            The estimated p-value."""

        if isinstance(distances, Variable):

            distances = distances.data

        return permutation_test_mat(distances.cpu().numpy(),

                                    self.n_1, self.n_2,

                                    n_permutations,

                                    a00=self.a00, a11=self.a11, a01=self.a01)

"""

This paper 

https://arxiv.org/pdf/1611.04488.pdf says that the most common way to 

calculate sigma is to use the median pairwise distances between the joint data.

"""

def pairwisedistances(X,Y,norm=2):

    dist = pdist(X,Y,norm)

    return np.median(dist.numpy())

""" 

Function for loading ECG Data 

"""

def GetECGData(source_file,class_id):

  compose = transforms.Compose(

        [PD_to_Tensor()

        ])

  return ECGData(source_file ,class_id = class_id, transform = compose)

"""

Creating the training set of sine/ECG signals

"""

#Taking normal ECG data for now

source_filename = './mitbih_train.csv'

ecg_data = GetECGData(source_file = source_filename,class_id = 0)

sample_size = 119 #batch size needed for Data Loader and the noise creator function.

# Create loader with data, so that we can iterate over it

data_loader = torch.utils.data.DataLoader(ecg_data, batch_size=sample_size, shuffle=True)

# Num batches

num_batches = len(data_loader)

print(num_batches)

"""Creating the Test Set"""

test_filename =  './mitbih_test.csv'

ecg_data_test = GetECGData(source_file = test_filename,class_id = 0)

data_loader_test = torch.utils.data.DataLoader(ecg_data_test[:18088], batch_size=sample_size, shuffle=True)

"""##Defining the noise creation function"""

def noise(batch_size, features):

  noise_vec = torch.randn(batch_size, features).to(device)

  return noise_vec

"""#Initialising Parameters"""

seq_length = ecg_data[0].size()[0] #Number of features

#Params for the generator

hidden_nodes_g = 50

layers = 2

tanh_layer = False

#No. of training rounds per epoch

D_rounds = 3

G_rounds = 1

num_epoch = 35

learning_rate = 0.0002

#Params for the Discriminator

minibatch_layer = 0

minibatch_normal_init_ = True

num_cvs = 2

cv1_out= 10

cv1_k = 3

cv1_s = 1

p1_k = 3

p1_s = 2

cv2_out = 10

cv2_k = 3

cv2_s = 1

p2_k = 3

p2_s = 2

"""# Evaluation of GAN with 2 CNN Layer in Discriminator

##Generator and Discriminator training phase

"""

minibatch_out = [0,3,5,8,10]

for minibatch_layer in minibatch_out:

  path = ".../your_path/Run_"+str(today.strftime("%d_%m_%Y"))+"_"+ str(datetime.datetime.now().time()).split('.')[0]

  os.mkdir(path)

  dict = {'data' : source_filename, 

          'sample_size' : sample_size, 

          'seq_length' : seq_length,

          'num_layers': layers, 

          'tanh_layer': tanh_layer,

          'hidden_dims_generator': hidden_nodes_g, 

          'minibatch_layer': minibatch_layer,

          'minibatch_normal_init_' : minibatch_normal_init_,

          'num_cvs':num_cvs,

          'cv1_out':cv1_out,

          'cv1_k':cv1_k,

          'cv1_s':cv1_s,

          'p1_k':p1_k,

          'p1_s':p1_s,

          'cv2_out':cv2_out,

          'cv2_k':cv2_k,

          'cv2_s':cv2_s,

          'p2_k':p2_k,

          'p2_s':p2_s,

          'num_epoch':num_epoch,

          'D_rounds': D_rounds,

          'G_rounds': G_rounds,  

          'learning_rate' : learning_rate

         }

  json = js.dumps(dict)

  f = open(path+"/settings.json","w")

  f.write(json)

  f.close()

  generator_1 = Generator(seq_length,sample_size,hidden_dim =  hidden_nodes_g, tanh_output = tanh_layer).to(device)

  discriminator_1 = Discriminator(seq_length, sample_size ,minibatch_normal_init = minibatch_normal_init_, minibatch = minibatch_layer,num_cv = num_cvs, cv1_out = cv1_out,cv1_k = cv1_k, cv1_s = cv1_s, p1_k = p1_k, p1_s = p1_s, cv2_out= cv2_out, cv2_k = cv2_k, cv2_s = cv2_s, p2_k = p2_k, p2_s = p2_s).to(device)

  #Loss function 

  loss_1 = torch.nn.BCELoss()

  generator_1.train()

  discriminator_1.train()

  d_optimizer_1 = torch.optim.Adam(discriminator_1.parameters(),lr = learning_rate)

  g_optimizer_1 = torch.optim.Adam(generator_1.parameters(),lr = learning_rate)

  G_losses = []

  D_losses = []

  mmd_list = []

  series_list = np.zeros((1,seq_length))

  for n in tqdm(range(num_epoch)):

     # for k in range(1):

      for n_batch, sample_data in enumerate(data_loader):

      ### TRAIN DISCRIMINATOR ON FAKE DATA

        for d in range(D_rounds):

          discriminator_1.zero_grad()

          h_g = generator_1.init_hidden()

          #Generating the noise and label data

          noise_sample = Variable(noise(len(sample_data),seq_length))

          #Use this line if generator outputs hidden states: dis_fake_data, (h_g_n,c_g_n) = generator.forward(noise_sample,h_g)

          dis_fake_data = generator_1.forward(noise_sample,h_g).detach()

          y_pred_fake = discriminator_1(dis_fake_data)

          loss_fake = loss_1(y_pred_fake,torch.zeros([len(sample_data),1]).to(device))

          loss_fake.backward()    

          #Train discriminator on real data   

          real_data = Variable(sample_data.float()).to(device)    

          y_pred_real  = discriminator_1.forward(real_data)

          loss_real = loss_1(y_pred_real,torch.ones([len(sample_data),1]).to(device))

          loss_real.backward()

          d_optimizer_1.step() #Updating the weights based on the predictions for both real and fake calculations.

        #Train Generator  

        for g in range(G_rounds):

          generator_1.zero_grad()

          h_g = generator_1.init_hidden()

          noise_sample = Variable(noise(len(sample_data), seq_length))

          #Use this line if generator outputs hidden states: gen_fake_data, (h_g_n,c_g_n) = generator.forward(noise_sample,h_g)

          gen_fake_data = generator_1.forward(noise_sample,h_g)

          y_pred_gen = discriminator_1(gen_fake_data)

          error_gen = loss_1(y_pred_gen,torch.ones([len(sample_data),1]).to(device))

          error_gen.backward()

          g_optimizer_1.step()         

      if n_batch ==( num_batches - 1):

          G_losses.append(error_gen.item())

          D_losses.append((loss_real+loss_fake).item())

          torch.save(generator_1.state_dict(), path+'/generator_state_'+str(n)+'.pt')

          torch.save(discriminator_1.state_dict(),path+ '/discriminator_state_'+str(n)+'.pt')

        # Check how the generator is doing by saving G's output on fixed_noise

          with torch.no_grad():

              h_g = generator_1.init_hidden()

              fake = generator_1(noise(len(sample_data), seq_length),h_g).detach().cpu()

              generated_sample = torch.zeros(1,seq_length).to(device)

              for iter in range(0,int(len(ecg_data_test[:18088])/sample_size)):

                noise_sample_test = noise(sample_size, seq_length)

                h_g = generator_1.init_hidden()

                generated_data = generator_1.forward(noise_sample_test,h_g).detach().squeeze()

                generated_sample = torch.cat((generated_sample,generated_data),dim = 0)

              # Getting the MMD Statistic for each Training Epoch

              generated_sample = generated_sample[1:][:]

              sigma = [pairwisedistances(ecg_data_test[:18088].type(torch.DoubleTensor),generated_sample.type(torch.DoubleTensor).squeeze())] 

              mmd = MMDStatistic(len(ecg_data_test[:18088]),generated_sample.size(0))

              mmd_eval = mmd(ecg_data_test[:18088].type(torch.DoubleTensor),generated_sample.type(torch.DoubleTensor).squeeze(),sigma, ret_matrix=False)

              mmd_list.append(mmd_eval.item())

              series_list = np.append(series_list,fake[0].numpy().reshape((1,seq_length)),axis=0)

  #Dumping the errors and mmd evaluations for each training epoch.

  with open(path+'/generator_losses.txt', 'wb') as fp:

      pickle.dump(G_losses, fp)

  with open(path+'/discriminator_losses.txt', 'wb') as fp:

      pickle.dump(D_losses, fp)   

  with open(path+'/mmd_list.txt', 'wb') as fp:

      pickle.dump(mmd_list, fp)

  #Plotting the error graph

  plt.plot(G_losses,'-r',label='Generator Error')

  plt.plot(D_losses, '-b', label = 'Discriminator Error')

  plt.title('GAN Errors in Training')

  plt.legend()

  plt.savefig(path+'/GAN_errors.png')

  plt.close() 

  #Plot a figure for each training epoch with the MMD value in the title

  i = 0

  while i < num_epoch:

    if i%3==0:

      fig, ax = plt.subplots(3,1,constrained_layout=True)

      fig.suptitle("Generated fake data")

    for j in range(0,3):

      ax[j].plot(series_list[i][:])

      ax[j].set_title('Epoch '+str(i)+ ', MMD: %.4f' % (mmd_list[i]))

      i = i+1

    plt.savefig(path+'/Training_Epoch_Samples_MMD_'+str(i)+'.png')

    plt.close(fig) 

  #Checking the diversity of the samples:

  generator_1.eval()

  h_g = generator_1.init_hidden()

  test_noise_sample = noise(sample_size, seq_length)

  gen_data= generator_1.forward(test_noise_sample,h_g).detach()

  plt.title("Generated ECG Waves")

  plt.plot(gen_data[random.randint(0,sample_size-1)].tolist(),'-b')

  plt.plot(gen_data[random.randint(0,sample_size-1)].tolist(),'-r')

  plt.plot(gen_data[random.randint(0,sample_size-1)].tolist(),'-g')

  plt.plot(gen_data[random.randint(0,sample_size-1)].tolist(),'-', color = 'orange')

  plt.savefig(path+'/Generated_Data_Sample1.png')

  plt.close()