#!/usr/bin/env python
import scipy.io
import numpy as np
import pywt
from scipy import signal
from biosppy.signals import ecg

# Global variable
Fs = 300

def mad(a, axis=None):
    """
    Compute *Median Absolute Deviation* of an array along given axis.
    """

    # Median along given axis, but *keeping* the reduced axis so that
    # result can still broadcast against a.
    med = np.median(a, axis=axis, keepdims=True)
    mad = np.median(np.absolute(a - med), axis=axis)  # MAD along given axis
    return mad


def ReadRecordMat(record):
    sample_len = 8192
    file_path = record + '.mat'
    mat = scipy.io.loadmat(file_path)
    ECG = signal.detrend(mat['val'][0])
    ECG = signal.resample(ECG,sample_len)
    return ECG


def HeartRate(ECG):
    Fs = int(len(ECG) / 30)
    try:
        out = ecg.ecg(signal=ECG, sampling_rate=Fs, show=False)
    except:
        out = ecg.christov_segmenter(signal=ECG,sampling_rate=Fs)
    r_peaks = out['rpeaks'] / float(Fs)
    HR = list()
    iternum = len(r_peaks)
    for idx in range(iternum):
        if idx > 0:
            prev_r = r_peaks[idx-1]
            current_r = r_peaks[idx]
            HR.append(current_r-prev_r)
    return HR

def PreHeartRate(HR):
    pre_HR = list()
    iternum = len(HR)
    for idx in range(iternum):
        if idx > 0:
            prev = HR[idx-1]
            cur = HR[idx]
            pre_HR.append(cur-prev)
    return pre_HR

def LogRatioHR(HR):
    iternum = len(HR)
    LogHR = list()
    for idx in range(iternum):
        if idx > 0:
            prev = np.log(HR[idx-1])
            cur = np.log(HR[idx])
            LogHR.append(cur-prev)
    return LogHR


def MedianTemplateECGsWCs(ECG,sample_temp_len = 256, wc_constant=3):
    Fs = int(len(ECG) / 30)
    try:
        out = ecg.ecg(signal=ECG, sampling_rate=Fs, show=False)
    except:
        out_RR = ecg.christov_segmenter(signal=ECG,sampling_rate=Fs)
        RR = out_RR['rpeaks']
        out = ecg.extract_heartbeats(signal=ECG, sampling_rate = Fs, rpeaks=RR)
    templates = out['templates']
    med_ECG = np.median(templates,axis=0)
    med_ECG = signal.resample(med_ECG,sample_temp_len)
    wc = WaveDecomp(med_ECG,dec_lvl=4)
    wc = WaveDenoise(wc,sample_temp_len,wc_constant)
    wc = RemoveNested(wc)
    return wc

def WaveDecomp(ECG, dec_lvl, wav_fun = 'db8'):
    new_len = int(2** np.ceil(np.log2(len(ECG))))
    ECG = signal.resample(ECG,new_len)
    wc = pywt.wavedec(data=ECG,wavelet=pywt.Wavelet(wav_fun),mode='per',level=dec_lvl)
    return wc

def WaveDenoise(wc,sample_len, constant=1):
    sigma = mad(wc[-1])
    uthresh = constant * sigma*np.sqrt(2*np.log(sample_len))
    denoised = wc[:]
    denoised[1:] = (pywt.threshold(i, uthresh,'soft') for i in denoised[1:])
    return denoised

def WaveRecov(wc, wav_fun = 'db8'):
    return pywt.waverec(wc,wav_fun,mode='per')

def OrderStat(myarray, ordernum):
    N = len(myarray)
    myarray = np.asarray(myarray)
    order_list = list()
    mean_val = np.mean(myarray)
    order_list.append(mean_val)
    if ordernum == 1:
        return order_list
    order_list.append(np.var(myarray))
    if ordernum == 2:
        return order_list

    if ordernum - 2 > 0:
        iternum = ordernum - 2
        for idx in range(iternum):
            val = np.sum((np.array(myarray) - mean_val)**(idx + 3))/float(N)
            order_list.append(val)
        return order_list

def RemoveNested(nested_list):
    return_list = list()
    for inner_list in nested_list:
        for elem in inner_list:
            return_list.append(elem)
    return return_list

