--- a
+++ b/python/load_MITBIH.py
@@ -0,0 +1,600 @@
+#!/usr/bin/env python
+
+"""
+load_MITBIH.py
+
+Download .csv files and annotations from:
+    kaggle.com/mondejar/mitbih-database
+
+VARPA, University of Coruna
+Mondejar Guerra, Victor M.
+23 Oct 2017
+"""
+
+import os
+import csv
+import gc
+import cPickle as pickle
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import medfilt
+import scipy.stats
+import pywt
+import time
+import sklearn
+from sklearn import decomposition
+from sklearn.decomposition import PCA, IncrementalPCA
+
+from features_ECG import *
+
+from numpy.polynomial.hermite import hermfit, hermval
+
+
+
+def create_features_labels_name(DS, winL, winR, do_preprocess, maxRR, use_RR, norm_RR, compute_morph, db_path, reduced_DS, leads_flag):
+    
+    features_labels_name = db_path + 'features/' + 'w_' + str(winL) + '_' + str(winR) + '_' + DS 
+
+    if do_preprocess:
+        features_labels_name += '_rm_bsline'
+
+    if maxRR:
+        features_labels_name += '_maxRR'
+
+    if use_RR:
+        features_labels_name += '_RR'
+    
+    if norm_RR:
+        features_labels_name += '_norm_RR'
+
+    for descp in compute_morph:
+        features_labels_name += '_' + descp
+
+    if reduced_DS:
+        features_labels_name += '_reduced'
+        
+    if leads_flag[0] == 1:
+        features_labels_name += '_MLII'
+
+    if leads_flag[1] == 1:
+        features_labels_name += '_V1'
+
+    features_labels_name += '.p'
+
+    return features_labels_name
+
+
+def save_wvlt_PCA(PCA, pca_k, family, level):
+    f = open('Wvlt_' + family + '_' + str(level) + '_PCA_' + str(pca_k) + '.p', 'wb')
+    pickle.dump(PCA, f, 2)
+    f.close
+
+
+def load_wvlt_PCA(pca_k, family, level):
+    f = open('Wvlt_' + family + '_' + str(level) + '_PCA_' + str(pca_k) + '.p', 'rb')
+    # disable garbage collector       
+    gc.disable()# this improve the required loading time!
+    PCA = pickle.load(f)
+    gc.enable()
+    f.close()
+
+    return PCA
+# Load the data with the configuration and features selected
+# Params:
+# - leads_flag = [MLII, V1] set the value to 0 or 1 to reference if that lead is used
+# - reduced_DS = load DS1, DS2 patients division (Chazal) or reduced version, 
+#                i.e., only patients in common that contains both MLII and V1 
+def load_mit_db(DS, winL, winR, do_preprocess, maxRR, use_RR, norm_RR, compute_morph, db_path, reduced_DS, leads_flag):
+
+    features_labels_name = create_features_labels_name(DS, winL, winR, do_preprocess, maxRR, use_RR, norm_RR, compute_morph, db_path, reduced_DS, leads_flag) 
+
+    if os.path.isfile(features_labels_name):
+        print("Loading pickle: " + features_labels_name + "...")
+        f = open(features_labels_name, 'rb')
+        # disable garbage collector       
+        gc.disable()# this improve the required loading time!
+        features, labels, patient_num_beats = pickle.load(f)
+        gc.enable()
+        f.close()
+
+
+    else:
+        print("Loading MIT BIH arr (" + DS + ") ...")
+
+        # ML-II
+        if reduced_DS == False:
+            DS1 = [101, 106, 108, 109, 112, 114, 115, 116, 118, 119, 122, 124, 201, 203, 205, 207, 208, 209, 215, 220, 223, 230]
+            DS2 = [100, 103, 105, 111, 113, 117, 121, 123, 200, 202, 210, 212, 213, 214, 219, 221, 222, 228, 231, 232, 233, 234]
+
+        # ML-II + V1
+        else:
+            DS1 = [101, 106, 108, 109, 112, 115, 118, 119, 201, 203, 205, 207, 208, 209, 215, 220, 223, 230]
+            DS2 = [105, 111, 113, 121, 200, 202, 210, 212, 213, 214, 219, 221, 222, 228, 231, 232, 233, 234]
+
+        mit_pickle_name = db_path + 'python_mit'
+        if reduced_DS:
+            mit_pickle_name = mit_pickle_name + '_reduced_'
+
+        if do_preprocess:
+            mit_pickle_name = mit_pickle_name + '_rm_bsline'
+
+        mit_pickle_name = mit_pickle_name + '_wL_' + str(winL) + '_wR_' + str(winR)
+        mit_pickle_name = mit_pickle_name + '_' + DS + '.p'
+
+        # If the data with that configuration has been already computed Load pickle
+        if os.path.isfile(mit_pickle_name):
+            f = open(mit_pickle_name, 'rb')
+            # disable garbage collector       
+            gc.disable()# this improve the required loading time!
+            my_db = pickle.load(f)
+            gc.enable()
+            f.close()
+        else: # Load data and compute de RR features 
+            if DS == 'DS1':
+                my_db = load_signal(DS1, winL, winR, do_preprocess)
+            else:
+                my_db = load_signal(DS2, winL, winR, do_preprocess)
+
+            print("Saving signal processed data ...")
+            # Save data
+            # Protocol version 0 itr_features_balanceds the original ASCII protocol and is backwards compatible with earlier versions of Python.
+            # Protocol version 1 is the old binary format which is also compatible with earlier versions of Python.
+            # Protocol version 2 was introduced in Python 2.3. It provides much more efficient pickling of new-style classes.
+            f = open(mit_pickle_name, 'wb')
+            pickle.dump(my_db, f, 2)
+            f.close
+
+        
+        features = np.array([], dtype=float)
+        labels = np.array([], dtype=np.int32)
+
+        # This array contains the number of beats for each patient (for cross_val)
+        patient_num_beats = np.array([], dtype=np.int32)
+        for p in range(len(my_db.beat)):
+            patient_num_beats = np.append(patient_num_beats, len(my_db.beat[p]))
+
+        # Compute RR features
+        if use_RR or norm_RR:
+            if DS == 'DS1':
+                RR = [RR_intervals() for i in range(len(DS1))]
+            else:
+                RR = [RR_intervals() for i in range(len(DS2))]
+
+            print("Computing RR intervals ...")
+
+            for p in range(len(my_db.beat)):
+                if maxRR:
+                    RR[p] = compute_RR_intervals(my_db.R_pos[p])
+                else:
+                    RR[p] = compute_RR_intervals(my_db.orig_R_pos[p])
+                    
+                RR[p].pre_R = RR[p].pre_R[(my_db.valid_R[p] == 1)]
+                RR[p].post_R = RR[p].post_R[(my_db.valid_R[p] == 1)]
+                RR[p].local_R = RR[p].local_R[(my_db.valid_R[p] == 1)]
+                RR[p].global_R = RR[p].global_R[(my_db.valid_R[p] == 1)]
+
+
+        if use_RR:
+            f_RR = np.empty((0,4))
+            for p in range(len(RR)):
+                row = np.column_stack((RR[p].pre_R, RR[p].post_R, RR[p].local_R, RR[p].global_R))
+                f_RR = np.vstack((f_RR, row))
+
+            features = np.column_stack((features, f_RR)) if features.size else f_RR
+        
+        if norm_RR:
+            f_RR_norm = np.empty((0,4))
+            for p in range(len(RR)):
+                # Compute avg values!
+                avg_pre_R = np.average(RR[p].pre_R)
+                avg_post_R = np.average(RR[p].post_R)
+                avg_local_R = np.average(RR[p].local_R)
+                avg_global_R = np.average(RR[p].global_R)
+
+                row = np.column_stack((RR[p].pre_R / avg_pre_R, RR[p].post_R / avg_post_R, RR[p].local_R / avg_local_R, RR[p].global_R / avg_global_R))
+                f_RR_norm = np.vstack((f_RR_norm, row))
+
+            features = np.column_stack((features, f_RR_norm))  if features.size else f_RR_norm
+
+        #########################################################################################
+        # Compute morphological features
+        print("Computing morphological features (" + DS + ") ...")
+
+        num_leads = np.sum(leads_flag)
+
+        # Raw
+        if 'resample_10' in compute_morph:
+            print("Resample_10 ...")
+            start = time.time()
+
+            f_raw = np.empty((0, 10 * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for beat in my_db.beat[p]:
+                    f_raw_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            resamp_beat = scipy.signal.resample(beat[s], 10)
+                            if f_raw_lead.size == 1:
+                                f_raw_lead =  resamp_beat
+                            else:
+                                f_raw_lead = np.hstack((f_raw_lead, resamp_beat))
+                    f_raw = np.vstack((f_raw, f_raw_lead))
+
+            features = np.column_stack((features, f_raw))  if features.size else f_raw
+
+            end = time.time()
+            print("Time resample: " + str(format(end - start, '.2f')) + " sec" )
+
+        if 'raw' in compute_morph:
+            print("Raw ...")
+            start = time.time()
+
+            f_raw = np.empty((0, (winL + winR) * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for beat in my_db.beat[p]:
+                    f_raw_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_raw_lead.size == 1:
+                                f_raw_lead =  beat[s]
+                            else:
+                                f_raw_lead = np.hstack((f_raw_lead, beat[s]))
+                    f_raw = np.vstack((f_raw, f_raw_lead))
+
+            features = np.column_stack((features, f_raw))  if features.size else f_raw
+
+            end = time.time()
+            print("Time raw: " + str(format(end - start, '.2f')) + " sec" )
+        # LBP 1D
+        # 1D-local binary pattern based feature extraction for classification of epileptic EEG signals: 2014, unas 55 citas, Q2-Q1 Matematicas
+        # https://ac.els-cdn.com/S0096300314008285/1-s2.0-S0096300314008285-main.pdf?_tid=8a8433a6-e57f-11e7-98ec-00000aab0f6c&acdnat=1513772341_eb5d4d26addb6c0b71ded4fd6cc23ed5
+
+        # 1D-LBP method, which derived from implementation steps of 2D-LBP, was firstly proposed by Chatlani et al. for detection of speech signals that is non-stationary in nature [23]
+
+        # From Raw signal
+
+        # TODO: Some kind of preprocesing or clean high frequency noise?
+
+        # Compute 2 Histograms: LBP or Uniform LBP
+        # LBP 8 = 0-255
+        # U-LBP 8 = 0-58
+        # Uniform LBP are only those pattern wich only presents 2 (or less) transitions from 0-1 or 1-0
+        # All the non-uniform patterns are asigned to the same value in the histogram
+
+        if 'u-lbp' in compute_morph:
+            print("u-lbp ...")
+
+            f_lbp = np.empty((0, 59 * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for beat in my_db.beat[p]:
+                    f_lbp_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_lbp_lead.size == 1:
+
+                                f_lbp_lead = compute_Uniform_LBP(beat[s], 8)
+                            else:
+                                f_lbp_lead = np.hstack((f_lbp_lead, compute_Uniform_LBP(beat[s], 8)))
+                    f_lbp = np.vstack((f_lbp, f_lbp_lead))
+
+            features = np.column_stack((features, f_lbp))  if features.size else f_lbp
+            print(features.shape)
+
+        if 'lbp' in compute_morph:
+            print("lbp ...")
+
+            f_lbp = np.empty((0, 16 * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for beat in my_db.beat[p]:
+                    f_lbp_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_lbp_lead.size == 1:
+
+                                f_lbp_lead = compute_LBP(beat[s], 4)
+                            else:
+                                f_lbp_lead = np.hstack((f_lbp_lead, compute_LBP(beat[s], 4)))
+                    f_lbp = np.vstack((f_lbp, f_lbp_lead))
+
+            features = np.column_stack((features, f_lbp))  if features.size else f_lbp
+            print(features.shape)
+
+
+        if 'hbf5' in compute_morph:
+            print("hbf ...")
+
+            f_hbf = np.empty((0, 15 * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for beat in my_db.beat[p]:
+                    f_hbf_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_hbf_lead.size == 1:
+
+                                f_hbf_lead = compute_HBF(beat[s])
+                            else:
+                                f_hbf_lead = np.hstack((f_hbf_lead, compute_HBF(beat[s])))
+                    f_hbf = np.vstack((f_hbf, f_hbf_lead))
+
+            features = np.column_stack((features, f_hbf))  if features.size else f_hbf
+            print(features.shape)
+
+        # Wavelets
+        if 'wvlt' in compute_morph:
+            print("Wavelets ...")
+
+            f_wav = np.empty((0, 23 * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for b in my_db.beat[p]:
+                    f_wav_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_wav_lead.size == 1:
+                                f_wav_lead =  compute_wavelet_descriptor(b[s], 'db1', 3)
+                            else:
+                                f_wav_lead = np.hstack((f_wav_lead, compute_wavelet_descriptor(b[s], 'db1', 3)))
+                    f_wav = np.vstack((f_wav, f_wav_lead))
+                    #f_wav = np.vstack((f_wav, compute_wavelet_descriptor(b,  'db1', 3)))
+
+            features = np.column_stack((features, f_wav))  if features.size else f_wav
+
+
+        # Wavelets
+        if 'wvlt+pca' in compute_morph:
+            pca_k = 7
+            print("Wavelets + PCA ("+ str(pca_k) + "...")
+            
+            family = 'db1'
+            level = 3
+
+            f_wav = np.empty((0, 23 * num_leads))
+
+            for p in range(len(my_db.beat)):
+                for b in my_db.beat[p]:
+                    f_wav_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_wav_lead.size == 1:
+                                f_wav_lead =  compute_wavelet_descriptor(b[s], family, level)
+                            else:
+                                f_wav_lead = np.hstack((f_wav_lead, compute_wavelet_descriptor(b[s], family, level)))
+                    f_wav = np.vstack((f_wav, f_wav_lead))
+                    #f_wav = np.vstack((f_wav, compute_wavelet_descriptor(b,  'db1', 3)))
+
+
+            if DS == 'DS1':
+                # Compute PCA
+                #PCA = sklearn.decomposition.KernelPCA(pca_k) # gamma_pca
+                IPCA = IncrementalPCA(n_components=pca_k, batch_size=10)# NOTE: due to memory errors, we employ IncrementalPCA
+                IPCA.fit(f_wav) 
+
+                # Save PCA
+                save_wvlt_PCA(IPCA, pca_k, family, level)
+            else:
+                # Load PCAfrom sklearn.decomposition import PCA, IncrementalPCA
+                IPCA = load_wvlt_PCA( pca_k, family, level)
+            # Extract the PCA
+            #f_wav_PCA = np.empty((0, pca_k * num_leads))
+            f_wav_PCA = IPCA.transform(f_wav)
+            features = np.column_stack((features, f_wav_PCA))  if features.size else f_wav_PCA
+
+        # HOS
+        if 'HOS' in compute_morph:
+            print("HOS ...")
+            n_intervals = 6
+            lag = int(round( (winL + winR )/ n_intervals))
+
+            f_HOS = np.empty((0, (n_intervals-1) * 2 * num_leads))
+            for p in range(len(my_db.beat)):
+                for b in my_db.beat[p]:
+                    f_HOS_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_HOS_lead.size == 1:
+                                f_HOS_lead =  compute_hos_descriptor(b[s], n_intervals, lag)
+                            else:
+                                f_HOS_lead = np.hstack((f_HOS_lead, compute_hos_descriptor(b[s], n_intervals, lag)))
+                    f_HOS = np.vstack((f_HOS, f_HOS_lead))
+                    #f_HOS = np.vstack((f_HOS, compute_hos_descriptor(b, n_intervals, lag)))
+
+            features = np.column_stack((features, f_HOS))  if features.size else f_HOS
+            print(features.shape)
+
+        # My morphological descriptor
+        if 'myMorph' in compute_morph:
+            print("My Descriptor ...")
+            f_myMorhp = np.empty((0,4 * num_leads))
+            for p in range(len(my_db.beat)):
+                for b in my_db.beat[p]:
+                    f_myMorhp_lead = np.empty([])
+                    for s in range(2):
+                        if leads_flag[s] == 1:
+                            if f_myMorhp_lead.size == 1:
+                                f_myMorhp_lead =  compute_my_own_descriptor(b[s], winL, winR)
+                            else:
+                                f_myMorhp_lead = np.hstack((f_myMorhp_lead, compute_my_own_descriptor(b[s], winL, winR)))
+                    f_myMorhp = np.vstack((f_myMorhp, f_myMorhp_lead))
+                    #f_myMorhp = np.vstack((f_myMorhp, compute_my_own_descriptor(b, winL, winR)))
+                    
+            features = np.column_stack((features, f_myMorhp))  if features.size else f_myMorhp
+
+        
+        labels = np.array(sum(my_db.class_ID, [])).flatten()
+        print("labels")
+
+        # Set labels array!
+        print('writing pickle: ' +  features_labels_name + '...')
+        f = open(features_labels_name, 'wb')
+        pickle.dump([features, labels, patient_num_beats], f, 2)
+        f.close
+
+    return features, labels, patient_num_beats
+
+
+# DS: contains the patient list for load
+# winL, winR: indicates the size of the window centred at R-peak at left and right side
+# do_preprocess: indicates if preprocesing of remove baseline on signal is performed
+def load_signal(DS, winL, winR, do_preprocess):
+
+    class_ID = [[] for i in range(len(DS))]
+    beat = [[] for i in range(len(DS))] # record, beat, lead
+    R_poses = [ np.array([]) for i in range(len(DS))]
+    Original_R_poses = [ np.array([]) for i in range(len(DS))]   
+    valid_R = [ np.array([]) for i in range(len(DS))]
+    my_db = mit_db()
+    patients = []
+
+    # Lists 
+    # beats = []
+    # classes = []
+    # valid_R = np.empty([])
+    # R_poses = np.empty([])
+    # Original_R_poses = np.empty([])
+
+    size_RR_max = 20
+
+    pathDB = '/home/mondejar/dataset/ECG/'
+    DB_name = 'mitdb'
+    fs = 360
+    jump_lines = 1
+
+    # Read files: signal (.csv )  annotations (.txt)    
+    fRecords = list()
+    fAnnotations = list()
+
+    lst = os.listdir(pathDB + DB_name + "/csv")
+    lst.sort()
+    for file in lst:
+        if file.endswith(".csv"):
+            if int(file[0:3]) in DS:
+                fRecords.append(file)
+        elif file.endswith(".txt"):
+            if int(file[0:3]) in DS:
+                fAnnotations.append(file)        
+
+    MITBIH_classes = ['N', 'L', 'R', 'e', 'j', 'A', 'a', 'J', 'S', 'V', 'E', 'F']#, 'P', '/', 'f', 'u']
+    AAMI_classes = []
+    AAMI_classes.append(['N', 'L', 'R'])                    # N
+    AAMI_classes.append(['A', 'a', 'J', 'S', 'e', 'j'])     # SVEB 
+    AAMI_classes.append(['V', 'E'])                         # VEB
+    AAMI_classes.append(['F'])                              # F
+    #AAMI_classes.append(['P', '/', 'f', 'u'])              # Q
+
+    RAW_signals = []
+    r_index = 0
+
+    #for r, a in zip(fRecords, fAnnotations):
+    for r in range(0, len(fRecords)):
+
+        print("Processing signal " + str(r) + " / " + str(len(fRecords)) + "...")
+
+        # 1. Read signalR_poses
+        filename = pathDB + DB_name + "/csv/" + fRecords[r]
+        print filename
+        f = open(filename, 'rb')
+        reader = csv.reader(f, delimiter=',')
+        next(reader) # skip first line!
+        MLII_index = 1
+        V1_index = 2
+        if int(fRecords[r][0:3]) == 114:
+            MLII_index = 2
+            V1_index = 1
+
+        MLII = []
+        V1 = []
+        for row in reader:
+            MLII.append((int(row[MLII_index])))
+            V1.append((int(row[V1_index])))
+        f.close()
+
+
+        RAW_signals.append((MLII, V1)) ## NOTE a copy must be created in order to preserve the original signal
+        # display_signal(MLII)
+
+        # 2. Read annotations
+        filename = pathDB + DB_name + "/csv/" + fAnnotations[r]
+        print filename
+        f = open(filename, 'rb')
+        next(f) # skip first line!
+
+        annotations = []
+        for line in f:
+            annotations.append(line)
+        f.close
+        # 3. Preprocessing signal!
+        if do_preprocess:
+            #scipy.signal
+            # median_filter1D
+            baseline = medfilt(MLII, 71) 
+            baseline = medfilt(baseline, 215) 
+
+            # Remove Baseline
+            for i in range(0, len(MLII)):
+                MLII[i] = MLII[i] - baseline[i]
+
+            # TODO Remove High Freqs
+
+            # median_filter1D
+            baseline = medfilt(V1, 71) 
+            baseline = medfilt(baseline, 215) 
+
+            # Remove Baseline
+            for i in range(0, len(V1)):
+                V1[i] = V1[i] - baseline[i]
+
+
+        # Extract the R-peaks from annotations
+        for a in annotations:
+            aS = a.split()
+            
+            pos = int(aS[1])
+            originalPos = int(aS[1])
+            classAnttd = aS[2]
+            if pos > size_RR_max and pos < (len(MLII) - size_RR_max):
+                index, value = max(enumerate(MLII[pos - size_RR_max : pos + size_RR_max]), key=operator.itemgetter(1))
+                pos = (pos - size_RR_max) + index
+
+            peak_type = 0
+            #pos = pos-1
+            
+            if classAnttd in MITBIH_classes:
+                if(pos > winL and pos < (len(MLII) - winR)):
+                    beat[r].append( (MLII[pos - winL : pos + winR], V1[pos - winL : pos + winR]))
+                    for i in range(0,len(AAMI_classes)):
+                        if classAnttd in AAMI_classes[i]:
+                            class_AAMI = i
+                            break #exit loop
+                    #convert class
+                    class_ID[r].append(class_AAMI)
+
+                    valid_R[r] = np.append(valid_R[r], 1)
+                else:
+                    valid_R[r] = np.append(valid_R[r], 0)
+            else:
+                valid_R[r] = np.append(valid_R[r], 0)
+            
+            R_poses[r] = np.append(R_poses[r], pos)
+            Original_R_poses[r] = np.append(Original_R_poses[r], originalPos)
+        
+        #R_poses[r] = R_poses[r][(valid_R[r] == 1)]
+        #Original_R_poses[r] = Original_R_poses[r][(valid_R[r] == 1)]
+
+        
+    # Set the data into a bigger struct that keep all the records!
+    my_db.filename = fRecords
+
+    my_db.raw_signal = RAW_signals
+    my_db.beat = beat # record, beat, lead
+    my_db.class_ID = class_ID
+    my_db.valid_R = valid_R
+    my_db.R_pos = R_poses
+    my_db.orig_R_pos = Original_R_poses
+
+    return my_db