#!/usr/bin/env python

"""

features_ECG.py

VARPA, University of Coruna

Mondejar Guerra, Victor M.

23 Oct 2017

"""

import numpy as np

from scipy.signal import medfilt

import scipy.stats

import pywt

import operator

from mit_db import *

# Input: the R-peaks from a signal

# Return: the features RR intervals 

#   (pre_RR, post_RR, local_RR, global_RR)

#    for each beat 

def compute_RR_intervals(R_poses):

    features_RR = RR_intervals()

    pre_R = np.array([], dtype=int)

    post_R = np.array([], dtype=int)

    local_R = np.array([], dtype=int)

    global_R = np.array([], dtype=int)

    # Pre_R and Post_R

    pre_R = np.append(pre_R, 0)

    post_R = np.append(post_R, R_poses[1] - R_poses[0])

    for i in range(1, len(R_poses)-1):

        pre_R = np.append(pre_R, R_poses[i] - R_poses[i-1])

        post_R = np.append(post_R, R_poses[i+1] - R_poses[i])

    pre_R[0] = pre_R[1]

    pre_R = np.append(pre_R, R_poses[-1] - R_poses[-2])  

    post_R = np.append(post_R, post_R[-1])

    # Local_R: AVG from last 10 pre_R values

    for i in range(0, len(R_poses)):

        num = 0

        avg_val = 0

        for j in range(-9, 1):

            if j+i >= 0:

                avg_val = avg_val + pre_R[i+j]

                num = num +1

        local_R = np.append(local_R, avg_val / float(num))

    # Global R AVG: from full past-signal

    # TODO: AVG from past 5 minutes = 108000 samples

    global_R = np.append(global_R, pre_R[0])    

    for i in range(1, len(R_poses)):

        num = 0

        avg_val = 0

        for j in range( 0, i):

            if (R_poses[i] - R_poses[j]) < 108000:

                avg_val = avg_val + pre_R[j]

                num = num + 1

        #num = i

        global_R = np.append(global_R, avg_val / float(num))

    for i in range(0, len(R_poses)):

        features_RR.pre_R = np.append(features_RR.pre_R, pre_R[i])

        features_RR.post_R = np.append(features_RR.post_R, post_R[i])

        features_RR.local_R = np.append(features_RR.local_R, local_R[i])

        features_RR.global_R = np.append(features_RR.global_R, global_R[i])

        #features_RR.append([pre_R[i], post_R[i], local_R[i], global_R[i]])

    return features_RR

# Compute the wavelet descriptor for a beat

def compute_wavelet_descriptor(beat, family, level):

    wave_family = pywt.Wavelet(family)

    coeffs = pywt.wavedec(beat, wave_family, level=level)

    return coeffs[0]

# Compute my descriptor based on amplitudes of several intervals

def compute_my_own_descriptor(beat, winL, winR):

    R_pos = int((winL + winR) / 2)

    R_value = beat[R_pos]

    my_morph = np.zeros((4))

    y_values = np.zeros(4)

    x_values = np.zeros(4)

    # Obtain (max/min) values and index from the intervals

    [x_values[0], y_values[0]] = max(enumerate(beat[0:40]), key=operator.itemgetter(1))

    [x_values[1], y_values[1]] = min(enumerate(beat[75:85]), key=operator.itemgetter(1))

    [x_values[2], y_values[2]] = min(enumerate(beat[95:105]), key=operator.itemgetter(1))

    [x_values[3], y_values[3]] = max(enumerate(beat[150:180]), key=operator.itemgetter(1))

    x_values[1] = x_values[1] + 75

    x_values[2] = x_values[2] + 95

    x_values[3] = x_values[3] + 150

    # Norm data before compute distance

    x_max = max(x_values)

    y_max = max(np.append(y_values, R_value))

    x_min = min(x_values)

    y_min = min(np.append(y_values, R_value))

    R_pos = (R_pos - x_min) / (x_max - x_min)

    R_value = (R_value - y_min) / (y_max - y_min)

    for n in range(0,4):

        x_values[n] = (x_values[n] - x_min) / (x_max - x_min)

        y_values[n] = (y_values[n] - y_min) / (y_max - y_min)

        x_diff = (R_pos - x_values[n]) 

        y_diff = R_value - y_values[n]

        my_morph[n] =  np.linalg.norm([x_diff, y_diff])

        # TODO test with np.sqrt(np.dot(x_diff, y_diff))

    if np.isnan(my_morph[n]):

        my_morph[n] = 0.0

    return my_morph

# Compute the HOS descriptor for a beat

# Skewness (3 cumulant) and kurtosis (4 cumulant)

def compute_hos_descriptor(beat, n_intervals, lag):

    hos_b = np.zeros(( (n_intervals-1) * 2))

    for i in range(0, n_intervals-1):

        pose = (lag * (i+1))

        interval = beat[(pose -(lag/2) ):(pose + (lag/2))]

        # Skewness  

        hos_b[i] = scipy.stats.skew(interval, 0, True)

        if np.isnan(hos_b[i]):

            hos_b[i] = 0.0

        # Kurtosis

        hos_b[(n_intervals-1) +i] = scipy.stats.kurtosis(interval, 0, False, True)

        if np.isnan(hos_b[(n_intervals-1) +i]):

            hos_b[(n_intervals-1) +i] = 0.0

    return hos_b

uniform_pattern_list = np.array([0, 1, 2, 3, 4, 6, 7, 8, 12, 14, 15, 16, 24, 28, 30, 31, 32, 48, 56, 60, 62, 63, 64, 96, 112, 120, 124, 126, 127, 128,

     129, 131, 135, 143, 159, 191, 192, 193, 195, 199, 207, 223, 224, 225, 227, 231, 239, 240, 241, 243, 247, 248, 249, 251, 252, 253, 254, 255])

# Compute the uniform LBP 1D from signal with neigh equal to number of neighbours

# and return the 59 histogram:

# 0-57: uniform patterns

# 58: the non uniform pattern

# NOTE: this method only works with neigh = 8

def compute_Uniform_LBP(signal, neigh=8):

    hist_u_lbp = np.zeros(59, dtype=float)

    avg_win_size = 2

    # NOTE: Reduce sampling by half

    #signal_avg = scipy.signal.resample(signal, len(signal) / avg_win_size)

    for i in range(neigh/2, len(signal) - neigh/2):

        pattern = np.zeros(neigh)

        ind = 0

        for n in range(-neigh/2,0) + range(1,neigh/2+1):

            if signal[i] > signal[i+n]:

                pattern[ind] = 1          

            ind += 1

        # Convert pattern to id-int 0-255 (for neigh == 8)

        pattern_id = int("".join(str(c) for c in pattern.astype(int)), 2)

        # Convert id to uniform LBP id 0-57 (uniform LBP)  58: (non uniform LBP)

        if pattern_id in uniform_pattern_list:

            pattern_uniform_id = int(np.argwhere(uniform_pattern_list == pattern_id))

        else:

            pattern_uniform_id = 58 # Non uniforms patternsuse

        hist_u_lbp[pattern_uniform_id] += 1.0

    return hist_u_lbp

def compute_LBP(signal, neigh=4):

    hist_u_lbp = np.zeros(np.power(2, neigh), dtype=float)

    avg_win_size = 2

    # TODO: use some type of average of the data instead the full signal...

    # Average window-5 of the signal?

    #signal_avg = average_signal(signal, avg_win_size)

    signal_avg = scipy.signal.resample(signal, len(signal) / avg_win_size)

    for i in range(neigh/2, len(signal) - neigh/2):

        pattern = np.zeros(neigh)

        ind = 0

        for n in range(-neigh/2,0) + range(1,neigh/2+1):

            if signal[i] > signal[i+n]:

                pattern[ind] = 1          

            ind += 1

        # Convert pattern to id-int 0-255 (for neigh == 8)

        pattern_id = int("".join(str(c) for c in pattern.astype(int)), 2)

        hist_u_lbp[pattern_id] += 1.0

    return hist_u_lbp

# https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.polynomials.hermite.html

# Support Vector Machine-Based Expert System for Reliable Heartbeat Recognition

# 15 hermite coefficients!

def compute_HBF(beat):

    coeffs_hbf = np.zeros(15, dtype=float)

    coeffs_HBF_3 = hermfit(range(0,len(beat)), beat, 3) # 3, 4, 5, 6?

    coeffs_HBF_4 = hermfit(range(0,len(beat)), beat, 4)

    coeffs_HBF_5 = hermfit(range(0,len(beat)), beat, 5)

    #coeffs_HBF_6 = hermfit(range(0,len(beat)), beat, 6)

    coeffs_hbf = np.concatenate((coeffs_HBF_3, coeffs_HBF_4, coeffs_HBF_5))

    return coeffs_hbf