text:
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 22 08:55:50 2019

@author: logiusti
"""

import pandas as pd
import numpy as np
from numpy.random import seed
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(color_codes=True)

# Common imports
import os
from sklearn import preprocessing
from sklearn.decomposition import KernelPCA
from sklearn.ensemble import IsolationForest

import tensorflow as tf
from tensorflow.keras.layers import Input, Lambda, LeakyReLU
from tensorflow.keras.layers import Dense, LSTM, RepeatVector, TimeDistributed, Bidirectional
from tensorflow.keras.models import Model, Sequential
from tensorflow.keras.losses import mse
from tensorflow.keras.callbacks import EarlyStopping

gpus = tf.config.experimental.list_physical_devices('GPU')
for gpu in gpus:
    tf.config.experimental.set_memory_growth(gpu, True)

from scipy.stats import wasserstein_distance

seed(10)
tf.random.set_seed(10)

import Loader
class Detector:
    r'''
    Anomaly detection (or outlier detection) is the identification of rare items,
    events or observations which raise suspicions by differing significantly from the majority of the data.
    Typically, anomalous data can be connected to some kind of problem or rare event such as:
        bank fraud, medical problems, structural defects, malfunctioning equipment etc.
    This connection makes it very interesting to be able to pick out which data points can be considered anomalies,
    as identifying these events are typically very interesting from a business perspective.

    In the case of two-dimensional data (X and Y), or at most three-dimensional data (X, Y and Z),
    it becomes quite easy to visually identify anomalies through data points located outside the typical distribution.
    However, often is not possible to identify the outlier directly from investigating one variable at the time:
        It is the combination of the X and Y (and Z if exist) variable that allows us to easily identify the anomaly.
    This complicates the matter substantially when we scale up from two variables to 10–100s of variables,
    which is often the case in practical applications of anomaly detection.

   '''

    def __init__(self):
        self.loader = Loader.Loader()

    def load_and_scale(self, data):
        r'''
        Before setting up the models,
        we need to define train/test data.
        To do this, we perform a simple split where we train on the 90% of the dataset
        (which should represent normal operating conditions),
        and test on the remaining parts of the dataset leading up to the eventual failure.
        '''
        test_size = np.ceil(len(data)*.15).astype(int)

        dataset_train = data[:-test_size]
        dataset_test = data[-test_size:]

        r'''
        Scale the input variables of the model.

        Standardize features by removing the mean and scaling to unit variance.
        Centering and scaling happen independently on each feature by computing the relevant statistics on the samples in the training set.
        Mean and standard deviation are then stored to be used on later data using the transform method.
        Standardization of a dataset is a common requirement for many machine learning estimators:
            they might behave badly if the individual features do not more or less look like standard normally distributed data (e.g. Gaussian with 0 mean and unit variance).
        For instance many elements used in the objective function of a learning algorithm (such as the RBF kernel of Support Vector Machines or the L1 and L2 regularizers of linear models) assume that all features are centered around 0 and have variance in the same order.
        If a feature has a variance that is orders of magnitude larger that others, it might dominate the objective function and make the estimator unable to learn from other features correctly as expected.
        '''
        scaler = preprocessing.MinMaxScaler()

        X_train = pd.DataFrame(scaler.fit_transform(dataset_train),
                               columns=dataset_train.columns,
                               index=dataset_train.index)#.sample(frac=1)# Random shuffle training data

        X_test = pd.DataFrame(scaler.transform(dataset_test),
                              columns=dataset_test.columns,
                              index=dataset_test.index)

        return (X_train, X_test)

    def scheduler(self):
        ups_to_functionals = self.loader.load_temperature_functionals()
        i = 0
        for ups in ups_to_functionals:
            i += 1
            if i % 20 == 0:
                print(i)
            if ups[0] == 'EBS11_SLASH_38':# or ups[0] not in ['ESS329_SLASH_7E']: #
                continue
            if not os.path.isdir(r'./data/anom/'+ups[0]):
                os.mkdir(r'./data/anom/'+ups[0])

            #X_train, X_test = self.load_and_scale(ups_to_functionals[ups])
            #md = self.mahalanobis_distances(X_train, X_test, ups)

            X_train, X_test = self.load_and_scale(ups[1])
            #autoenc = self.nonlinear_autoencoder_detect(X_train, X_test, ups)
            LSTM_autoenc = self.nonlinear_LSTM_autoencoder_detect(X_train, X_test, ups[0])

            fig, axes = plt.subplots(nrows=1)

            #md.plot(logy=True, figsize=(30, 18), ylim=[1e-3, 1e3], color=['green', 'red'], title='Mahalanobis Anomaly' ,ax=axes[0])
            #autoenc.plot(logy=True, figsize=(30, 18), ylim=[1e-3, 1e3], color=['green', 'red'], title='Autoencoder Anomaly' , ax=axes[1])
            LSTM_autoenc.plot(logy=True, figsize=(30, 18), ylim=[1e-3, 1e3], color=['green', 'red'],  title='LSTM Autoencoder Anomaly',  ax=axes)
            plt.savefig(r'./data/anom/'+ups[0]+'/merged.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')

            plt.close('all')

    def nonlinear_LSTM_autoencoder_detect(self, X_train, X_test, ups):

            def temporalize(X, lookback):
                output_X = []
                for i in range(len(X)-lookback-1):
                    t = []
                    for j in range(1,lookback+1):
                        # Gather past records upto the lookback period
                        t.append(X[[(i+j+1)]])
                    output_X.append(t)
                return output_X

            def flatten(X):
                '''
                Flatten a 3D array.

                Input
                X            A 3D array for lstm, where the array is sample x timesteps x features.

                Output
                flattened_X  A 2D array, sample x features.
                '''
                flattened_X = np.empty((X.shape[0], X.shape[2]))  # sample x features array.
                for i in range(X.shape[0]):
                    flattened_X[i] = X[i, (X.shape[1]-1), :]
                return(flattened_X)

            lookback = 3
            features = 6

            df = X_train.append(X_test, sort=False)
            t = df.index
            test_size = np.ceil(len(df)*.15).astype(int)

            frame = df.to_numpy()

            tmp =  np.apply_along_axis(temporalize, 0, frame[:,0], lookback=lookback)
            for i in range(1,features):
                tmp = np.append(tmp, np.apply_along_axis(temporalize, 0, frame[:,i], lookback=lookback),axis=2)

            train = tmp[:-test_size]
            test = tmp[-test_size:]

            #corrupted_train = train.copy()

            #random_rows_to_corrupt = np.random.randint(0, corrupted_train.shape[0], size=(int(corrupted_train.size * 0.05),))
            #random_cols_to_corrupt = np.random.randint(0, corrupted_train.shape[1], size=(int(corrupted_train.size * 0.05),))

            #corrupted_train[random_rows_to_corrupt, random_cols_to_corrupt] = 0

            model = Sequential()

            model.add(Bidirectional(LSTM(8, activation='relu',
                                         kernel_initializer='lecun_normal',
                                         #kernel_regularizer=tf.keras.regularizers.l2(0.001),
                                         input_shape=(lookback,features), return_sequences=True)))
            model.add(Bidirectional(LSTM(4, activation='relu',
                                         kernel_initializer='lecun_normal',
                                         #kernel_regularizer=tf.keras.regularizers.l2(0.001),
                                         return_sequences=False)))
            model.add(RepeatVector(lookback))
            model.add(Bidirectional(LSTM(4, activation='relu',
                                         kernel_initializer='lecun_normal',
                                         #kernel_regularizer=tf.keras.regularizers.l2(0.001),
                                         return_sequences=True)))
            model.add(Bidirectional(LSTM(8, activation='relu',
                                         kernel_initializer='lecun_normal',
                                         #kernel_regularizer=tf.keras.regularizers.l2(0.001),
                                         return_sequences=True)))
            model.add(TimeDistributed(Dense(features)))

            adam = tf.keras.optimizers.Adam(learning_rate=0.003, amsgrad=True)
            model.compile(optimizer=adam, loss='mse')
            #model.summary()
            # fit model
            es = EarlyStopping(monitor='val_loss', mode='min', min_delta=0.001, verbose=0, patience=15)
            batch_size = int(np.floor(len(df)*0.03))
            if batch_size == 0 :
                batch_size = 1
            history = model.fit(train+np.random.normal(0,0.1), train , epochs=150, batch_size=batch_size, verbose=1, steps_per_epoch=None, validation_split=0.05, callbacks=[es])

            plt.plot(history.history['loss'],
                                 'b',
                                 label='Training loss')
            plt.plot(history.history['val_loss'],
                     'r',
                     label='Validation loss')
            plt.legend(loc='upper right')
            plt.xlabel('Epochs')
            plt.ylabel('Loss, [mse]')
            ymax = max(max(history.history['loss']), max(history.history['val_loss']))+0.5
            ymin = min(min(history.history['loss']), min(history.history['val_loss']))-0.5
            plt.ylim([ymin, ymax])
            plt.savefig(r'./data/anom/'+ups+'/LSTM_Autoencoder_loss.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')

            X_pred = model.predict(train).reshape(train.shape[0], lookback, train.shape[-1])

            scored_train = pd.DataFrame(index= t[lookback+1:lookback+train.shape[0]+1])
            scored_train['Loss_mae'] = np.mean(flatten(np.abs(X_pred-train)), axis = 1)

            plt.figure()
            sns.distplot(scored_train['Loss_mae'],
                                     bins = 100,
                                     kde= True,
                                    color = 'darkblue');
            plt.xlim([min(scored_train['Loss_mae']), max(scored_train['Loss_mae'])])
            plt.savefig(r'./data/anom/'+ups+'/LSTM_Autoencoder_loss_dist_train.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
            plt.close()

            threshold = scored_train['Loss_mae'].quantile(.95) + 1*(scored_train['Loss_mae'].quantile(.95) - scored_train['Loss_mae'].quantile(.05))
            scored_train['Threshold'] = threshold
            scored_train['Anomaly'] = scored_train['Loss_mae'] > scored_train['Threshold']

            X_test_pred = model.predict(test).reshape(test.shape[0], lookback, test.shape[-1])
            scored_test = pd.DataFrame(index=t[lookback+train.shape[0]+1:])
            scored_test['Loss_mae'] = np.mean(flatten(np.abs(X_test_pred-test)), axis=1)
            scored_test['Threshold'] = threshold
            scored_test['Anomaly'] = scored_test['Loss_mae'] > scored_test['Threshold']

            plt.figure()
            sns.distplot(scored_test['Loss_mae'],
                                     bins =35,
                                     kde= True,
                                    color = 'darkorange');
            plt.xlim([min(scored_test['Loss_mae']), max(scored_test['Loss_mae'])])
            plt.savefig(r'./data/anom/'+ups+'/LSTM_Autoencoder_loss_dist_test.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
            plt.close()

            plt.figure()
            sns.distplot(scored_train['Loss_mae'],
                                     bins = 150,
                                     kde= True,
                                    color = 'darkblue',
                                    label='Train Density');
            sns.distplot(scored_test['Loss_mae'],
                                     bins = 35,
                                     kde= True,
                                    color = 'darkorange',
                                    label='Test Density');

            plt.axvline(x=threshold, alpha=0.88, color='black', linestyle='--', label='Threshold')
            lb = min(threshold, min(scored_test['Loss_mae']), min(scored_train['Loss_mae']))
            ub = max(threshold, max(scored_test['Loss_mae']), max(scored_train['Loss_mae']))
            plt.xlim([ lb*1.1, ub*1.1 ])
            plt.legend(loc="upper right")
            plt.savefig(r'./data/anom/'+ups+'/LSTM_Autoencoder_loss_dist_overlap.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
            plt.close()

            scored = pd.concat([scored_train, scored_test], sort=False)

            ymax = 100*max(scored['Loss_mae'])
            ymin = 0.001*min(scored['Loss_mae'])
            fig = scored.plot(logy=True, figsize=(15, 9), ylim=[ymin, ymax], color=['blue', 'red'])
            fig.get_figure().savefig(r'./data/anom/'+ups+'/LSTM_Autoencoder_test_threshold.jpeg',
                                     quality=95, optimize=True, progressive=True, format='jpeg')

            plt.close('all')

            scored.sort_index(inplace=True)
            scored.to_csv(r'./data/anom/'+ups+'/LSTM_autoencoder_distance.csv')

            return scored

    def performance_measurements(self):
        r"""
        True Positive (TP) if the device is labelled as anomalous and in the test period we observe a value that goes over the threshold and continue over it until the end of its period even with decreases.
    	False Positive (FP) if the device is labelled as anomalous and in the test period we don't observe values over the threshold or values over it for a limited period of time which doesn't include the ending of the time series which can represent an anomaly of the surrounding environment.
    	True Negative (TN) if the device is labelled as normal and in the test period we don't observe values over the threshold or values over it for a limited period of time which doesn't include the ending of the time series.
    	False Negative (FN) if the device is labelled as normal and in the test period we observe a value that goes over the threshold and continue over it until the end of its period even with decreases.
        """
        pass

    def variational_autoencoder(self):

        ups_to_functionals = self.loader.load_temperature_functionals()

        for ups in ups_to_functionals:
            if ups == 'EBS11_SLASH_38'  or len(ups_to_functionals[ups]) < 50:
                continue
            if not os.path.isdir(r'./data/anom/'+ups):
                os.mkdir(r'./data/anom/'+ups)

            r'''
            Before setting up the models,
            we need to define train/test data.
            To do this, we perform a simple split where we train on the 90% of the dataset
            (which should represent normal operating conditions),
            and test on the remaining parts of the dataset leading up to the eventual failure.
            '''
            data = ups_to_functionals[ups]
            test_size = np.ceil(len(data)*.1).astype(int)

            dataset_train = data[:-test_size]
            dataset_test = data[-test_size:]

            r'''
            Scale the input variables of the model.

            Standardize features by removing the mean and scaling to unit variance.
            Centering and scaling happen independently on each feature by computing the relevant statistics on the samples in the training set.
            Mean and standard deviation are then stored to be used on later data using the transform method.
            Standardization of a dataset is a common requirement for many machine learning estimators:
                they might behave badly if the individual features do not more or less look like standard normally distributed data (e.g. Gaussian with 0 mean and unit variance).
            For instance many elements used in the objective function of a learning algorithm (such as the RBF kernel of Support Vector Machines or the L1 and L2 regularizers of linear models) assume that all features are centered around 0 and have variance in the same order.
            If a feature has a variance that is orders of magnitude larger that others, it might dominate the objective function and make the estimator unable to learn from other features correctly as expected.
            '''
            scaler = preprocessing.StandardScaler()

            X_train = pd.DataFrame(scaler.fit_transform(dataset_train),
                                   columns=dataset_train.columns,
                                   index=dataset_train.index).sample(frac=1)# Random shuffle training data

            X_test = pd.DataFrame(scaler.transform(dataset_test),
                                  columns=dataset_test.columns,
                                  index=dataset_test.index)

            #act_func = 'relu'
            train = X_train.to_numpy()
            #test = X_test.to_numpy()
            intermediate_dim = 8
            batch_size = 64
            latent_dim = 2
            epochs = 150

            def sampling(args):
                """Reparameterization trick by sampling from an isotropic unit Gaussian.

                # Arguments
                    args (tensor): mean and log of variance of Q(z|X)

                # Returns
                    z (tensor): sampled latent vector
                """

                z_mean, z_log_var = args
                batch = tf.shape(z_mean)[0]
                dim = tf.shape(z_mean)[1]
                # by default, random_normal has mean = 0 and std = 1.0
                epsilon = tf.random.normal(shape=(batch, dim))
                return z_mean + tf.math.exp(0.5 * z_log_var) * epsilon

            # VAE model = encoder + decoder
            # build encoder model
            inputs = Input(shape=(train.shape[1],), name='encoder_input')
            x = Dense(intermediate_dim, activation='relu')(inputs)
            z_mean = Dense(latent_dim, name='z_mean')(x)
            z_log_var = Dense(latent_dim, name='z_log_var')(x)

            # use reparameterization trick to push the sampling out as input
            # note that "output_shape" isn't necessary with the TensorFlow backend
            z = Lambda(sampling, output_shape=(latent_dim,), name='z')([z_mean, z_log_var])

            # instantiate encoder model
            encoder = Model(inputs, [z_mean, z_log_var, z], name='encoder')

            # build decoder model
            latent_inputs = Input(shape=(latent_dim,), name='z_sampling')
            x = Dense(intermediate_dim, activation='relu')(latent_inputs)
            outputs = Dense(train.shape[1], activation='relu')(x)

            # instantiate decoder model
            decoder = Model(latent_inputs, outputs, name='decoder')

            # instantiate VAE model
            outputs = decoder(encoder(inputs)[2])
            vae = Model(inputs, outputs, name='vae_mlp')

            reconstruction_loss = mse(inputs, outputs)
            reconstruction_loss *= train.shape[1]
            kl_loss = 1 + z_log_var - tf.math.square(z_mean) - tf.math.exp(z_log_var)
            kl_loss = tf.reduce_sum(kl_loss, axis=-1)
            kl_loss *= -0.5
            vae_loss = tf.reduce_mean(reconstruction_loss + kl_loss)
            vae.add_loss(vae_loss)
            vae.compile(optimizer='adam')
            history = vae.fit(train,
                    epochs=epochs,
                    batch_size=batch_size,
                  validation_split=0.05,
                  verbose=0)

            # Visualize training/validation loss
            plt.plot(history.history['loss'], 'b', label='Training loss')
            plt.plot(history.history['val_loss'], 'r', label='Validation loss')
            plt.legend(loc='upper right')
            plt.xlabel('Epochs')
            plt.ylabel('Loss, [mse]')
            ylim = max(max(history.history['loss']), max(history.history['val_loss']))
            plt.ylim([0, ylim])
            plt.savefig(r'./data/anom/'+ups+'/VAE_loss.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
            plt.close()

            r'''
            Distribution of loss function in the training set:

            By plotting the distribution of the calculated loss in the training set,
            one can use this to identify a suitable threshold value for identifying an anomaly.
            In doing this, one can make sure that this threshold is set above the “noise level”,
            and that any flagged anomalies should be statistically significant above the noise background.
            '''
            X_pred = vae.predict(np.array(X_train))
            X_pred = pd.DataFrame(X_pred,
                                  columns=X_train.columns)
            X_pred.index = X_train.index

            scored = pd.DataFrame(index=X_train.index)
            scored['Loss_mae'] = np.mean(np.abs(X_pred-X_train), axis=1)
            plt.figure()
            sns.distplot(scored['Loss_mae'],
                         bins=50,
                         kde=True,
                         color='blue')
            plt.xlim([0.0, max(scored['Loss_mae'])])
            plt.savefig(r'./data/anom/'+ups+'/VAE_loss_dist.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
            plt.close()

            r'''
            From the above loss distribution,
            let us try a threshold based on one of the last quantiles for flagging an anomaly.
            We can then calculate the loss in the test set, to check when the output crosses the anomaly threshold.
            '''

            threshold = scored['Loss_mae'].quantile(.95) + 1*(scored['Loss_mae'].quantile(.95) - scored['Loss_mae'].quantile(.05))

            X_pred = vae.predict(np.array(X_test))
            X_pred = pd.DataFrame(X_pred,
                                  columns=X_test.columns)
            X_pred.index = X_test.index

            scored = pd.DataFrame(index=X_test.index)
            scored['Loss_mae'] = np.mean(np.abs(X_pred-X_test), axis=1)
            scored['Threshold'] = threshold
            scored['Anomaly'] = scored['Loss_mae'] > scored['Threshold']

            ymax = 100*max(scored['Loss_mae'])
            ymin = 0.001*min(scored['Loss_mae'])
            fig = scored.plot(logy=True, figsize=(30, 18), ylim=[ymin, ymax], color=['blue', 'red'])
            fig.get_figure().savefig(r'./data/anom/'+ups+'/VAE_test_threshold.jpeg',
                                     quality=95, optimize=True, progressive=True, format='jpeg')

            plt.close('all')

            r'''
            We then calculate the same metrics also for the training set, and merge all data in a single dataframe.
            '''
            X_pred = vae.predict(np.array(X_test))
            X_pred = pd.DataFrame(X_pred,
                                  columns=X_test.columns)
            X_pred.index = X_test.index

            scored_test = pd.DataFrame(index=X_test.index)
            scored_test['Loss_mae'] = np.mean(np.abs(X_pred-X_test), axis=1)
            scored_test['Threshold'] = threshold
            scored_test['Anomaly'] = scored_test['Loss_mae'] > scored_test['Threshold']

            X_pred_train = vae.predict(np.array(X_train))
            X_pred_train = pd.DataFrame(X_pred_train,
                                        columns=X_train.columns)
            X_pred_train.index = X_train.index

            scored_train = pd.DataFrame(index=X_train.index)
            scored_train['Loss_mae'] = np.mean(np.abs(X_pred_train-X_train), axis=1)
            scored_train['Threshold'] = threshold
            scored_train['Anomaly'] = scored_train['Loss_mae'] > scored_train['Threshold']

            scored = pd.concat([scored_train, scored_test], sort=False)

            r'''
            Isolation forest prediction based on anomality detected by the autoencoder
            '''
            contamination = (scored.Anomaly == True).sum()/len(scored)

            clf = IsolationForest(behaviour='new', n_estimators=500, max_samples='auto', contamination=contamination, \
                                    max_features=1.0, bootstrap=True, n_jobs=-1, random_state=42, verbose=0)

            clf.fit(scored.Loss_mae.to_numpy().reshape(-1, 1))

            pred = clf.predict(scored.Loss_mae.to_numpy().reshape(-1, 1))
            dct = {-1:threshold, 1:np.nan}
            pred = [dct[k] for k in pred]

            r'''
            Results from Autoencoder model:
            Having calculated the loss distribution and the anomaly threshold, we can visualize the model output in the time leading up to the bearing failure:
            '''
            scored.plot(logy=True, figsize=(30, 18), ylim=[1e-3, 1e3], color=['green', 'red'])
            isolation_forest = plt.scatter(scored.index, pred, s=150, c='b', marker=".")
            plt.legend((isolation_forest,), ('Isolation Forest',), scatterpoints=1, loc='upper right', ncol=1, fontsize=8)
            plt.savefig(r'./data/anom/'+ups+'/VAE_full_threshold.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
            plt.close('all')

    def false_positive_detect(self, crit_seqs_path=r"C:\Users\logiusti\Lorenzo\PyWorkspace\scripts\Wrapper\data\critical_dates",
                              anomalities_path=r'C:\Users\logiusti\Lorenzo\PyWorkspace\scripts\Wrapper\data\anom'):

        ups_to_false_positive = dict()
        def pick_the_peak(df):
            r'''
            We pick the time in which the loss is maximum in a t+-dt period.
            '''
            return df.loc[df['Loss_mae'] == max(df['Loss_mae'])]

        os.chdir(anomalities_path)
        for ups in os.listdir():
            ups_to_false_positive[ups] = []
            # read the critical sequences timestamps of each ups
            try:
                crit_seq = pd.to_datetime(pd.Series(pd.read_csv(crit_seqs_path+"\\"+ups+'_crit_seq_dates.csv', header=None)[0])).sort_values().reset_index(drop=True)
            except FileNotFoundError:
                crit_seq = pd.Series()
            #read the anomaly alterts, we choose the LSTM denoising autoencore since it seems to have better performances.
            anomalities = pd.read_csv("./"+ups+'/LSTM_autoencoder_distance.csv')
            anomalities = anomalities.loc[anomalities['Anomaly'] == True] #use only values flagged as anomaly
            anomalities['Time'] = pd.to_datetime(anomalities['Time'], format="%Y-%m-%d %H:%M:%S.%f")

            if anomalities.Time.empty:
                ups_to_false_positive[ups] = []
                continue
            #tim off, first 10% of the period should not be considered in the false positive count
            dayoff = ((anomalities.Time.iloc[-1] - anomalities.Time.iloc[0]).days)*.1
            #cut-off the DF
            anomalities = anomalities.loc[anomalities['Time']  >= (pd.to_datetime(anomalities.Time.iloc[0], format="%Y-%m-%d") + pd.DateOffset(days=dayoff))]

            #we want to group an entire period in which we've spotted an anomaly (7 days as upper bound)
            anomalities['groups'] = (anomalities.Time.diff().fillna(pd.Timedelta(seconds=0))/np.timedelta64(7, 'D'))\
                                            .gt(1).cumsum().add(1).astype(str)

            anomalities.set_index('groups', drop=True, inplace=True)
            danger = anomalities.groupby(level=0, group_keys=False).apply(pick_the_peak) #take the max of that period
            danger['Time'] = danger['Time'].dt.date

            #if nothing happens in the time(max[])+4days
            for idx,row in danger.iterrows():
                start = pd.to_datetime(row.Time, format="%Y-%m-%d")
                end = start + pd.DateOffset(days=4)
                if crit_seq.loc[(crit_seq >= start ) & (crit_seq <= end) ].empty and row['Loss_mae'] >= 1.5*row['Threshold']:
                    print("Falso Positivo", ups, "\n", row)
                    ups_to_false_positive[ups].append(row['Time'])
        return ups_to_false_positive

    def min_cost_transport(self, anomalities_path=r'C:\Users\logiusti\Lorenzo\PyWorkspace\scripts\Wrapper\data\anom'):
        os.chdir(anomalities_path)
        ups_to_dist = dict()
        for ups in os.listdir():
            df = pd.read_csv("./"+ups+'/LSTM_autoencoder_distance.csv')
            df['Time'] = pd.to_datetime(df['Time'], format="%Y-%m-%d %H:%M:%S.%f")
            df.sort_values(by='Time', inplace=True)
            df.reset_index(drop=True, inplace=True)
            test_size = np.ceil(len(df)*.5).astype(int)
            train = df[:-test_size]
            test = df[-test_size:]
            wa, a = np.histogram(test['Loss_mae'], density=True)
            wb, b = np.histogram(train['Loss_mae'], density=True)
            ups_to_dist[ups] = wasserstein_distance(a,b)

        return ups_to_dist

    def mahalanobis_distances(self, X_train, X_test, ups):
        r'''
        The Mahalanobis distance is widely used in cluster analysis and classification techniques.
        In order to use the Mahalanobis distance to classify a test point as belonging to one of N classes,
        one first estimates the covariance matrix of each class,
        usually based on samples known to belong to each class.
        In our case, as we are only interested in classifying “normal” vs “anomaly”,
        we use training data that only contains normal operating conditions to calculate the covariance matrix.
        Then, given a test sample, we compute the Mahalanobis distance to the “normal” class,
        and classifies the test point as an “anomaly” if the distance is above a certain threshold.
        '''

        def is_pos_def(A):
            r'''
            Check if matrix is positive definite.
            '''
            if np.allclose(A, A.T):
                try:
                    np.linalg.cholesky(A)
                    return True
                except np.linalg.LinAlgError:
                    return False
            else:
                return False

        def cov_matrix_f(data):
            r'''
            Calculate the covariance matrix.
            '''
            covariance_matrix = np.cov(data, rowvar=False)
            if is_pos_def(covariance_matrix):
                inv_covariance_matrix = np.linalg.inv(covariance_matrix)
                if is_pos_def(inv_covariance_matrix):
                    return covariance_matrix, inv_covariance_matrix
                else:
                    print("Error: Inverse of Covariance Matrix is not positive definite!")
            else:
                print("Error: Covariance Matrix is not positive definite!")

        def MahalanobisDist(inv_cov_matrix, mean_distr, data):
            r'''
            Calculate the Mahalanobis distance.
            '''
            inv_covariance_matrix = inv_cov_matrix
            vars_mean = mean_distr
            diff = data - vars_mean
            md = []
            for i in range(len(diff)):
                md.append(np.sqrt(diff[i].dot(inv_covariance_matrix).dot(diff[i])))
            return md

        def MD_detectOutliers(dist, extreme=False):
            r'''
            Detecting outliers.
            '''
            k = 3. if extreme else 2.
            threshold = np.mean(dist) * k
            outliers = []
            for i in range(len(dist)):
                if dist[i] >= threshold:
                    outliers.append(i)  # index of the outlier
            return np.array(outliers)

        def MD_threshold(dist, extreme=False):
            r'''
            Calculate threshold value for classifying datapoint as anomaly.
            '''
            k = 3. if extreme else 2.
            threshold = np.mean(dist) * k
            return threshold

        r'''
        As dealing with high dimensional sensor data is often challenging,
        there are several techniques to reduce the number of variables (dimensionality reduction).
        One of the main techniques is principal component analysis (PCA),
        which performs a linear mapping of the data to a lower-dimensional space
        in such a way that the variance of the data in the low-dimensional representation is maximized.
        In practice, the covariance matrix of the data is constructed and the eigenvectors of this matrix are computed.
        The eigenvectors that correspond to the largest eigenvalues (the principal components) can now be used to reconstruct a large fraction of the variance of the original data.
        The original feature space has now been reduced (with some data loss, but hopefully retaining the most important variance) to the space spanned by a few eigenvectors.

        We use Kernel Principal component analysis (KPCA) since we are almost sure that the patterns are NON-Linear.

        Non-linear dimensionality reduction through the use of kernels.
        '''

        pca = KernelPCA(n_components=2, kernel='poly', degree=3, coef0=0, fit_inverse_transform=False,
                        remove_zero_eig=True, n_jobs=-1, eigen_solver='auto')

        X_train_PCA = pca.fit_transform(X_train)
        X_train_PCA = pd.DataFrame(X_train_PCA)
        X_train_PCA.index = X_train.index

        X_test_PCA = pca.transform(X_test)
        X_test_PCA = pd.DataFrame(X_test_PCA)
        X_test_PCA.index = X_test.index

        # Define train/test set from the two main principal components.
        data_train = np.array(X_train_PCA.values)
        data_test = np.array(X_test_PCA.values)

        # Calculate the covariance matrix and its inverse, based on data in the training set.
        cov_matrix, inv_cov_matrix = cov_matrix_f(data_train)

        r'''
        We also calculate the mean value for the input variables in the training set,
        as this is used later to calculate the Mahalanobis distance to datapoints in the test set
        '''
        mean_distr = data_train.mean(axis=0)

        r'''
        Using the covariance matrix and its inverse,
        we can calculate the Mahalanobis distance for the training data defining “normal conditions”,
        and find the threshold value to flag datapoints as an anomaly.
        One can then calculate the Mahalanobis distance for the datapoints in the test set
        and compare that with the anomaly threshold.
        '''

        dist_test = MahalanobisDist(inv_cov_matrix, mean_distr, data_test)
        dist_train = MahalanobisDist(inv_cov_matrix, mean_distr, data_train)
        threshold = MD_threshold(dist_train, extreme=True)

        r'''
        Threshold value for flagging an anomaly:

        The square of the Mahalanobis distance to the centroid of the distribution should follow a χ2 distribution
        if the assumption of normal distributed input variables is fulfilled.
        This is also the assumption behind the above calculation of the “threshold value” for flagging an anomaly.
        As this assumption is not necessarily fulfilled in our case,
        it is beneficial to visualize the distribution of the Mahalanobis distance to set a good threshold value for flagging anomalies.

        We start by visualizing the square of the Mahalanobis distance, which should then ideally follow a χ2 distribution.
        '''

        plt.figure()
        sns.distplot(np.square(dist_train),
                     bins=75,
                     kde=False);
        plt.xlim([0.0, 15])
        plt.xlabel('Square of the Mahalanobis distance')
        plt.savefig(r'./data/anom/'+ups+'/MahalanobisDistHistSquared.jpeg',
                    quality=95, optimize=True, progressive=True, format='jpeg')

        # Then visualize the Mahalanobis distance itself
        plt.figure()
        sns.distplot(dist_train,
                     bins=75,
                     kde=True,
                     color='green')
        plt.xlim([0.0, 5])
        plt.xlabel('Mahalanobis dist')
        plt.savefig(r'./data/anom/'+ups+'/MahalanobisDistHist.jpeg',
                    quality=95, optimize=True, progressive=True, format='jpeg')

        r'''
        From the above distributions,
        the calculated threshold value is defined as 3 standard deviations from the center of the distribution.
        '''

        r'''
        We can then save the Mahalanobis distance,
        as well as the threshold value and “anomaly flag” variable for both train and test data in a dataframe.
        '''
        anomaly_train = pd.DataFrame()
        anomaly_train['Mob dist'] = dist_train
        anomaly_train['Thresh'] = threshold
        # If Mob dist above threshold: Flag as anomaly
        anomaly_train['Anomaly'] = anomaly_train['Mob dist'] > anomaly_train['Thresh']
        anomaly_train.index = X_train_PCA.index
        anomaly = pd.DataFrame()
        anomaly['Mob dist'] = dist_test
        anomaly['Thresh'] = threshold
        # If Mob dist above threshold: Flag as anomaly
        anomaly['Anomaly'] = anomaly['Mob dist'] > anomaly['Thresh']
        anomaly.index = X_test_PCA.index

        r'''
        Based on the calculated statistics, any distance above the threshold value will be flagged as an anomaly.
        We can now merge the data in a single dataframe and save it as a .csv file
        '''

        anomaly_alldata = pd.concat([anomaly_train, anomaly], sort=False)
        anomaly_alldata.sort_index(inplace=True)
        anomaly_alldata.to_csv(r'./data/anom/'+ups+'/Anomaly_distance.csv')

        r'''
        Verifying PCA model on test data:
        We can now plot the calculated anomaly metric (Mob dist),
        and check when it crosses the anomaly threshold (note the logarithmic y-axis).
        '''

        ymax = 100*max(anomaly['Mob dist'])
        ymin = 0.001*min(anomaly['Mob dist'])
        fig = anomaly_alldata.plot(logy=True, figsize=(30, 18), ylim=[ymin, ymax], color=['green', 'red'])
        fig.get_figure().savefig(r'./data/anom/'+ups+'/Anomaly_distance.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
        plt.close('all')
        return anomaly_alldata

    def nonlinear_autoencoder_detect(self, X_train, X_test, ups):
        r'''
        The basic idea here is to use an autoencoder neural network to “compress” the sensor readings to a low dimensional representation,
        which captures the correlations and interactions between the various variables.
        (Essentially the same principle as the PCA model, but here we also allow for non-linearities among the input variables).

        We use a 3 layer neural network:
            First layer has 6 nodes, middle layer has 2 nodes, and third layer has 6 nodes.
        We use the mean square error as loss function, and train the model using the SGD optimizer with 0.8 as momentum value.
        '''

        train = X_train.to_numpy()

        corrupted_train = train.copy()

        random_rows_to_corrupt = np.random.randint(0, corrupted_train.shape[0], size=(int(corrupted_train.size * 0.4),))
        random_cols_to_corrupt = np.random.randint(0, corrupted_train.shape[1], size=(int(corrupted_train.size * 0.4),))

        corrupted_train[random_rows_to_corrupt, random_cols_to_corrupt] = 0

        act_func = LeakyReLU(alpha=0.01)

        # Input layer:
        model = Sequential()
        # First hidden layer, connected to input vector X.
        model.add(Dense(6, activation=act_func,
                        kernel_initializer='glorot_uniform',
                        input_shape=(X_train.shape[1],)))

        model.add(Dense(2, activation=act_func,
                        kernel_initializer='glorot_uniform'))

        model.add(Dense(6, activation=act_func,
                        kernel_initializer='glorot_uniform'))

        model.add(Dense(X_train.shape[1],
                        kernel_initializer='glorot_uniform'))

        adam = tf.keras.optimizers.Adam(amsgrad=True)
        model.compile(loss='mae', optimizer=adam)

        # Train model for 200 epochs, batch size of 32:
        NUM_EPOCHS = 350
        BATCH_SIZE = 64

        r'''
        To keep track of the accuracy during training,
        we use 5% of the training data for validation after each epoch (validation_split = 0.05).
        '''

        es = EarlyStopping(monitor='val_loss', mode='min', min_delta=0, verbose=0, patience=25)
        history = model.fit(corrupted_train + np.random.normal(), np.array(X_train),
                          batch_size=BATCH_SIZE,
                          epochs=NUM_EPOCHS,
                          validation_split=0.05,
                          verbose=0,
                          callbacks=[es])

        # Visualize training/validation loss
        plt.plot(history.history['loss'], 'b', label='Training loss')
        plt.plot(history.history['val_loss'], 'r', label='Validation loss')
        plt.legend(loc='upper right')
        plt.xlabel('Epochs')
        plt.ylabel('Loss, [mse]')
        ylim = max(max(history.history['loss']), max(history.history['val_loss']))
        plt.ylim([0, ylim])
        plt.savefig(r'./data/anom/'+ups+'/Autoencoder_loss.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
        plt.close()

        r'''
        Distribution of loss function in the training set:

        By plotting the distribution of the calculated loss in the training set,
        one can use this to identify a suitable threshold value for identifying an anomaly.
        In doing this, one can make sure that this threshold is set above the “noise level”,
        and that any flagged anomalies should be statistically significant above the noise background.
        '''
        X_pred = model.predict(np.array(X_train))
        X_pred = pd.DataFrame(X_pred,
                              columns=X_train.columns)
        X_pred.index = X_train.index

        scored = pd.DataFrame(index=X_train.index)
        scored['Loss_mae'] = np.mean(np.abs(X_pred-X_train), axis=1)
        plt.figure()
        sns.distplot(scored['Loss_mae'],
                     bins=50,
                     kde=True,
                     color='blue')
        plt.xlim([0.0, max(scored['Loss_mae'])])
        plt.savefig(r'./data/anom/'+ups+'/Autoencoder_loss_dist.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
        plt.close()

        r'''
        From the above loss distribution,
        let us try a threshold based on one of the last quantiles for flagging an anomaly.
        We can then calculate the loss in the test set, to check when the output crosses the anomaly threshold.
        '''

        threshold = scored['Loss_mae'].quantile(.95) + 1*(scored['Loss_mae'].quantile(.95) - scored['Loss_mae'].quantile(.05))
        X_pred = model.predict(np.array(X_test))
        X_pred = pd.DataFrame(X_pred,
                              columns=X_test.columns)
        X_pred.index = X_test.index

        scored = pd.DataFrame(index=X_test.index)
        scored['Loss_mae'] = np.mean(np.abs(X_pred-X_test), axis=1)
        scored['Threshold'] = threshold
        scored['Anomaly'] = scored['Loss_mae'] > scored['Threshold']

        fig = scored.plot(logy=True, figsize=(30, 18), ylim=[1e-3, 1e3], color=['blue', 'red'])
        fig.get_figure().savefig(r'./data/anom/'+ups+'/Autoencoder_test_threshold.jpeg',
                                 quality=95, optimize=True, progressive=True, format='jpeg')

        plt.close('all')

        r'''
        We then calculate the same metrics also for the training set, and merge all data in a single dataframe.
        '''
        X_pred = model.predict(np.array(X_test))
        X_pred = pd.DataFrame(X_pred,
                              columns=X_test.columns)
        X_pred.index = X_test.index

        scored_test = pd.DataFrame(index=X_test.index)
        scored_test['Loss_mae'] = np.mean(np.abs(X_pred-X_test), axis=1)
        scored_test['Threshold'] = threshold
        scored_test['Anomaly'] = scored_test['Loss_mae'] > scored_test['Threshold']

        X_pred_train = model.predict(np.array(X_train))
        X_pred_train = pd.DataFrame(X_pred_train,
                                    columns=X_train.columns)
        X_pred_train.index = X_train.index

        scored_train = pd.DataFrame(index=X_train.index)
        scored_train['Loss_mae'] = np.mean(np.abs(X_pred_train-X_train), axis=1)
        scored_train['Threshold'] = threshold
        scored_train['Anomaly'] = scored_train['Loss_mae'] > scored_train['Threshold']

        scored = pd.concat([scored_train, scored_test], sort=False)

        r'''
        Isolation forest prediction based on anomality detected by the autoencoder
        '''
        contamination = (scored.Anomaly == True).sum()/len(scored)

        clf = IsolationForest(behaviour='new', n_estimators=500, max_samples='auto', contamination=contamination, \
                                max_features=1.0, bootstrap=True, n_jobs=-1, random_state=42, verbose=0)

        clf.fit(scored.Loss_mae.to_numpy().reshape(-1, 1))

        pred = clf.predict(scored.Loss_mae.to_numpy().reshape(-1, 1))
        dct = {-1:threshold, 1:np.nan}
        pred = [dct[k] for k in pred]

        r'''
        Results from Autoencoder model:
        Having calculated the loss distribution and the anomaly threshold, we can visualize the model output in the time leading up to the bearing failure:
        '''

        ymax = 100*max(scored['Loss_mae'])
        ymin = 0.001*min(scored['Loss_mae'])
        scored.plot(logy=True, figsize=(30, 18), ylim=[ymin, ymax], color=['green', 'red'])
        isolation_forest = plt.scatter(scored.index, pred, s=150, c='b', marker=".")
        plt.legend((isolation_forest,), ('Isolation Forest',), scatterpoints=1, loc='upper right', ncol=1, fontsize=8)
        plt.savefig(r'./data/anom/'+ups+'/Autoencoder_full_threshold.jpeg', quality=95, optimize=True, progressive=True, format='jpeg')
        plt.close('all')

        scored.sort_index(inplace=True)
        scored.to_csv(r'./data/anom/'+ups+'/Dense_autoencoder_distance.csv')

        return scored

r'''

Summary:
Both modeling approaches give similar results,
where they are able to flag the upcoming bearing malfunction well in advance of the actual failure.
The main difference is essentially how to define a suitable threshold value for flagging anomalies,
to avoid to many false positives during normal operating conditions
'''

Detector().scheduler()
#iaia = Detector().min_cost_transport()

#Detector().nonlinear_autoencoder_detect()
#Detector().nonlinear_LSTM_autoencoder_detect()
#Detector().variational_autoencoder()

max_stars_repo_name: lrnzgiusti/Extreme-Rare-Event-Detector

max_issues_repo_name: lrnzgiusti/Extreme-Rare-Event-Detector

max_forks_repo_name: lrnzgiusti/Extreme-Rare-Event-Detector

reason: import numpy,from numpy,from scipy


text:
#!/usr/bin/env python

#########################################################################
#
#  Code to detect base pairs from correlated chemical probing data
#
#  See README for further details
#  Run with -h flag for arguments
#
#  Lead developer:  Anthony Mustoe
#  Contributors: Nicole Lama, Steve Busan
#
#  This file is licensed under the terms of the MIT license
#
#  Version 1
#  December 2018
#
#########################################################################

import argparse, sys, itertools
import numpy as np

from ringmapper import RINGexperiment
from ReactivityProfile import ReactivityProfile

class PairMapper(object):
    """ class to hold matrix of pairs """

    def __init__(self, ringexp, profile,
                       chi2cut=23.9,
                       primary_reactivity=0.2, primary_zscore=2.0,
                       secondary_reactivity=0.5, secondary_zscore=2.0,
                       maxGU=None, maxNC=0):
        """Initialize PairMapper from ringexp object"""

        self.parent = ringexp
        self.profile = profile

        # first check whether filters are passed
        self.passfilter = self.checkMutationRates(primary_reactivity = primary_reactivity)

        self.complementarycorrs = None
        self.primary = None
        self.secondary = None

        # compute list of complementary correlated pairs
        self.complementarycorrs = self.computeComplementaryCorrs(chi2cut, maxGU=maxGU, maxNC=maxNC)

        self.computePrimaryCorrs(primary_reactivity, primary_zscore)
        self.computeSecondaryCorrs(secondary_reactivity, secondary_zscore)

    def computePrimaryCorrs(self, primary_reactivity, primary_zscore):
        """Compute primary pair signals"""

        corrs = self.getStrongest()

        # now need to filter by reactivity
        corrs = self.filterReactive(corrs, primary_reactivity)

        # filter by zave
        corrs = self.filterZscore(corrs, primary_zscore)

        self.primary = corrs

    def computeSecondaryCorrs(self, secondary_reactivity, secondary_zscore):
        """Compute secondary pair signals"""

        # start from all complementary corrs
        corrs = self.complementarycorrs

        # filter by reactivity
        corrs = self.filterReactive(corrs, secondary_reactivity)

        # filter by zave
        corrs = self.filterZscore(corrs, secondary_zscore)

        # eliminate corrs that are primary
        corrs = list( set(corrs)-set(self.primary) )
        corrs.sort()

        self.secondary = corrs

    def filterZscore(self, corrlist, cutoff):
        """Filter the set of correlations in corrlist by average zscore"""

        filtered = []
        for i,j in corrlist:
            if self.parent.getMeanZ(i,j) >= cutoff:
                filtered.append((i,j))

        return filtered

    def filterReactive(self, corrlist, r_cutoff):
        """Filter the set of correlations in corrlist by reactivity """

        filtered = []
        for c in corrlist:
            r0 = safeReactivityMean(self.profile.normprofile[c[0]:c[0]+self.parent.window])
            r1 = safeReactivityMean(self.profile.normprofile[c[1]:c[1]+self.parent.window])

            if r0<=r_cutoff and r1<=r_cutoff:
                filtered.append(c)

        return filtered

    def getStrongest(self):
        """Return non-conflicting strongest correlations"""

        # get list of strongest correlation at each nt
        corrdict = {}
        for i,j in self.complementarycorrs:
            if i not in corrdict:
                corrdict[i] = (i,j, self.parent.ex_correlations[i,j])
            elif self.parent.ex_correlations[i,j] > corrdict[i][2]:
                corrdict[i] = (i,j, self.parent.ex_correlations[i,j])

            if j not in corrdict:
                corrdict[j] = (i,j, self.parent.ex_correlations[i,j])
            elif self.parent.ex_correlations[i,j] > corrdict[j][2]:
                corrdict[j] = (i,j, self.parent.ex_correlations[i,j])

        # eliminate correlations that are not mutually strongest
        ntkeys = corrdict.keys()
        for k in ntkeys:
            # need to check since we are deleting things
            if k not in corrdict:
                continue

            i,j,c = corrdict[k]
            if k == i and (j not in corrdict or corrdict[j] != corrdict[k]):
                del corrdict[k]
            elif k == j and (i not in corrdict or corrdict[i] != corrdict[k]):
                del corrdict[k]

        # function to test if two correlations are parallel
        def parallel(c1, c2):
            if c1[0]-c2[0] == -1*(c1[1]-c2[1]):
                return True
            return False

        # build final list while filtering out conflicting windows
        output = []
        for k in corrdict:

            # only consider/add each correlation once
            if k != corrdict[k][0]:
                continue

            # scan up and downstream of each correlation leg to see if there is
            # a stronger conflicting correlation
            keep = True
            for i in range(1-self.parent.window,0)+range(1,self.parent.window):

                if k+i in corrdict and not parallel(corrdict[k], corrdict[k+i]) \
                                   and corrdict[k+i][2] > corrdict[k][2]:
                    keep = False
                    break

                # repeat for other side of the correlation
                j = corrdict[k][1]
                if j+i in corrdict and not parallel(corrdict[k], corrdict[j+i]) \
                                   and corrdict[j+i][2] > corrdict[k][2]:
                    keep = False
                    break

            if keep:
                output.append(corrdict[k][:2])

        output.sort()

        return output

    def computeComplementaryCorrs(self, chi2cut, maxGU=None, maxNC=0):
        """Return list of complementary correlations"""

        def isComplement(seq1, seq2, maxGU, maxNC):
            counts = [0,0,0]

            # Not considered complementary if different lengths
            if len(seq1) != len(seq2):
                return False

            for i in range(len(seq1)):
                pair = seq1[i] + seq2[-(i+1)]

                if pair in ('AU', 'UA', 'TA','AT', 'GC', 'CG'):
                    counts[0]+=1
                elif pair in ('GU', 'UG', 'GT', 'TG'):
                    counts[1]+=1
                else:
                    counts[2]+=1

            if counts[2] > maxNC or counts[1] > maxGU:
                return False

            return True

        # default is all GU pairs are valid
        if maxGU is None:
            maxGU = self.parent.window

        # get all positive correlations
        corrs = self.parent.significantCorrelations('ex', chi2cut, sign=1)

        compcorrs = []

        for i,j in corrs:

            if isComplement(self.parent.sequence[i:i+self.parent.window],
                    self.parent.sequence[j:j+self.parent.window], maxGU=maxGU, maxNC=maxNC):
                compcorrs.append((i,j))

        return compcorrs

    def writePairs(self, outname):

        def writeline(OUT, i, j, code):
            #function to write a single line of output
            OUT.write("{0}\t{1}\t".format(i+1, j+1))
            OUT.write("{0:.2f}\t".format(self.parent.ex_correlations[i,j]))
            OUT.write("{0}\t".format(code))
            OUT.write("{0:.2f}\t".format(self.parent.getMeanZ(i,j)))
            OUT.write("{0:.2f}\t{1:.2f}\t".format(self.parent.ex_zscores[i,j],
                                                  self.parent.ex_zscores[j,i]))
            OUT.write("{0}\t{1}\t".format(self.parent.ex_readarr[i,j],
                                          self.parent.ex_comutarr[i,j]))
            OUT.write("{0}\t{1}\n".format(self.parent.bg_readarr[i,j],
                                        self.parent.bg_comutarr[i,j]))

        # determine complementary correlations that aren't primary or secondary
        remainder = set(self.complementarycorrs)-set(self.primary)-set(self.secondary)
        remainder = list(remainder)
        remainder.sort()

        with open(outname,'w') as outf:

            outf.write("{0}\tWindow={1}\tMetric=PAIRMAPPER\t".format(len(self.parent.sequence), self.parent.window))
            outf.write("\n".format())

            outf.write("i\tj\tSig.\tClass\tZij\tZi\tZj\tMod_Depth\tMod_Comuts\tUnt_Depth\tUnt_Comuts\n")

            for i,j in self.primary:
                writeline(outf, i, j, 1)

            for i,j in self.secondary:
                writeline(outf, i, j, 2)

            for i,j in remainder:
                writeline(outf, i, j, 0)

    def writePairBonusFile(self, filepath, chi2cut=23.9, maxNC=1, scale=0.5, intercept=0, fileformat=1):
        """Write matrix of pairing bonuses for use in RNAstructure folding"""

        seqlen = len(self.parent.sequence)
        pairmat = np.zeros((seqlen, seqlen))

        corrs = self.computeComplementaryCorrs(chi2cut, maxNC=maxNC)

        for i,j in corrs:

            meanz = self.parent.getMeanZ(i,j)
            if meanz > 1:
                energy = -0.616*scale*np.log(meanz) - intercept

                # add to each base pair
                for w in range(self.parent.window):
                    pairmat[i+w, j+self.parent.window-1-w] += energy

        # make symmetric
        pairmat += pairmat.transpose()

        if fileformat==0:
            np.savetxt(filepath, pairmat, fmt='%.4f')

        else:
            with open(filepath,'w') as out:
                out.write('; i j bonus\n')
                for i in xrange(pairmat.shape[0]):
                    for j in xrange(i+1, pairmat.shape[0]):
                        if pairmat[i,j] != 0:
                            out.write('{0} {1} {2:.7f}\n'.format(i+1, j+1, pairmat[i,j]))

    def checkMutationRates(self, depthcut=10000, primary_reactivity = 0.2):
        """Check comutation rates and read depths to make sure they are high enough to accurately measure pairs"""

        seqlen = self.parent.ex_readarr.shape[0]

        # initiate the masked array
        comuts = np.ma.masked_array(self.parent.ex_comutarr, copy=True, dtype=np.float64)

        # mask out the lower diagonal (which is void of info)
        comuts[np.tril_indices(seqlen)] = np.ma.masked

        # mask out nts with no mutation signal in the modifed
        nts = self.parent.getUnreactiveNts(5e-4, prefix='ex')
        for i in nts:
            comuts[i,:] = np.ma.masked
            comuts[:,i] = np.ma.masked

        # mask out high background positions
        nts = self.parent.getReactiveNts(self.parent.highbgrate, prefix='bg')
        for i in nts:
            comuts[i,:] = np.ma.masked
            comuts[:,i] = np.ma.masked

        # compute the comut rates
        for i in range(seqlen):

            # mask out near diagonals
            for j in range(i+1, min(i+self.parent.corrbuffer, seqlen)):
                comuts[i,j] = np.ma.masked

            # now go through remaining (valid) nts
            for j in range(i+self.parent.corrbuffer, seqlen):

                # convert to rate
                if not comuts.mask[i,j] and self.parent.ex_readarr[i,j] > depthcut:
                    comuts[i,j] /= self.parent.ex_readarr[i,j]

                    # do background subtraction
                    if self.parent.bg_readarr[i,j] > depthcut:
                        comuts[i,j] -= float(self.parent.bg_comutarr[i,j])/self.parent.bg_readarr[i,j]
                    else:
                        comuts[i,j] = np.ma.masked

                else:
                    comuts[i,j] = np.ma.masked

        # now compute the median over valid (non-masked) values
        median_comut = np.ma.median(comuts)
        median_count = comuts.count()

        # now compute median comut rate of unreactive positions by masking out reactive nts
        for i in range(comuts.shape[0]):

            react = safeReactivityMean( self.profile.normprofile[i:i+self.parent.window] )

            if react is np.nan or react > primary_reactivity:
                comuts[i,:] = np.ma.masked
                comuts[:,i] = np.ma.masked

        median_unreact_comut = np.ma.median(comuts)
        median_unreact_count = comuts.count()

        allpassed = True

        # mutation rate checks
        if median_comut < 5e-4 or median_unreact_comut < 1e-4:
            sys.stdout.write('******************************\n')
            sys.stdout.write('WARNING: Low comutation rates!\n')
            sys.stdout.write('\tMedian for all nts = {0:.2e}  ({1} pairs)\n'.format(median_comut, median_count))
            sys.stdout.write('\tMedian for unreactive nts = {0:.2e}  ({1} pairs)\n'.format(median_unreact_comut, median_unreact_count))
            sys.stdout.write('PAIR-MaP data likely untrustworthy\n')

            allpassed = False

        else:
            sys.stdout.write("Reactivity rate quality checks passed\n")
            sys.stdout.write('\tMedian for all nts = {0:.2e}  ({1} pairs)\n'.format(median_comut, median_count))
            sys.stdout.write('\tMedian for unreactive nts = {0:.2e}  ({1} pairs)\n'.format(median_unreact_comut, median_unreact_count))

        # read depth checks
        median_depth = np.median(np.diag(self.parent.ex_readarr))

        if median_depth<400000:
            sys.stdout.write('******************************\n')
            sys.stdout.write("WARNING: Low read depths!\n")
            sys.stdout.write("Median read depth = {}!\n".format(median_depth))
            allpassed = False

        return allpassed

    def plot(self, output, minz=2, maxz=6, ct=None):
        """Plot reactivity data and primary+secondary pairmap correlations"""

        from pmanalysis import PairMap
        from arcPlot import ArcPlot

        plot = ArcPlot()

        plot.reactprofile = self.profile.normprofile
        plot.reactprofileType = 'DMS'
        plot.seq = ''.join(self.profile.sequence)

        pmobj = PairMap()
        pmobj.primary = [(i+1,j+1,0, self.parent.getMeanZ(i,j)) for i,j in self.primary]
        pmobj.secondary = [(i+1,j+1,0, self.parent.getMeanZ(i,j)) for i,j in self.secondary]
        pmobj.window = self.parent.window

        plot.addPairMap( pmobj, panel=-1)

        msg=None
        if ct is not None:
            plot.addCT(ct)

            p,s = pmobj.ppvsens_duplex(ct, ptype=1, profile=self.profile.normprofile)
            msg = "PPV={0:.0f}  Sens={1:.0f}".format(p*100, s*100)

        plot.writePlot(output, msg=msg)

#########################################################################################

def safeReactivityMean(rarray):
    """Compute the mean reactivity of the provided array, setting negative values to 0"""

    c = 0
    n = 0

    for r in rarray:
        if r == r and r>-10: # make sure defined
            n += 1
            c += max(r, 0) # don't let negatives skew

    if n==0:
        return np.nan
    else:
        return c/n

def parseArguments():

    parser = argparse.ArgumentParser(description = "Detect base pairs from correlated DMS data")

    optional = parser._action_groups.pop()
    required = parser.add_argument_group('required arguments')

    required.add_argument('--modified_parsed', help='Path to Modified parsed.mut')
    required.add_argument('--untreated_parsed', help='Path to Untreated parsed.mut')
    required.add_argument('--profile', help='Path to profile.txt')
    required.add_argument('--out', help='output prefix')

    ######################################################
    # optional arguments

    optional.add_argument('--undersample', type=int, default=-1, help="""Randomly undersample specified number of reads
                                           from modified file (default=-1 [disabled]).""")

    optional.add_argument('--ignorents', help="""A list of comma-separated (1-indexed) nts to ignore (e.g. 10,12,35)""")

    optional.add_argument('--mincoverage', type=float, default=0, help="""Quality filter reads by requiring a minimum
                                           number of positional matches to reference sequence.
                                           This parameter can be set as a fraction of the total molecule length
                                           (i.e. if 0.8 is passed, reads will be required to have at least 0.8*len(fasta)
                                           valid matches) Alternatively, an integer value can be passed
                                           (i.e. if 150 is passed, reads will be required to have at least 150 valid matches).
                                           By default, this filter is disabled.""")

    optional.add_argument('--chisq_cut', type=float, default=23.9, help="Set chisq cutoff (default = 23.9)")

    optional.add_argument('--highbg_rate', type=float, default = 0.02, help="""Ignore nt position with bg reactivity above
                                           this value. For PAIR-MaP analysis, windows with bg 3X this rate will be ignored""")

    optional.add_argument('--highbg_corr', type=float, default = 10.83, help="""Ignore nt pairs correlated in the bg
                                           sample, with correlation determined via this significance value
                                           (default=10.83 --> P=1e-3)""")

    optional.add_argument('--mincorrdistance', type=int, default=6, help="""Minimum distance allowed between correlations
                                               (default = 6)""")

    optional.add_argument('--mindepth', type=int, default=10000, help="""Minimum pairwise read depth allowed for calculating
                                        correlations (default = 10000)""")

    optional.add_argument('--mincount', type=int, default=50, help="""Minimum required count in contigency table
                                      (default = 50). Nt pairs with fewer than this number of comutations are ignored""")

    optional.add_argument('--primary_reactivity', type=float, default=0.2, help="""Reactivity cutoff for primary correlations
                                      (default = 0.2)""")
    optional.add_argument('--primary_zscore', type=float, default=2.0, help="""Zscore cutoff for primary correlations
                                      (default = 2.0)""")
    optional.add_argument('--secondary_reactivity', type=float, default=0.5, help="""Reactivity cutoff for secondary correlations
                                      (default = 0.5)""")
    optional.add_argument('--secondary_zscore', type=float, default=2.0, help="""Zscore cutoff for secondary correlations
                                      (default = 2.0)""")

    optional.add_argument('--override_qualcheck', action='store_true', help="""Override quality checks and perform PAIR-MaP analysis despite poor data quality""")

    parser._action_groups.append(optional)

    args = parser.parse_args()

    # check that required args are defined
    if not args.modified_parsed or not args.untreated_parsed \
            or not args.profile or not args.out:

        sys.stderr.write("! Required arguments not provided !\n\n")
        parser.print_help(sys.stderr)
        sys.exit(1)

    # parse ignorents argument
    if args.ignorents:
        spl = args.ignorents.split(',')
        ig = []
        try:
            for x in spl:
                if ':' in x:
                    xspl = x.split(':')
                    ig.extend(range(int(xspl[0])-1, int(xspl[1])))
                else:
                    ig.append(int(x)-1)  # correct for 1-indexing of input
        except ValueError:
            raise ValueError('Invalid ignorents option: {0}'.format(args.ignorents))

        args.ignorents = ig
    else:
        args.ignorents = []

    return args

if __name__ == '__main__':

    verbal = True

    args = parseArguments()

    # read in the reactivity profile information
    profile = ReactivityProfile(args.profile, bg=args.highbg_rate, ignorents=args.ignorents)
    profile.normalize(DMS=True)
    # write normed reactivities to file
    profile.writeReactivity('{}.dms'.format(args.out))

    # initialize ringmapper and data structures
    ringexp = RINGexperiment(corrtype = 'apc', arraysize=len(profile.sequence))
    ringexp.sequence = profile.sequence

    ringexp.initDataMatrices('ex', args.modified_parsed, window=3,
                             mincoverage=args.mincoverage, undersample = args.undersample,
                             verbal = verbal)

    ringexp.initDataMatrices('bg', args.untreated_parsed, window=3,
                             mincoverage=args.mincoverage, undersample = args.undersample,
                             verbal = verbal)

    # compute the correlation matrices
    ringexp.computeCorrelationMatrix(corrbuffer = args.mincorrdistance,
                                     mindepth = args.mindepth,
                                     mincount = args.mincount,
                                     ignorents = args.ignorents,
                                     highbgrate = args.highbg_rate,
                                     highbgcorr = args.highbg_corr,
                                     verbal = verbal)

    # at this point, all the major calculations have been done

    # write all correlations to file
    ringexp.writeCorrelations('{0}-allcorrs.txt'.format(args.out), chi2cut=args.chisq_cut)

    # filter out pairs
    pairs = PairMapper(ringexp, profile, chi2cut=args.chisq_cut,
                       primary_reactivity=args.primary_reactivity,
                       primary_zscore=args.primary_zscore,
                       secondary_reactivity=args.secondary_reactivity,
                       secondary_zscore=args.secondary_zscore)

    if pairs.passfilter or args.override_qualcheck:

        # write out pairmap data
        pairs.writePairs('{}-pairmap.txt'.format(args.out))

        # write out folding restraint matrix
        pairs.writePairBonusFile('{}.bp'.format(args.out), chi2cut=args.chisq_cut, maxNC=1)

    # make plot
    #try:
    #    pairs.plot('{}-pairmapper.pdf'.format(args.out))
    #except:
    #    print("Was unable to generate plot")

# vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4

ext: py

max_stars_repo_name: chaseaw/RingMapper

max_issues_repo_name: chaseaw/RingMapper

max_forks_repo_name: chaseaw/RingMapper

reason: import numpy


text:
"""
Utils to replicate Twins experiments (Appendix E.2)
"""
# Author: Alicia Curth
import csv
import os
from pathlib import Path

from sklearn import clone
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score, average_precision_score
from sklearn.preprocessing import label_binarize
import numpy as np

from catenets.datasets import load
from catenets.models.jax.base import check_shape_1d_data
from catenets.experiment_utils.base import eval_root_mse
from catenets.models.jax import TARNet, TNet, OffsetNet, FlexTENet, DRNet
from catenets.models.jax import TARNET_NAME, T_NAME, OFFSET_NAME, FLEXTE_NAME, DRNET_NAME

RESULT_DIR = Path("results/experiments_inductive_bias/twins")
SEP = "_"

PARAMS_DEPTH = {'n_layers_r': 1, 'n_layers_out': 1}
PARAMS_DEPTH_2 = {'n_layers_r': 1, 'n_layers_out': 1, 'n_layers_r_t': 1, 'n_layers_out_t': 1}
PENALTY_DIFF = 0.01
PENALTY_ORTHOGONAL = 0.1

ALL_MODELS = {
    T_NAME: TNet(**PARAMS_DEPTH),
    T_NAME
    + "_reg": TNet(train_separate=False, penalty_diff=PENALTY_DIFF, **PARAMS_DEPTH),
    TARNET_NAME: TARNet(**PARAMS_DEPTH),
    TARNET_NAME
    + "_reg": TARNet(
        reg_diff=True, penalty_diff=PENALTY_DIFF, same_init=True, **PARAMS_DEPTH
    ),
    OFFSET_NAME: OffsetNet(penalty_l2_p=PENALTY_DIFF, **PARAMS_DEPTH),
    FLEXTE_NAME: FlexTENet(
        penalty_orthogonal=PENALTY_ORTHOGONAL, penalty_l2_p=PENALTY_DIFF, **PARAMS_DEPTH
    ),
    FLEXTE_NAME + "_noortho_reg_same": FlexTENet(penalty_orthogonal=0, **PARAMS_DEPTH),
    DRNET_NAME: DRNet(**PARAMS_DEPTH_2),
    DRNET_NAME + "_TAR": DRNet(first_stage_strategy="Tar", **PARAMS_DEPTH_2),
}

def do_twins_experiment_loop(n_train_loop=[500, 1000, 2000, 5000, None],
                             prop_loop=[0.1, 0.25, 0.5, 0.75, 0.9],
                             n_exp: int = 10, file_name: str = 'twins', models: dict = None,
                             test_size=0.5):
    for n in n_train_loop:
        for prop in prop_loop:
            print('Running twins experiment for {} training samples with {} treated.'.format(n,
                                                                                prop))
            do_twins_experiments(n_exp=n_exp, file_name=file_name, models=models, subset_train=n,
                                 prop_treated=prop, test_size=test_size)

def do_twins_experiments(n_exp: int = 10, file_name: str = 'twins',
                         models: dict = None, subset_train: int = None, prop_treated=0.5,
                         test_size=0.5):
    if models is None:
        models = ALL_MODELS

    # get file to write in
    if not os.path.isdir(RESULT_DIR):
        os.makedirs(RESULT_DIR)

    out_file = open(RESULT_DIR / (file_name + SEP + str(prop_treated) + SEP + str(subset_train) +
                    '.csv'), 'w', buffering=1)
    writer = csv.writer(out_file)
    header = [name + '_cate' for name in models.keys()] + \
             [name + '_auc_ite' for name in models.keys() if 'R' not in name and 'X' not in
              name] + \
             [name + '_auc_mu0' for name in models.keys() if 'R' not in name and 'X' not in
              name] + \
             [name + '_auc_mu1' for name in models.keys() if 'R' not in name and 'X' not in
              name] + [name + '_ap_mu0' for name in models.keys() if 'R' not in name and 'X'
                       not in name] + \
             [name + '_ap_mu1' for name in models.keys() if 'R' not in name and 'X' not in name]

    writer.writerow(header)

    for i_exp in range(n_exp):
        pehe_out = []
        auc_ite = []
        auc_mu0 = []
        auc_mu1 = []
        ap_mu0 = []
        ap_mu1 = []

        # get data
        x, w, y, pos, _, _ = load('twins', seed=i_exp, treat_prop=prop_treated, train_ratio=1)

        # split data
        X, X_t, y, y_t, w, w_t, y0_in, y0_out, y1_in, y1_out = split_data(x, y, w, pos,
                                                                          random_state=i_exp,
                                                                          subset_train=subset_train,
                                                                          test_size=test_size)

        ite_out = y1_out - y0_out

        ite_out_encoded = label_binarize(ite_out, [-1, 0, 1])

        n_test = X_t.shape[0]

        # split data
        for model_name, estimator in models.items():
            print("Experiment {} with {}".format(i_exp, model_name))
            estimator_temp = clone(estimator)
            estimator_temp.set_params(**{'binary_y': True})

            # fit estimator
            estimator_temp.fit(X=X, y=y, w=w)

            if 'DR' not in model_name and 'R' not in model_name and 'X' not in model_name:
                cate_pred_out, mu0_pred, mu1_pred = estimator_temp.predict(X_t, return_po=True)

                # create probabilities for each possible level of ITE
                probs = np.zeros((n_test, 3))
                probs[:, 0] = (mu0_pred * (1 - mu1_pred)).reshape((-1,))  # P(Y1-Y0=-1)
                probs[:, 1] = ((mu0_pred * mu1_pred) + ((1 - mu0_pred) * (1 - mu1_pred))).reshape(
                    (-1,))  # P(Y1-Y0=0)
                probs[:, 2] = (mu1_pred * (1 - mu0_pred)).reshape((-1,))  # P(Y1-Y0=1)
                auc_ite.append(roc_auc_score(ite_out_encoded, probs))

                # evaluate performance on potential outcomes
                auc_mu0.append(eval_roc_auc(y0_out, mu0_pred))
                auc_mu1.append(eval_roc_auc(y1_out, mu1_pred))
                ap_mu0.append(eval_ap(y0_out, mu0_pred))
                ap_mu1.append(eval_ap(y1_out, mu1_pred))
            else:
                cate_pred_out = estimator_temp.predict(X_t)

            pehe_out.append(eval_root_mse(cate_pred_out, ite_out))

        writer.writerow(pehe_out + auc_ite + auc_mu0 + auc_mu1 + ap_mu0 + ap_mu1)

    out_file.close()

# utils -------
def split_data(X, y, w, pos, test_size=0.5, random_state=42,
               subset_train: int = None):
    X, X_t, y, y_t, w, w_t, y0_in, \
        y0_out, y1_in, y1_out = train_test_split(X, y, w, pos[:, 0], pos[:, 1],
                                                 test_size=test_size,
                                                 random_state=random_state)
    if subset_train is not None:
        X, y, w, y0_in, y1_in = X[:subset_train, :], y[:subset_train], w[:subset_train], \
                                y0_in[:subset_train], y1_in[:subset_train]

    return X, X_t, y, y_t, w, w_t, y0_in, y0_out, y1_in, y1_out

def eval_roc_auc(targets, preds):
    preds = check_shape_1d_data(preds)
    targets = check_shape_1d_data(targets)
    return roc_auc_score(targets, preds)

def eval_ap(targets, preds):
    preds = check_shape_1d_data(preds)
    targets = check_shape_1d_data(targets)
    return average_precision_score(targets, preds)

max_stars_repo_name: AliciaCurth/CATENets

max_issues_repo_name: vanderschaarlab/CATENets

max_forks_repo_name: vanderschaarlab/CATENets


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------
from pathlib import Path
from typing import Any, Dict, List, Tuple

import numpy as np
import torch
from pytorch_lightning import LightningDataModule, LightningModule
from pytorch_lightning.metrics import MeanAbsoluteError
from torch.optim import Adam, Optimizer
from torch.optim.lr_scheduler import StepLR, _LRScheduler
from torch.utils.data import DataLoader, Dataset

from InnerEye.Common import fixed_paths
from InnerEye.ML.lightning_container import LightningContainer

class HelloDataset(Dataset):
    """
    A simple 1dim regression task, read from a data file stored in the test data folder.
    """

    # Creating the data file:
    # import numpy as np
    # import torch
    #
    # N = 100
    # x = torch.rand((N, 1)) * 10
    # y = 0.2 * x + 0.1 * torch.randn(x.size())
    # xy = torch.cat((x, y), dim=1)
    # np.savetxt("InnerEye/ML/configs/other/hellocontainer.csv", xy.numpy(), delimiter=",")
    def __init__(self, root_folder: Path, start_index: int, end_index: int) -> None:
        """
        Creates the 1-dim regression dataset.
        :param root_folder: The folder in which the data file lives ("hellocontainer.csv")
        :param start_index: The first row to read.
        :param end_index: The last row to read (exclusive)
        """
        super().__init__()
        raw_data = np.loadtxt(str(root_folder / "hellocontainer.csv"), delimiter=",")[start_index:end_index]
        self.data = torch.tensor(raw_data, dtype=torch.float)

    def __len__(self) -> int:
        return self.data.shape[0]

    def __getitem__(self, item: int) -> Dict[str, torch.Tensor]:
        return {'x': self.data[item][0:1], 'y': self.data[item][1:2]}

class HelloDataModule(LightningDataModule):
    """
    A data module that gives the training, validation and test data for a simple 1-dim regression task.
    """

    def __init__(self, root_folder: Path) -> None:
        super().__init__()
        self.train = HelloDataset(root_folder, start_index=0, end_index=50)
        self.val = HelloDataset(root_folder, start_index=50, end_index=70)
        self.test = HelloDataset(root_folder, start_index=70, end_index=100)

    def train_dataloader(self, *args: Any, **kwargs: Any) -> DataLoader:
        return DataLoader(self.train, batch_size=5)

    def val_dataloader(self, *args: Any, **kwargs: Any) -> DataLoader:
        return DataLoader(self.val, batch_size=5)

    def test_dataloader(self, *args: Any, **kwargs: Any) -> DataLoader:
        return DataLoader(self.test, batch_size=5)

class HelloRegression(LightningModule):
    """
    A simple 1-dim regression model.
    """

    def __init__(self) -> None:
        super().__init__()
        self.model = torch.nn.Linear(in_features=1, out_features=1, bias=True)
        self.test_mse: List[torch.Tensor] = []
        self.test_mae = MeanAbsoluteError()

    def forward(self, x: torch.Tensor) -> torch.Tensor:  # type: ignore
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        It runs a forward pass of a tensor through the model.
        :param x: The input tensor(s)
        :return: The model output.
        """
        return self.model(x)

    def training_step(self, batch: Dict[str, torch.Tensor], *args: Any, **kwargs: Any) -> torch.Tensor:  # type: ignore
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        It consumes a minibatch of training data (coming out of the data loader), does forward propagation, and
        computes the loss.
        :param batch: The batch of training data
        :return: The loss value with a computation graph attached.
        """
        loss = self.shared_step(batch)
        self.log("loss", loss, on_epoch=True, on_step=False)
        return loss

    def validation_step(self, batch: Dict[str, torch.Tensor], *args: Any,  # type: ignore
                        **kwargs: Any) -> torch.Tensor:
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        It consumes a minibatch of validation data (coming out of the data loader), does forward propagation, and
        computes the loss.
        :param batch: The batch of validation data
        :return: The loss value on the validation data.
        """
        loss = self.shared_step(batch)
        self.log("val_loss", loss, on_epoch=True, on_step=False)
        return loss

    def shared_step(self, batch: Dict[str, torch.Tensor]) -> torch.Tensor:
        """
        This is a convenience method to reduce code duplication, because training, validation, and test step share
        large amounts of code.
        :param batch: The batch of data to process, with input data and targets.
        :return: The MSE loss that the model achieved on this batch.
        """
        input = batch["x"]
        target = batch["y"]
        prediction = self.forward(input)
        return torch.nn.functional.mse_loss(prediction, target)

    def configure_optimizers(self) -> Tuple[List[Optimizer], List[_LRScheduler]]:
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        It returns the PyTorch optimizer(s) and learning rate scheduler(s) that should be used for training.
=        """
        optimizer = Adam(self.parameters(), lr=1e-1)
        scheduler = StepLR(optimizer, step_size=20, gamma=0.5)
        return [optimizer], [scheduler]

    def on_test_epoch_start(self) -> None:
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        In this method, you can prepare data structures that need to be in place before evaluating the model on the
        test set (that is done in the test_step).
        """
        self.test_mse = []
        self.test_mae.reset()

    def test_step(self, batch: Dict[str, torch.Tensor], batch_idx: int) -> torch.Tensor:  # type: ignore
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        It evaluates the model in "inference mode" on data coming from the test set. It could, for example,
        also write each model prediction to disk.
        :param batch: The batch of test data.
        :param batch_idx: The index (0, 1, ...) of the batch when the data loader is enumerated.
        :return: The loss on the test data.
        """
        input = batch["x"]
        target = batch["y"]
        prediction = self.forward(input)
        # This illustrates two ways of computing metrics: Using standard torch
        loss = torch.nn.functional.mse_loss(prediction, target)
        self.test_mse.append(loss)
        # Metrics computed using PyTorch Lightning objects. Note that these will, by default, attempt
        # to synchronize across GPUs.
        self.test_mae.update(preds=prediction, target=target)
        return loss

    def on_test_epoch_end(self) -> None:
        """
        This method is part of the standard PyTorch Lightning interface. For an introduction, please see
        https://pytorch-lightning.readthedocs.io/en/stable/starter/converting.html
        In this method, you can finish off anything to do with evaluating the model on the test set,
        for example writing aggregate metrics to disk.
        """
        average_mse = torch.mean(torch.stack(self.test_mse))
        Path("test_mse.txt").write_text(str(average_mse.item()))
        Path("test_mae.txt").write_text(str(self.test_mae.compute()))

class HelloContainer(LightningContainer):
    """
    An example for using the InnerEye functionality to "bring your own lightning model". This container has methods
    to generate the actual Lightning model, and read out the datamodule that will be used for training.
    The number of training epochs is controlled at container level.
    You can train this model by running `python InnerEye/ML/runner.py --model=HelloContainer` on the local box,
    or via `python InnerEye/ML/runner.py --model=HelloContainer --azureml=True` in AzureML
    """

    def __init__(self) -> None:
        super().__init__()
        self.local_dataset = fixed_paths.repository_root_directory() / "InnerEye" / "ML" / "configs" / "other"
        self.num_epochs = 20

    # This method must be overridden by any subclass of LightningContainer. It returns the model that you wish to
    # train, as a LightningModule
    def create_model(self) -> LightningModule:
        return HelloRegression()

    # This method must be overridden by any subclass of LightningContainer. It returns a data module, which
    # in turn contains 3 data loaders for training, validation, and test set.
    def get_data_module(self) -> LightningDataModule:
        assert self.local_dataset is not None
        return HelloDataModule(root_folder=self.local_dataset)  # type: ignore

    # This is an optional override: This report creation method can read out any files that were written during
    # training, and cook them into a nice looking report. Here, the report is a simple text file.
    def create_report(self) -> None:
        # This just prints out the test MSE, but you could also generate a Jupyter notebook here, for example.
        test_mse = Path("test_mse.txt").read_text().strip()
        test_mae = Path("test_mae.txt").read_text().strip()
        report = f"Performance on test set: MSE = {test_mse}, MAE = {test_mae}"
        print(report)
        Path("report.txt").write_text(report)

max_stars_repo_name: albernsurya/InnerEye-DeepLearning

max_issues_repo_name: albernsrya/InnerEye-DeepLearning

max_forks_repo_name: albernsrya/InnerEye-DeepLearning


text:
'''
@author:     Zhengguang Zhao
@copyright:  Copyright 2016-2019, Zhengguang Zhao.
@license:    MIT
@contact:    zg.zhao@outlook.com

'''

import os
import sys

from sklearn import svm, preprocessing
from sklearn.metrics import confusion_matrix
from sklearn.preprocessing import normalize, scale

from scipy import stats

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.axes_grid1 import make_axes_locatable

import seaborn as sns

import numpy as np

import pandas as pd
from pandas import set_option

import pickle

from utils.io import load_data, load_csv, display_adj_cm, display_cm, read_data, feature_normalize, wave_norm, create_directory
from utils.featextractor import FeatureExtractor
from utils.plot import visualize_ml_result, plot_coefficients, crossplot_features, crossplot_dual_features,\
     crossplot_pca, heatplot_pca, plot_correlations, compare_classifiers, plot_ROC, calc_auc

def create_classifier(classifier_name):

    if classifier_name=='SVM':

        from classifiers import svm
        return svm.ClassifierSVM()

    elif classifier_name=='KNN':

        from classifiers import knn
        return knn.ClassifierKNN()

    elif classifier_name=='DT':

        from classifiers import dtree
        return dtree.ClassifierDecisionTree()

    elif classifier_name=='LogReg':

        from classifiers import logreg
        return logreg.ClassifierLogisticRegression()

    elif classifier_name=='LDA':

        from classifiers import lda
        return lda.ClassifierLDA()

    elif classifier_name=='QDA':

        from classifiers import qda
        return qda.ClassifierQDA()

def main():

    ############################################### main ###########################################

    # change this directory for your machine

    # it should contain the archive folder containing both univariate and multivariate archives

    root_dir = 'F:\\datafolder\\dl4ms_data\\dataset'

    archive_name = 'UTS' #

    dataset_name = 'TX_DEMO_P_256' #'MP_NOISE_PSN_ZZ_256'#MP_PSN_ZZ_256'#'MULTIWELL_B_P_512' #'MP_NOISE_P_128' #  #'TX_PSN_128'# 'TX_DEMO_P_256'  ##

    segment_size = int(dataset_name.split('_')[-1])

    classifier_name= 'SVM'#

    model_dir = os.path.join(root_dir, 'models', classifier_name, archive_name, dataset_name)

    output_directory = os.path.join(root_dir, 'results', classifier_name, archive_name, dataset_name)

    create_directory(output_directory)
    create_directory(model_dir)

    print('\nInfo: ',archive_name, dataset_name, segment_size, classifier_name, '\n')

    file_name = os.path.join(root_dir,'archives', archive_name, dataset_name, dataset_name +'.csv')

    column_name = ['ID', 'FileName', 'Class'] + list(range(segment_size))
    hd = 0 # .csv has header
    datasets_df = load_csv(file_name, hd, column_name)

    if dataset_name.split('_')[-2] == 'PS' or dataset_name.split('_')[-2] == 'P':
        class_labels = ['Event', 'Noise']
    else:
        class_labels = ['P-wave Event', 'S-wave Event','Noise']

    ## Extract features
    file_name = os.path.join(root_dir,'archives', archive_name, dataset_name, dataset_name +'_features.csv')
    if not os.path.exists(file_name):
        extractor = FeatureExtractor()
        extractor.set_dataset(datasets_df)
        extractor.set_class_labels(class_labels)

        fs = 500 # unit is Hz
        window_length = segment_size  # a wavelength is usually 30 samples, we choose 2*wavelength
        overlap_length = int(window_length/2)
        signal_length = segment_size

        extractor.extract_features(fs, window_length, overlap_length, signal_length)
        extractor.save_features(file_name)
        training_data = extractor.feature_data

    else:
        training_data = pd.read_csv(file_name, header= 0, index_col= False)

    ## Features Vector Conditioning
    print(training_data.describe())
    data = training_data.copy()
    for columname in data.columns:
        if data[columname].count() != len(data):
            loc = data[columname][data[columname].isnull().values==True].index.tolist()
            print('Column Name："{}", Row #{} has null value.'.format(columname,loc))

    #blind = training_data[training_data['FileName'] == '41532_140407_224100_1000']
    #training_data = training_data[training_data['FileName'] != '41532_140407_224100_1000']
    training_data['FileName'] = training_data['FileName'].astype('category')
    print(training_data['FileName'].unique())

    # class_colors = ['#196F3D', '#F5B041']
    # class_labels = ['Event', 'Noise']
    # #class_color_map is a dictionary that maps class labels
    # #to their respective colors
    # class_color_map = {}
    # for ind, label in enumerate(class_labels):
    #     class_color_map[label] = class_colors[ind]

    # training_data.dropna(axis=0,how='any') #drop all rows that have any NaN values
    #training_data.fillna(0)

    #count the number of unique entries for each class, sort them by
    #class number (instead of by number of entries)
    class_counts = training_data['Class'].value_counts().sort_index()
    #use class labels to index each count
    class_counts.index = class_labels
    # class_counts.plot(kind='bar',color=class_colors,
    #                 title='Distribution of Training Data by Class')
    print(class_counts)

    # ## save plot display settings to change back to when done plotting with seaborn
    # inline_rc = dict(mpl.rcParams)

    # sns.set()
    # sns.pairplot(training_data.drop(['FileName','Class'],axis=1), hue='ClassLabels', palette=class_color_map, hue_order=list(reversed(class_labels)))

    # #switch back to default matplotlib plot style
    # mpl.rcParams.update(inline_rc)

    ## Visualize dataset and gain insight
    # crossplot_dual_features(data, ['Peak Frequency', 'Shannon Entropy'])
    #crossplot_features(training_data, ['Event', 'Noise'])

    feature_names = training_data.columns.values.tolist()[3:-1]

    # plot_correlations(training_data, ['Mean', 'Peak Frequency', 'Shannon Entropy', 'MFCC 1'])
    # plot_correlations(training_data, feature_names[0:10])
    # plot_correlations(training_data, feature_names[11:25])
    # plot_correlations(training_data, feature_names[26:46])
    # plot_correlations(training_data, feature_names)

    ## Compare different classifiers, including
        # ["Nearest Neighbors", "Linear SVM", "RBF SVM", "Gaussian Process",
        # "Decision Tree", "Random Forest", "Neural Net", "AdaBoost",
        # "Naive Bayes", "QDA"]

    numeric_class_labels = training_data['Class'].values
    feature_labels = training_data['ClassLabels'].values

    feature_vector = training_data.drop(['FeatureID', 'FileName','Class','ClassLabels'], axis=1)
    feature_vector.describe()

    count = len(numeric_class_labels) # 300
    dataset = [(feature_vector[['Energy', 'Spectral Spread']][0:count], numeric_class_labels[0:count]),
                (feature_vector[['Peak Frequency', 'STALTA']][0:count], numeric_class_labels[0:count]),
                (feature_vector[['Variance', 'Kurtosis']][0:count], numeric_class_labels[0:count])]

    #compare_classifiers(dataset)

    ## Train Classifier
    print('Training ['+ classifier_name + '] classifer in process[Waiting...]')
    classifier = create_classifier(classifier_name)
    classifier.set_data(training_data, class_labels)
    classifier.split_dataset(classifier.scaled_features)
    classifier.fit()
    classifier.predict()

    y_pred = classifier.clf.predict_proba(classifier.X_test)[:,1]
    plot_ROC(classifier.y_test, y_pred, pos_label = 2)
    #calc_auc(classifier.y_test, y_pred, pos_label = 2)

    # save trained classifier to disk
    file_name = os.path.join(model_dir, classifier_name + '_model_raw.sav')
    pickle.dump(classifier, open(file_name, 'wb'))

    # # visualize result
    # feature_names = ['Peak Frequency', 'Shannon Entropy']
    # visualize_ml_result(classifier_name, training_data, feature_names, class_labels, count = None)

    # PCA analysis
    n_components = 3
    print('\nPCA analysis results:\n')
    crossplot_pca(classifier.pca(n_components)[1], classifier.feature_labels)
    heatplot_pca(classifier.pca(n_components)[0], classifier.feature_names)

    classifier.split_dataset(classifier.pca(n_components)[1])
    classifier.fit()
    classifier.predict()
    y_pred = classifier.clf.predict_proba(classifier.X_test)[:,1]
    plot_ROC(classifier.y_test, y_pred, pos_label = 2)
    #calc_auc(classifier.y_test, y_pred, pos_label = 2)

    # save trained classifier to disk
    file_name = os.path.join(model_dir, classifier_name + '_model_pca.sav')
    pickle.dump(classifier, open(file_name, 'wb'))

    if classifier_name == 'SVM':

        # Train SVM classifier
        classifier.visualize_ml_result(['Peak Frequency', 'Shannon Entropy'])

        classifier.fit(kernel = 'linear')
        classifier.predict()
        # Plot coefficients
        classifier.plot_coefficients()

        # C_range = np.array([.01, 1, 5, 10, 20, 50, 100, 1000, 5000, 10000])
        # gamma_range = np.array([0.0001, 0.001, 0.01, 0.1, 1, 10])
        # classifier.model_param_selection(C_range, gamma_range)
        # classifier.fit_with_selected_model_param(10, 'auto')

    print('\nTraining completed[OK]')

# This will actually run this code if called stand-alone
if __name__ == '__main__':
    main()

reason: import numpy,from scipy