text:
import sys, os
import keras.backend.tensorflow_backend as ktf
import keras
import tensorflow as tf
import pandas as pd
import numpy as np
import time
import datetime
from scipy import interp
from datetime import datetime
import datetime as dt
from Segmentation.model_keras import *
from Segmentation.load_datasets import load_filenames_2nd, load_data, keep_t2
from sklearn.metrics import auc, roc_curve
from glob2 import glob
import nibabel as nib
import seaborn as sns
import json
import warnings
warnings.filterwarnings("ignore")

RAND_SEED = 42

# Disable stdout
def blockPrint():
    sys.stdout = open(os.devnull, 'w')

# Restore stdout
def enablePrint():
    sys.stdout = sys.__stdout__

def train_model(networks, spaths, only_t2):
    """
    Train specified model
    Args:
        networks (list): list of keras networks to train
        spaths (list): list of output directories
        only_t2 (bool): whether or not to only use T2 data

    Returns:

    """

    epochs = 600
    batch_size = 20
    block_size = [18, 142, 142]
    oversamp = 1.0
    oversamp_test = 1.0
    lab_trun = 2
    im_freq = 50
    val_split = 0.2
    test_split = 0.1
    lr = 1e-4

    adaptive_hist = False

    # Load training data
    image_base_path = '/media/matt/Seagate Expansion Drive/MR Data/MR_Images_Sarcoma'
    # image_base_path = 'E:/MR Data/MR_Images_Sarcoma'

    # Load training data
    filenames = load_filenames_2nd(base_path=image_base_path)
    nfiles = len(filenames)

    # Remove T2 images
    if only_t2:
        filenames = keep_t2(filenames)

    # Remove validation and test set
    inds = np.array((range(nfiles)), dtype=int)
    np.random.seed(RAND_SEED)
    np.random.shuffle(inds)

    # Validation data
    val_inds = inds[:round(val_split*nfiles)]
    val_file = [filenames[i] for i in val_inds]

    # Test data
    test_inds = inds[-round(test_split*nfiles):]
    test_file = [filenames[i] for i in test_inds]

    # Delete all data
    filenames = [filename for i, filename in enumerate(filenames) if i not in
                 list(val_inds) + list(test_inds)]

    # Load data
    print('Loading data')
    x, y, orig_size = load_data(filenames, block_size, oversamp,
                                lab_trun, adaptive_hist)

    # val_file = val_file[:1]
    x_val, Y_val, orig_size_val = load_data(val_file, block_size, oversamp,
                                            lab_trun, adaptive_hist)

    # test_file = test_file[:1]
    x_test, y_test, orig_size_test = load_data(test_file, block_size,
                                               oversamp_test,
                                               lab_trun, adaptive_hist)

    # shuffle training data
    inds = np.arange(0, x.shape[0])
    np.random.seed(5)
    np.random.shuffle(inds)
    x = x[inds]
    y = y[inds]

    print('Size of training set:\t\t', x.shape)
    print('Size of validation set: \t', x_val.shape)
    print('Size of test set: \t\t', x_test.shape)

    sz_patch = x.shape

    for network, spath in zip(networks, spaths):

        # Set up save path
        spath = spath % datetime.strftime(datetime.now(), '%Y_%m_%d_%H-%M-%S')

        if not os.path.exists(spath):
            os.mkdir(spath)

        # Display which network is training
        _, net_name = os.path.split(spath)
        if only_t2:
            net_name += '_t2'
        print('\n\n\n')
        print('Training: %s' % net_name)
        print('-' * 80 + '\n')

        # Save a copy of this code
        save_code(spath, os.path.realpath(__file__))

        # Load model
        model, opt = network(pretrained_weights=None,
                             input_size=(sz_patch[1],
                                         sz_patch[2],
                                         sz_patch[3],
                                         sz_patch[4]),
                             lr=lr)

        # Set up callbacks
        tensorboard = keras.callbacks.TensorBoard(log_dir="logs/%s_%s"
                                                  % net_name,
                                                  write_graph=False,
                                                  write_grads=False,
                                                  write_images=False,
                                                  histogram_freq=0)
        ckpoint_weights = keras.callbacks.ModelCheckpoint(os.path.join(spath, 'ModelCheckpoint.h5'),
                                                          monitor='val_loss',
                                                          verbose=1,
                                                          save_best_only=True,
                                                          save_weights_only=True,
                                                          mode='auto',
                                                          period=10)

        image_recon_callback = image_callback_val(x_val, Y_val, spath,
                                                  orig_size_val,
                                                  block_size=block_size,
                                                  oversamp=oversamp_test,
                                                  lab_trun=lab_trun,
                                                  im_freq=im_freq,
                                                  batch_size=batch_size)

        # Train model
        model.fit(x=x, y=y,
                  epochs=epochs,
                  batch_size=batch_size,
                  validation_data=(x_val, Y_val),
                  callbacks=[tensorboard,
                             ckpoint_weights,
                             image_recon_callback
                             ]
                  )

        model.save(os.path.join(spath, 'Trained_model.h5'))

        # Calculate best threshold from training data
        threshold = training_threshold(model, spath, X=x, Y=y)
        # threshold = 0.5

        # Evaluate test data
        test_set_3D(model, x_test, y_test, spath, orig_size_test, block_size,
                    oversamp_test, lab_trun, batch_size=1, threshold=threshold)

def run_model_test(spaths, only_t2):
    """
    Test trained network
    Args:
        spaths (list): list of paths which contain the trained models
        only_t2 (bool): whether or not to only use T2 data

    Returns:

    """

    epochs = 600
    batch_size = 20
    block_size = [18, 142, 142]
    oversamp = 1.0
    oversamp_test = 1.0
    lab_trun = 2
    im_freq = 50
    val_split = 0.2
    test_split = 0.1
    lr = 2e-4

    adaptive_hist = False

    # Load training data
    image_base_path = '/media/matt/Seagate Expansion Drive/MR Data/MR_Images_Sarcoma'

    # Load training data
    filenames = load_filenames_2nd(base_path=image_base_path)
    nfiles = len(filenames)

    # Remove all but T2 images
    if only_t2:
        filenames = keep_t2(filenames)

    # Remove validation and test set
    inds = np.array((range(nfiles)), dtype=int)
    np.random.seed(RAND_SEED)
    np.random.shuffle(inds)

    # Validation data
    val_inds = inds[:round(val_split*nfiles)]
    val_file = [filenames[i] for i in val_inds]

    # Test data
    test_inds = inds[-round(test_split*nfiles):]
    test_file = [filenames[i] for i in test_inds]

    # Delete all data
    filenames = [filename for i, filename in enumerate(filenames) if i not in
                 list(val_inds) + list(test_inds)]

    # Load data
    x, y, orig_size = load_data(filenames, block_size, oversamp,
                                lab_trun, adaptive_hist)

    x_test, y_test, orig_size_test = load_data(test_file, block_size,
                                               oversamp_test,
                                               lab_trun, adaptive_hist)
    print('Size of training set:\t\t', x.shape)
    print('Size of test set: \t\t', x_test.shape)

    for spath in spaths:

        # Display which network is training
        _, net_name = os.path.split(spath)
        print('\n\n\n')
        print('Testing: %s' % net_name)
        print('-' * 80 + '\n')

        # Load trained model
        model_path = os.path.join(spath, 'Trained_model.h5')
        model = keras.models.load_model(model_path,
                                        custom_objects={'dice_loss': dice_loss,
                                                        'dice_metric': dice_metric})

        # Calculate best threshold from training data
        threshold = training_threshold(model, spath, X=x, Y=y)
        # threshold = 0.5

        # Evaluate test data
        test_set_3D(model, x_test, y_test, spath, orig_size_test, block_size,
                    oversamp_test, lab_trun, batch_size=20, threshold=threshold, vols=len(test_file), continuous=True)

def run_model_test_best_weigts(spaths, only_t2):

    epochs = 600
    batch_size = 20
    block_size = [18, 142, 142]
    oversamp = 1.0
    oversamp_test = 1.0
    lab_trun = 2
    im_freq = 50
    val_split = 0.2
    test_split = 0.1
    lr = 2e-4

    adaptive_hist = False

    # Load training data
    image_base_path = '/media/matt/Seagate Expansion Drive/MR Data/MR_Images_Sarcoma'

    # Load training data
    filenames = load_filenames_2nd(base_path=image_base_path)
    nfiles = len(filenames)

    # Remove all but T2 images
    if only_t2:
        filenames = keep_t2(filenames)

    # Remove validation and test set
    inds = np.array((range(nfiles)), dtype=int)
    np.random.seed(RAND_SEED)
    np.random.shuffle(inds)

    # Validation data
    val_inds = inds[:round(val_split*nfiles)]
    val_file = [filenames[i] for i in val_inds]

    # Test data
    test_inds = inds[-round(test_split*nfiles):]
    test_file = [filenames[i] for i in test_inds]

    # Delete all data
    filenames = [filename for i, filename in enumerate(filenames) if i not in
                 list(val_inds) + list(test_inds)]

    # Load data
    x, y, orig_size = load_data(filenames, block_size, oversamp,
                                lab_trun, adaptive_hist)

    x_test, y_test, orig_size_test = load_data(test_file, block_size,
                                               oversamp_test,
                                               lab_trun, adaptive_hist)
    print('Size of training set:\t\t', x.shape)
    print('Size of test set: \t\t', x_test.shape)

    for spath in spaths:

        # Display which network is training
        _, net_name = os.path.split(spath)
        print('\n\n\n')
        print('Testing: %s' % net_name)
        print('-' * 80 + '\n')

        # Load trained model
        model_path = os.path.join(spath, 'Trained_model.h5')
        model = keras.models.load_model(model_path,
                                        custom_objects={'dice_loss': dice_loss,
                                                        'dice_metric': dice_metric})

        # Load best (from validation set)
        weights_file = os.path.join(spath, 'ModelCheckpoint.h5')
        model.load_weights(weights_file)

        # Calculate best threshold from training data
        threshold = training_threshold(model, spath, X=x, Y=y)
        # threshold = 0.5

        # Evaluate test data
        test_set_3D(model, x_test, y_test, spath, orig_size_test, block_size,
                    oversamp_test, lab_trun, batch_size=20, threshold=threshold, vols=len(test_file), continuous=True)

def train_networks():
    """
    Sets up all networks to train
    Returns:

    """

    # Save locations
    save_base_path = '/media/matt/Seagate Expansion Drive/MR Data/ML_Results'

    # List of networks to train
    networks = [cnn_model_3D_3lyr_relu_dice,
                cnn_model_3D_3lyr_do_relu_dice_skip,
                #cnn_model_3D_FlDense_dice,
                #cnn_model_3D_3lyr_do_relu_hing,
                cnn_model_3D_3lyr_do_relu_xentropy,
                cnn_model_3D_3lyr_do_relu_xentropy_skip
                ]

    spaths = ['%s_' + i.__name__ for i in networks]
    spaths = [os.path.join(save_base_path, path) for path in spaths]

    # Train the model
    train_model(networks, spaths, only_t2=False)

    # T2 - only
    # List of networks to train
    networks = [cnn_model_3D_3lyr_relu_dice,
                cnn_model_3D_3lyr_do_relu_dice_skip,
                # cnn_model_3D_FlDense_dice,
                # cnn_model_3D_3lyr_do_relu_hing,
                cnn_model_3D_3lyr_do_relu_xentropy,
                cnn_model_3D_3lyr_do_relu_xentropy_skip
                ]

    spaths = ['%s_t2_' + i.__name__ for i in networks]
    spaths = [os.path.join(save_base_path, path) for path in spaths]

    # Train the model
    train_model(networks, spaths, only_t2=True)

def run_networks_test():
    """
    Tests all trained networks.
    Returns:

    """
    # Paths to trained networks
    spaths = ['/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_08_14-36-46_cnn_model_3D_3lyr_relu_dice',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_08_21-50-21_cnn_model_3D_3lyr_do_relu_dice_skip',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_09_06-49-45_cnn_model_3D_3lyr_do_relu_xentropy',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_09_14-12-47_cnn_model_3D_3lyr_do_relu_xentropy_skip',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_09_23-04-28_t2_cnn_model_3D_3lyr_relu_dice',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_10_04-50-05_t2_cnn_model_3D_3lyr_do_relu_dice_skip',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_10_12-28-23_t2_cnn_model_3D_3lyr_do_relu_xentropy',
              '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_10_18-43-24_t2_cnn_model_3D_3lyr_do_relu_xentropy_skip'
              ]

    # run_model_test(spaths[:4], only_t2=False)
    # run_model_test(spaths[4:], only_t2=True)

    run_model_test_best_weigts(spaths[:4], only_t2=False)
    run_model_test_best_weigts(spaths[4:], only_t2=True)

def load_local_test_data(filenames):

    X = np.empty(shape=(0, 0, 0, 0, 0))
    for file in filenames:

        tmpx = nib.load(file).get_data().astype(np.float).squeeze()

        try:
            X = np.concatenate((X, tmpx), axis=2)
        except ValueError:
            X = tmpx

    return X

def load_train_volumes(only_t2=False, adaptive_hist=False):

    # Set up image path
    image_base_path = '/media/matt/Seagate Expansion Drive/MR Data/MR_Images_Sarcoma'

    # Set up data constants
    block_size = [18, 142, 142]
    oversamp_test = 1.0
    lab_trun = 2
    test_split = 0.1
    val_split = 0.2

    # Get filenames
    filenames = load_filenames_2nd(base_path=image_base_path)
    nfiles = len(filenames)

    # Yield the number of sets in the generator
    yield round((1 - (val_split + test_split) ) * nfiles)

    if only_t2:
        filenames = keep_t2(filenames)

    # Remove validation and test set
    inds = np.array((range(nfiles)), dtype=int)
    np.random.seed(RAND_SEED)
    np.random.shuffle(inds)
    mask = np.ones(inds.shape, dtype=bool)

    # Test data
    mask[:round(val_split * nfiles)] = 0
    mask[-round(test_split*nfiles):] = 0
    train_files = [filenames[i] for i in train_inds]

    while True:

        for train_file in train_files:

            X_train, Y_train, _ = load_data([train_file], block_size,
                                                       oversamp_test,
                                                       lab_trun, adaptive_hist)

            yield [X_train, Y_train]

def load_models(paths):
    """
    Loads a list of models
    Args:
        paths (list): list of paths to models (not including the filename)

    Returns:

    """

    model = []

    for path in paths:

        model_name = os.path.join(path, 'Trained_model.h5')

        model.append(keras.models.load_model(model_name,
                                             custom_objects=
                                             {'dice_loss': dice_loss}))

    return model

def read_tensorboard(path, fields, label):
    """
    Reads tensorboard files and returns data requests as a pandas array
    Args:
        path (str): path to tensorboard file
        fields (list of str): fields to return
        label (str): dataset identifier

    Returns:
        pandas dataframe
    """

    # Update fields with validation data
    fields = [[i] + ['val_{}'.format(i)] for i in fields]
    fields = sum(fields, [])

    # Set up output dictionary
    df = dict()
    for field in fields:
        df[field] = list()

    # Read the file
    for e in tf.train.summary_iterator(path):
        for v in e.summary.value:
            for field in fields:
                if field == v.tag:
                    df[field].append(v.simple_value)

    # Convert to pandas
    df = pd.DataFrame.from_dict(df)

    # Append set label
    df['label'] = label
    df['Epoch'] = list(range(1, len(df)+1))

    return df

def training_curves(spath):
    """
    Load and plot training curves.
    Args:
        spath (str): path on which to save output

    Returns:

    """

    # Location with saved curves
    dat_path = ['/home/matt/Documents/SegSarcoma/logs/2019_11_08_14-36-47_2019_11_08_14-36-46_cnn_model_3D_3lyr_relu_dice',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_08_21-50-22_2019_11_08_21-50-21_cnn_model_3D_3lyr_do_relu_dice_skip',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_09_06-49-45_2019_11_09_06-49-45_cnn_model_3D_3lyr_do_relu_xentropy',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_09_14-12-48_2019_11_09_14-12-47_cnn_model_3D_3lyr_do_relu_xentropy_skip',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_09_23-04-28_2019_11_09_23-04-28_t2_cnn_model_3D_3lyr_relu_dice_t2',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_10_04-50-05_2019_11_10_04-50-05_t2_cnn_model_3D_3lyr_do_relu_dice_skip_t2',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_10_12-28-24_2019_11_10_12-28-23_t2_cnn_model_3D_3lyr_do_relu_xentropy_t2',
                '/home/matt/Documents/SegSarcoma/logs/2019_11_10_18-43-24_2019_11_10_18-43-24_t2_cnn_model_3D_3lyr_do_relu_xentropy_skip_t2'
                ]
    labels = ['DICE', 'DICE - skip', 'XEntr', 'XEntr - skip', 'Dice T2', 'Dice T2 - skip', 'XEntr T2', 'XEntr T2 - skip']
    fields = ['binary_accuracy', 'dice_metric', 'loss']

    # Initialize dataframe
    df = pd.DataFrame()

    # Load data
    for dpath, label in zip(dat_path, labels):

        # Get tensorboard file
        file = glob(os.path.join(dpath, 'events*'))[0]

        # Read data
        tmp = read_tensorboard(file, fields, label)

        # Append data
        df = pd.concat((df, tmp))

    # Save path
    path = '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/JournalPaper/TrainingCurves'

    # Plot with and without skip connections
    sns.set('paper')
    sns.set_style('whitegrid')
    fs = 12
    lw = 2.5

    # Labels
    labs = ['DICE', 'DICE Skip', 'XEntropy', 'XEntropy - Skip']
    colors = ['b', 'r', 'g', 'k']

    snames = os.path.join(path, 'Acc.svg')

    plt.close()
    fig, ax = plt.subplots()
    for z in range(len(labels) // 2):
        inds = df['label'] == labels[z]
        ax.plot(df.loc[inds]['binary_accuracy'], '-', color=colors[z], lw=lw, label=labels[z])
        ax.plot(df.loc[inds]['val_binary_accuracy'], ':', color=colors[z], lw=lw, label='Validation')

        ax.set_xlabel('Epoch', fontsize=fs)
        ax.set_ylabel('Accuracy', fontsize=fs)
        ax.set_ylim([0.9, 1.0])
        # ax.set_yticks(np.linspace(0.9, 1.0, 3))
        ax.legend(fontsize=fs-1)
        for item in list(ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fs-2)
    fig.savefig(os.path.join(path, snames))

    # DICE
    snames = os.path.join(path, 'Dice_metric.svg')

    fig, ax = plt.subplots()
    for z in range(len(labels) // 2):
        inds = df['label'] == labels[z]
        ax.plot(df.loc[inds]['dice_metric'], '-', color=colors[z], lw=lw, label=labels[z])
        ax.plot(df.loc[inds]['val_dice_metric'], ':', color=colors[z], lw=lw, label='Validation')
        ax.set_xlabel('Epoch', fontsize=fs)
        ax.set_ylabel('DICE', fontsize=fs)
        ax.set_ylim([0.35, 1.0])
        ax.legend(fontsize=fs-1, loc=4)
        for item in list(ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fs-2)
    fig.savefig(os.path.join(path, snames))

    # Loss
    snames = os.path.join(path, 'Loss.svg')

    fig, ax = plt.subplots()
    for z in range(len(labels) // 2):
        inds = df['label'] == labels[z]
        ax.plot(df.loc[inds]['loss'], '-', color=colors[z], lw=lw, label=labels[z])
        ax.plot(df.loc[inds]['val_loss'], ':', color=colors[z], lw=lw, label='Validation')
        ax.set_xlabel('Epoch', fontsize=fs)
        ax.set_ylabel('Loss', fontsize=fs)
        ax.set_ylim([0.0, 0.65])
        ax.legend(fontsize=fs-1, loc=1)
        for item in list(ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fs-2)
    fig.savefig(os.path.join(path, snames))

    ## T2 ##
    snames = os.path.join(path, 'Acc_t2.svg')

    plt.close()
    fig, ax = plt.subplots()
    for z in range(len(labels) // 2, len(labels)):
        inds = df['label'] == labels[z]
        ax.plot(df.loc[inds]['binary_accuracy'], '-', color=colors[z-len(labels)//2], lw=lw, label=labels[z])
        ax.plot(df.loc[inds]['val_binary_accuracy'], ':', color=colors[z-len(labels)//2], lw=lw, label='Validation')

        ax.set_xlabel('Epoch', fontsize=fs)
        ax.set_ylabel('Accuracy', fontsize=fs)
        ax.set_ylim([0.9, 1.0])
        # ax.set_yticks(np.linspace(0.9, 1.0, 3))
        ax.legend(fontsize=fs - 1)
        for item in list(ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fs - 2)
    fig.savefig(os.path.join(path, snames))

    # DICE
    snames = os.path.join(path, 'Dice_metric_t2.svg')

    fig, ax = plt.subplots()
    for z in range(len(labels) // 2, len(labels)):
        inds = df['label'] == labels[z]
        ax.plot(df.loc[inds]['dice_metric'], '-', color=colors[z-len(labels)//2], lw=lw, label=labels[z])
        ax.plot(df.loc[inds]['val_dice_metric'], ':', color=colors[z-len(labels)//2], lw=lw, label='Validation')
        ax.set_xlabel('Epoch', fontsize=fs)
        ax.set_ylabel('DICE', fontsize=fs)
        ax.set_ylim([0.35, 1.0])
        ax.legend(fontsize=fs - 1, loc=4)
        for item in list(ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fs - 2)
    fig.savefig(os.path.join(path, snames))

    # Loss
    snames = os.path.join(path, 'Loss_t2.svg')

    fig, ax = plt.subplots()
    for z in range(len(labels) // 2, len(labels)):
        inds = df['label'] == labels[z]
        ax.plot(df.loc[inds]['loss'], '-', color=colors[z-len(labels)//2], lw=lw, label=labels[z])
        ax.plot(df.loc[inds]['val_loss'], ':', color=colors[z-len(labels)//2], lw=lw, label='Validation')
        ax.set_xlabel('Epoch', fontsize=fs)
        ax.set_ylabel('Loss', fontsize=fs)
        ax.set_ylim([0.0, 0.65])
        ax.legend(fontsize=fs - 1, loc=1)
        for item in list(ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(fs - 2)
    fig.savefig(os.path.join(path, snames))

def compare_dice():
    """
    Load and make plots of Dice score for each network trained.
    Returns:

    """

    folders = ['/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_08_14-36-46_cnn_model_3D_3lyr_relu_dice',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_08_21-50-21_cnn_model_3D_3lyr_do_relu_dice_skip',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_09_06-49-45_cnn_model_3D_3lyr_do_relu_xentropy',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_09_14-12-47_cnn_model_3D_3lyr_do_relu_xentropy_skip',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_09_23-04-28_t2_cnn_model_3D_3lyr_relu_dice',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_10_04-50-05_t2_cnn_model_3D_3lyr_do_relu_dice_skip',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_10_12-28-23_t2_cnn_model_3D_3lyr_do_relu_xentropy',
               '/media/matt/Seagate Expansion Drive/MR Data/ML_Results/2019_11_10_18-43-24_t2_cnn_model_3D_3lyr_do_relu_xentropy_skip']

    labels = ['DICE', 'DICE - skip', 'XEntropy', 'XEntropy - skip', 'Dice T2', 'Dice T2 - skip', 'XEntropy T2', 'Xentropy T2 - skip']

    file_name = 'metrics2.txt'

    df = {'Net':[], 'Threshold': [], 'Precision': [], 'Recall': [], 'ROC': [], 'DICE': [], 'VOE': []}

    keys = ['Threshold', 'ROC', 'Precision', 'Recall', 'DICE', 'VOE']

    for folder, label in zip(folders, labels):

        file = os.path.join(folder, file_name)

        df['Net'].append(label)

        # Read meetrics file
        with open(file, 'r') as f:
            dat = f.readlines()

        # Append values
        for ind, k in zip(range(-7, -1), keys):

            tmp = [i for i in dat[ind] if i.isdigit() or i == '.']
            df[k].append(float(''.join(tmp)))

    # Compute VOE
    from sklearn.metrics import jaccard_similarity_score

    # Convert to Pandas
    df = pd.DataFrame.from_dict(df)
    df = df.reindex([0, 1, 4, 5, 2, 3, 6, 7])

    # Write to CSV
    df.to_csv('/media/matt/Seagate Expansion Drive/MR Data/ML_Results/JournalPaper/metrics-val.csv')

    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_axes([0.1, 0.2, 0.80, 0.47])

    g = sns.barplot(x='Net', y='DICE', data=df, ax=ax, palette='Set1')
    plt.xticks(rotation=45)
    fig.savefig('/media/matt/Seagate Expansion Drive/MR Data/ML_Results/JournalPaper/DICE_nets.svg')

    plt.close(fig)

def statistical_metrics(spaths):

    epochs = 600
    batch_size = 20
    block_size = [18, 142, 142]
    oversamp = 1.0
    oversamp_test = 1.0
    lab_trun = 2
    im_freq = 50
    val_split = 0.2
    test_split = 0.1
    lr = 2e-4

    # Load training data
    image_base_path = '/media/matt/Seagate Expansion Drive/MR Data/MR_Images_Sarcoma'

    # Load training data
    filenames = load_filenames_2nd(base_path=image_base_path)
    nfiles = len(filenames)

    # Remove all but T2 images
    if only_t2:
        filenames = keep_t2(filenames)

    # Remove validation and test set
    inds = np.array((range(nfiles)), dtype=int)
    np.random.seed(RAND_SEED)
    np.random.shuffle(inds)

    # Validation data
    val_inds = inds[:round(val_split*nfiles)]
    val_file = [filenames[i] for i in val_inds]

    # Test data
    test_inds = inds[-round(test_split*nfiles):]
    test_file = [filenames[i] for i in test_inds]

    # Delete all data
    filenames = [filename for i, filename in enumerate(filenames) if i not in
                 list(val_inds) + list(test_inds)]

    # Load data
    x_test, y_test, orig_size_test = load_data(test_file, block_size,
                                               oversamp_test,
                                               lab_trun, adaptive_hist)
    print('Size of test set: \t\t', x_test.shape)

    for spath in spaths:

        # Display which network is training
        _, net_name = os.path.split(spath)
        print('\n\n\n')
        print('Testing: %s' % net_name)
        print('-' * 80 + '\n')

        # Load trained model
        model_path = os.path.join(spath, 'Trained_model.h5')
        model = keras.models.load_model(model_path,
                                        custom_objects={'dice_loss': dice_loss,
                                                        'dice_metric': dice_metric})

        # Load best threshold
        file = os.path.join(spath, 'metrics2.txt')

        # Read meetrics file
        with open(file, 'r') as f:
            dat = f.readlines()

        # Append values
        ind = -7
        tmp = [i for i in dat[ind] if i.isdigit() or i == '.']
        threshold = float(''.join(tmp))

        # Evaluate test data

if __name__ == '__main__':
    tstart = time.time()
    train_networks()
    run_networks_test()
    training_curves(spath='/media/matt/Seagate Expansion Drive/MR Data/ML_Results/JournalPaper/TrainingCurves')
    compare_dice()

    print('\tTotal time (HH:MM:SS): %s\n\n' % (str(dt.timedelta(seconds=round(time.time() - tstart)))))


text:
#-*- coding: utf-8 -*-
"""
All modification made by Cambricon Corporation: © 2018 Cambricon Corporation
All rights reserved.
All other contributions:
Copyright (c) 2014--2018, the respective contributors
All rights reserved.
For the list of contributors go to https://github.com/BVLC/caffe/blob/master/CONTRIBUTORS.md
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""

import unittest
import tempfile
import os
import numpy as np
import six

import caffe
from test_net import simple_net_file

class TestSolver(unittest.TestCase):
    def setUp(self):
        self.num_output = 13
        net_f = simple_net_file(self.num_output)
        f = tempfile.NamedTemporaryFile(mode='w+', delete=False)
        f.write("""net: '""" + net_f + """'
        test_iter: 10 test_interval: 10 base_lr: 0.01 momentum: 0.9
        weight_decay: 0.0005 lr_policy: 'inv' gamma: 0.0001 power: 0.75
        display: 100 max_iter: 100 snapshot_after_train: false
        snapshot_prefix: "model" """)
        f.close()
        self.solver = caffe.SGDSolver(f.name)
        # also make sure get_solver runs
        caffe.get_solver(f.name)
        caffe.set_mode_cpu()
        # fill in valid labels
        self.solver.net.blobs['label'].data[...] = \
                np.random.randint(self.num_output,
                    size=self.solver.net.blobs['label'].data.shape)
        self.solver.test_nets[0].blobs['label'].data[...] = \
                np.random.randint(self.num_output,
                    size=self.solver.test_nets[0].blobs['label'].data.shape)
        os.remove(f.name)
        os.remove(net_f)

    def test_solve(self):
        self.assertEqual(self.solver.iter, 0)
        self.solver.solve()
        self.assertEqual(self.solver.iter, 100)

    def test_net_memory(self):
        """Check that nets survive after the solver is destroyed."""

        nets = [self.solver.net] + list(self.solver.test_nets)
        self.assertEqual(len(nets), 2)
        del self.solver

        total = 0
        for net in nets:
            for ps in six.itervalues(net.params):
                for p in ps:
                    total += p.data.sum() + p.diff.sum()
            for bl in six.itervalues(net.blobs):
                total += bl.data.sum() + bl.diff.sum()

    def test_snapshot(self):
        self.solver.snapshot()
        # Check that these files exist and then remove them
        files = ['model_iter_0.caffemodel', 'model_iter_0.solverstate']
        for fn in files:
            assert os.path.isfile(fn)
            os.remove(fn)

ext: py

max_stars_repo_name: Julian-He/caffe

max_issues_repo_name: Julian-He/caffe

max_forks_repo_name: Julian-He/caffe

reason: import numpy


text:
import numpy as np

def cross_entropy_error(y, t):
    if y.ndim == 1:
        t = t.reshape(1, t.size)
        y = y.reshape(1, y.size)

        batch_size = y.shape[0]
        return -np.sum(t * np.log(y)) / batch_size

max_stars_repo_name: BlueWay-KU/Study

max_issues_repo_name: BlueWay-KU/Study

max_forks_repo_name: BlueWay-KU/Study


text:
import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.nn.functional as F
import random
import numpy as np

from layers import GCN, AvgReadout, Discriminator, Discriminator2, Discriminator_jsd, MeanAggregator, MeanAggregator_ml, Encoder
from estimator.estimator import mi_loss
from attacker.attacker import Attacker

class DGI(nn.Module):
    def __init__(self, n_in, n_h, activation, critic="bilinear", dataset=None, attack_model=True):
        super(DGI, self).__init__()
        self.gcn = GCN(n_in, n_h, activation)
        self.read = AvgReadout()

        self.sigm = nn.Sigmoid()

        self.disc = Discriminator(n_h, critic=critic, dataset=dataset, attack_model=attack_model)

    def forward(self, seq1, seq2, adj, sparse, msk, samp_bias1, samp_bias2):
        h_1 = self.gcn(seq1, adj, sparse)   # (1, 2708, 512)

        c = self.read(h_1, msk)   # (1, 512)
        c = self.sigm(c)

        h_2 = self.gcn(seq2, adj, sparse)

        ret = self.disc(c, h_1, h_2, samp_bias1, samp_bias2)   # (1, 5416)

        return ret

    # Detach the return variables
    def embed(self, seq, adj, sparse, msk, grad=False):
        h_1 = self.gcn(seq, adj, sparse)
        c = self.read(h_1, msk)

        if grad:
            return h_1, c
        else:
            return h_1.detach(), c.detach()

class DGI_local(nn.Module):
    def __init__(self, n_in, n_h, activation):
        super(DGI_local, self).__init__()
        self.gcn = GCN(n_in, n_h, activation)
        self.read = AvgReadout()

        self.sigm = nn.Sigmoid()

        self.disc1 = Discriminator(n_h)
        self.disc2 = Discriminator2(n_in, n_h)

    def forward(self, seq1, seq2, adj, sparse, msk, samp_bias1, samp_bias2):
        h_1 = self.gcn(seq1, adj, sparse)   # (1, 2708, 512)

        c = self.read(h_1, msk)   # (1, 512)
        c = self.sigm(c)

        h_2 = self.gcn(seq2, adj, sparse)

        ret1 = self.disc1(c, h_1, h_2, samp_bias1, samp_bias2)   # (1, 5416)
        # ret2 = self.disc2(h_1, seq1, seq2, samp_bias1, samp_bias2)  # (1, 5416)
        # ret = ret1 + ret2
        return ret1

    # Detach the return variables
    def embed(self, seq, adj, sparse, msk):
        h_1 = self.gcn(seq, adj, sparse)
        c = self.read(h_1, msk)

        return h_1.detach(), c.detach()

class DGI_ind(nn.Module):
    def __init__(self, features, adj_lists, ft_size, n_h, activation, num_sample=[10, 10], skip_connection=False, gcn=True):
        super(DGI_ind, self).__init__()
        self.features = features
        self.skip_connection = skip_connection
        self.agg1 = MeanAggregator(features, cuda=torch.cuda.is_available(), gcn=gcn, name='l1')
        self.enc1 = Encoder(features, ft_size, n_h, adj_lists, self.agg1, num_sample=num_sample[0], gcn=gcn,
                            cuda=torch.cuda.is_available(), activation=activation, skip_connection=skip_connection, name='l2')
        self.agg2 = MeanAggregator(lambda nodes: self.enc1(nodes), cuda=torch.cuda.is_available(), gcn=gcn, name='l3')
        self.enc2 = Encoder(lambda nodes: self.enc1(nodes), self.enc1.embed_dim, n_h, adj_lists, self.agg2, num_sample=num_sample[1],
                            base_model=self.enc1, gcn=gcn, cuda=torch.cuda.is_available(),
                            activation=activation, skip_connection=skip_connection, name='l4')
        self.read = AvgReadout()
        self.sigm = nn.Sigmoid()

        if skip_connection:
            self.disc = Discriminator(2*n_h)
        else:
            self.disc = Discriminator(n_h)
        # self.disc2 = Discriminator2(ft_size, 2*n_h)
        # self.disc = Discriminator_jsd(2*n_h)

    def forward(self, nodes, msk, samp_bias1, samp_bias2, make_adv=True):
        h_1, h_2 = self.enc2(nodes)

        h_1 = h_1.transpose(0,1).unsqueeze(0)
        h_2 = h_2.transpose(0,1).unsqueeze(0)

        c = self.read(h_1, msk)   # (1, 512)
        c = self.sigm(c)

        # x_1 = self.features[nodes]
        # idx = np.random.permutation(x_1.shape[0])
        # x_2 = x_1[idx, :]

        ret = self.disc(c, h_1, h_2, samp_bias1, samp_bias2)   # (1, 5416)
        # ret2 = self.disc2(h_1, x_1.unsqueeze(0), x_2.unsqueeze(0), samp_bias1, samp_bias2)  # (1, 5416)
        # ret = ret1 + ret2
        # ret = self.disc()
        return ret

    # Detach the return variables
    def embed(self, nodes, msk):
        h_1, _ = self.enc2(nodes)
        c = self.read(h_1, msk)

        return h_1.detach(), c.detach()

class DGI_ind_ml(nn.Module):
    def __init__(self, features, adj_lists, ft_size, n_h, activation, num_sample=[10, 10, 25], skip_connection=False, gcn=True):
        super(DGI_ind_ml, self).__init__()
        self.features = features
        self.skip_connection = skip_connection
        self.agg1 = MeanAggregator_ml(features, cuda=torch.cuda.is_available(), gcn=gcn, name='l1')
        self.enc1 = Encoder(features, ft_size, n_h, adj_lists, self.agg1, num_sample=num_sample[0], gcn=gcn,
                            cuda=torch.cuda.is_available(), activation=activation, skip_connection=skip_connection, name='l2')
        self.agg2 = MeanAggregator_ml(lambda nodes: self.enc1(nodes), cuda=torch.cuda.is_available(), gcn=gcn, name='l3')
        self.enc2 = Encoder(lambda nodes: self.enc1(nodes), self.enc1.embed_dim, n_h, adj_lists, self.agg2, num_sample=num_sample[1],
                            base_model=self.enc1, gcn=gcn, cuda=torch.cuda.is_available(),
                            activation=activation, skip_connection=skip_connection, name='l4')
        self.agg3 = MeanAggregator_ml(lambda nodes: self.enc2(nodes), cuda=torch.cuda.is_available(), gcn=gcn, name='l5')
        self.enc3 = Encoder(lambda nodes: self.enc2(nodes), self.enc2.embed_dim, n_h, adj_lists, self.agg3, num_sample=num_sample[2],
                            base_model=self.enc2, gcn=gcn, cuda=torch.cuda.is_available(),
                            activation=activation, skip_connection=skip_connection, name='l6')
        self.read = AvgReadout()
        self.sigm = nn.Sigmoid()

        if skip_connection:
            self.disc = Discriminator(2*n_h)
        else:
            self.disc = Discriminator(n_h)
        # self.disc2 = Discriminator2(ft_size, 2*n_h)
        # self.disc = Discriminator_jsd(2*n_h)

    def forward(self, nodes, msk, samp_bias1, samp_bias2, make_adv=True):
        h_1, h_2 = self.enc3(nodes)

        h_1 = h_1.transpose(0,1).unsqueeze(0)
        h_2 = h_2.transpose(0,1).unsqueeze(0)

        c = self.read(h_1, msk)   # (1, 512)
        c = self.sigm(c)

        ret = self.disc(c, h_1, h_2, samp_bias1, samp_bias2)   # (1, 5416)
        return ret

    # Detach the return variables
    def embed(self, nodes, msk):
        h_1, _ = self.enc2(nodes)
        c = self.read(h_1, msk)

        return h_1.detach(), c.detach()

max_stars_repo_name: Guo-lab/Graph

max_issues_repo_name: Guo-lab/Graph

max_forks_repo_name: Guo-lab/Graph


text:
import argparse

import cv2
import numpy as np
import torch
from numpy.linalg import norm as l2norm
from sklearn import preprocessing

from backbones import get_model

@torch.no_grad()
def inference(weight, name, img):
    device = torch.device("cpu")
    if img is None:
        img = np.random.randint(0, 255, size=(112, 112, 3), dtype=np.uint8)
    else:
        img = cv2.imread(img)
        img = cv2.resize(img, (112, 112))

    start = torch.cuda.Event(enable_timing=True)
    end = torch.cuda.Event(enable_timing=True)

    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    img = np.transpose(img, (2, 0, 1))
    img = torch.from_numpy(img).unsqueeze(0).float()
    img.div_(255).sub_(0.5).div_(0.5)
    net = get_model(name, fp16=False)
    net.load_state_dict(torch.load(weight))
    net.eval()
    net.to(device)
    img = img.to(device)
    for _ in range(10):
        _ = net(img)
    start.record()
    feat = net(img).cpu().numpy()
    end.record()
    torch.cuda.synchronize()

    print(start.elapsed_time(end))
    return feat

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='PyTorch ArcFace Training')
    parser.add_argument('--network', type=str, default='r50', help='backbone network')
    parser.add_argument('--weight', type=str, default='')

    args = parser.parse_args()
    feat1 = inference(args.weight, args.network, '/app/dataset_112x112/valid/85/85_51.jpg')
    feat2 = inference(args.weight, args.network, '/app/dataset_112x112/valid/85/85_566.jpg')
    norm_feat1 = feat1[0]/l2norm(feat1[0])
    norm_feat2 = feat2[0]/l2norm(feat2[0])
    print(np.dot(norm_feat1,norm_feat2))

reason: import numpy,from numpy


text:
import cv2
import numpy as np

import pyzbar.pyzbar as pyzbar

import cloudmersive_barcode_api_client
from cloudmersive_barcode_api_client.rest import ApiException

def decode(im) :
    # Find barcodes and QR codes
    decodedObjects = pyzbar.decode(im)
    return decodedObjects

def get_barcode(filepath):
    im              = cv2.imread(filepath)

    decodedObjects = None
    try:
        decodedObjects  = decode(im)
    except Exception as ex:
        print(ex)

    if decodedObjects is None:
        return None

    if len(decodedObjects) != 0:
        code        = str(decodedObjects[0].data)
        return code
    return None

# Display barcode and QR code location
def display(im, decodedObjects):

    # Loop over all decoded objects
    for decodedObject in decodedObjects:
        points = decodedObject.polygon

        # If the points do not form a quad, find convex hull
        if len(points) > 4 :
            hull = cv2.convexHull(np.array([point for point in points], dtype=np.float32))
            hull = list(map(tuple, np.squeeze(hull)))
        else :
            hull = points;

        # Number of points in the convex hull
        n = len(hull)

def decode_code(value):
    configuration = cloudmersive_barcode_api_client.Configuration()
    configuration.api_key['Apikey'] = '5f46158c-c65f-4199-bae4-df2b83b8c739'

    # create an instance of the API class
    api_instance = cloudmersive_barcode_api_client.BarcodeLookupApi(cloudmersive_barcode_api_client.ApiClient(configuration))
    #value = '5449000232250' # str | Barcode value

    api_response = None
    try:
        # Lookup a barcode value and return product data
        api_response = api_instance.barcode_lookup_ean_lookup(value)
    except ApiException as e:
        #print("Exception when calling BarcodeLookupApi->barcode_lookup_ean_lookup: %s\n" % e)
        pass
    return api_response

max_stars_repo_name: ZAurele/alpha-py

max_issues_repo_name: ZAurele/alpha-py

max_forks_repo_name: ZAurele/alpha-py


text:
import matplotlib.pyplot as plt
import numpy as np

from PIL import Image
from torchvision.transforms import transforms

def save_gradcam(filename, gcam, raw_image):
    gcam = gcam.squeeze().detach().cpu().numpy()
    cmap = plt.cm.jet_r(gcam)[..., :3] * 255.0
    cmap = cmap[:, :, ::-1]
    cmap = Image.fromarray(cmap.astype(np.uint8))
    raw_image = raw_image.permute(1, 2, 0).cpu().numpy() * 255
    transform = transforms.ToPILImage()
    raw_image = transform(raw_image.astype(np.uint8))
    result = Image.blend(raw_image, cmap, 0.5)
    result.save(filename)

max_stars_repo_name: ErikStammes/GradCAM-PyTorch

max_issues_repo_name: ErikStammes/GradCAM-PyTorch

max_forks_repo_name: ErikStammes/GradCAM-PyTorch


text:
import numpy as np
import matplotlib.pylab as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

datos = np.genfromtxt('datosonda.dat')

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1,2,1, projection='3d')
xx, yy = np.mgrid[0:datos.shape[0], 0:datos.shape[1]]
ax.plot_surface(xx, yy, datos ,rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0)

plt.xlabel("Tiempo[s]")
plt.ylabel("Posicion[m]")
plt.subplot(122)
plt.plot(datos[0:,], label="Tiempo incial")
plt.plot(datos[-1:,], label="Tiempo final")
plt.xlabel("Posicion[m]")
plt.ylabel("Desplazamiento[m]")
plt.legend(loc= 10)
plt.savefig("plotonda.png")


text:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

class data_transform:
    """A class similar to sklearn.utils.Bunch that I made for convenience.

    Examples:
        >>> import numpy as np
        >>> import pandas as pd
        >>> x0 = np.array(['a', 'b', 'a', 'a', 'c']).reshape(5, 1)
        >>> other_columns = np.random.rand(20).reshape(5, 4)
        >>> data = pd.DataFrame(np.hstack((x0, other_columns)),
                    columns=['x0', 'x1', 'x2', 'y0', 'y1'])
        >>> dt = data_transform(
                    df=data,
                    feature_names=['x0', 'x1', 'x2'],
                    target_names=['y0', 'y1'])
        >>> dt.target_names
        ['y0', 'y1']
        >>> dt.num_features = ['x0']
        >>> dt.get_cat_features()
        ['x1', 'x2']

    """

    def __init__(self, df, feature_names=None, target_names=None,
            num_features=None, cat_features=None):
        """[summary]

        Args:
            df ([type]): [description]
            feature_names (list, optional): [description]. Defaults to None.
            target_names (list, optional): [description]. Defaults to None.
            num_features (list, optional): [description]. Defaults to None.
            cat_features (list, optional): [description]. Defaults to None.
        """
        self.df = df
        self.feature_names = feature_names
        self.target_names = target_names
        self.num_features = num_features
        self.cat_features = cat_features

    def get_feature_names(self):
        """[summary]

        Returns:
            [type]: [description]
        """
        if self.target_names != None:
            cols = list(self.df.columns)
            for i in self.target_names:
                cols.remove(i)
            features = cols
            return features

        elif self.cat_features and self.num_features != None:
            features = self.num_features + self.cat_features
            return features

    def get_target_names(self):
        """[summary]

        Returns:
            [type]: [description]
        """
        if self.feature_names != None:
            cols = list(self.df.columns)
            for i in self.feature_names:
                cols.remove(i)
            self.targets = cols
            return self.targets

        elif self.cat_features and self.num_features != None:
            feature_names = self.get_feature_names()
            cols = list(self.df.columns)
            for j in feature_names:
                cols.remove(j)
            targets = cols
            return targets

    def get_num_features(self):
        """[summary]

        Returns:
            [type]: [description]
        """
        features = list(self.df.columns)
        # Remove target columns from feature list
        for j in self.target_names:
            features.remove(j)
        # Remove categircal features from feature list
        for i in self.cat_features:
            features.remove(i)
        # What's left over is the numerical features.
        return features

    def get_cat_features(self):
        """[summary]

        Returns:
            [type]: [description]
        """
        features = list(self.df.columns)
        # Remove target columns from feature list
        for j in self.target_names:
            features.remove(j)
        # Remove categircal features from feature list
        for i in self.num_features:
            features.remove(i)
        # What's left over is the categorical features.
        return features

    def one_hot_encode(self):
        """Converts an unstructued 2-D array consisting of categorical values
        into a 2-D array consisting of one-hot-encoded vectors.

        """
        X = self.df
        X_cat = X[self.cat_features]
        for cat in self.cat_features[:]:
            X = X.drop(cat, axis=1)

        # Replace the nonnumerical columns with one-hot encoded ones.
        for name in self.cat_features[:]:
            hot_one = pd.get_dummies(X_cat[name], prefix=name)
            X = pd.concat([X, hot_one.set_index(X.index)], axis=1)
        return X


text:
import torch
import numpy as np

class ToTensorGjz(object):
    def __call__(self, pic):
        img = torch.from_numpy(pic.transpose((2, 0, 1)))
        return img.float()

    def __repr__(self):
        return self.__class__.__name__ + '()'

class NormalizeGjz(object):
    def __init__(self, mean, std):
        self.mean = mean
        self.std = std

    def __call__(self, tensor):
        tensor.sub_(self.mean).div_(self.std)
        return tensor

def crop_img(img, roi_box):
    h, w = img.shape[:2]

    sx, sy, ex, ey = [int(round(_)) for _ in roi_box]
    dh, dw = ey - sy, ex - sx
    if len(img.shape) == 3:
        res = np.zeros((dh, dw, 3), dtype=np.uint8)
    else:
        res = np.zeros((dh, dw), dtype=np.uint8)
    if sx < 0:
        sx, dsx = 0, -sx
    else:
        dsx = 0

    if ex > w:
        ex, dex = w, dw - (ex - w)
    else:
        dex = dw

    if sy < 0:
        sy, dsy = 0, -sy
    else:
        dsy = 0

    if ey > h:
        ey, dey = h, dh - (ey - h)
    else:
        dey = dh

    res[dsy:dey, dsx:dex] = img[sy:ey, sx:ex]
    return res


text:
import os
import numpy as np

def save_model_generic(epoch, save_top_dir, saver, sess):
    # save_top_dir = os.path.join('../saved_models', model_name)
    if not os.path.exists(save_top_dir):
        os.makedirs(save_top_dir)
    ckpt_path = os.path.join(save_top_dir,str(epoch) + 'model.ckpt')
    if epoch == 0:
        saver.save(sess, ckpt_path, write_meta_graph=True)
    else:
        saver.save(sess, ckpt_path, write_meta_graph=False)
    print("Saving model at epoch %d to %s" % (epoch, ckpt_path))

def restore_model_generic(epoch, save_top_dir, saver, sess):
    ckpt_path = os.path.join(save_top_dir,str(epoch)+'model.ckpt')
    print("restoring from %s" % ckpt_path)
    saver.restore(sess, ckpt_path)

def write_tb(loss_dict, writer, cat_name, total_ct):
    for loss_name, loss_val in loss_dict.items():
        writer.add_scalar('{}/{}'.format(cat_name, loss_name), loss_val, total_ct)

def loss_dict_to_str(loss_dict):
    ret_str = ''
    for loss_name in sorted(loss_dict):
        loss_val = loss_dict[loss_name]
        ret_str += '{} {:.4f} '.format(loss_name, loss_val)
    return ret_str

def loss_dict_obj_cat_to_str(loss_dict, result_file_name_per_object_cat=None):
    ret_arr = []
    for loss_name in sorted(loss_dict):
        ret_str = loss_name + ' '
        for obj_cat in sorted(loss_dict[loss_name]):
            if not (result_file_name_per_object_cat is None):
                n_tmp = len(result_file_name_per_object_cat[obj_cat])
                ret_str += '{} {:.3f}({}) '.format(obj_cat, loss_dict[loss_name][obj_cat], n_tmp)
            else:
                ret_str += '{} {:.3f}'.format(obj_cat, loss_dict[loss_name][obj_cat])

        ret_arr.append(ret_str)
    return ret_arr

def get_acc(gt_label, pred_cla):
    pos_label_idx = np.where(gt_label == 1)[0]
    neg_label_idx = np.where(gt_label == 0)[0]
    acc_cla_pos = np.mean(gt_label[pos_label_idx] == pred_cla[pos_label_idx])
    acc_cla_neg = np.mean(gt_label[neg_label_idx] == pred_cla[neg_label_idx])

    true_pos = np.sum(gt_label[pos_label_idx] == pred_cla[pos_label_idx])
    true_neg = np.sum(gt_label[neg_label_idx] == pred_cla[neg_label_idx])
    false_pos = np.sum(gt_label[neg_label_idx] != pred_cla[neg_label_idx])
    false_neg = np.sum(gt_label[pos_label_idx] != pred_cla[pos_label_idx])

    return {
        'acc': np.mean(gt_label == pred_cla),
        'acc_cla_pos': 1. * true_pos / (true_pos + false_neg + 1e-6),
        'acc_cla_neg': 1. * true_neg / (true_neg + false_pos + 1e-6),
        'precision': 1. * true_pos / (true_pos + false_pos + 1e-6),
        'recall': 1. * true_pos / (true_pos + false_neg + 1e-6),
    }

import sys
import os
import argparse

BASE_DIR = os.path.dirname(os.path.abspath(__file__))
UTILS_DIR = os.path.abspath(os.path.join(BASE_DIR, '..', 'utils'))
sys.path.append(UTILS_DIR)
from data_helper import *

class LossTracker:
    def __init__(self):
        self.loss_dict = {}

    def add(self, loss_val, loss_name):
        if not loss_name in self.loss_dict:
            self.loss_dict[loss_name] = []

        self.loss_dict[loss_name].append(loss_val)

    def add_dict(self, loss_dict):
        for loss_name in loss_dict:
            self.add(loss_dict[loss_name], loss_name)

    def stat(self):
        mean_dict = {}
        for loss_name in self.loss_dict:
            mean_dict[loss_name] = np.mean(self.loss_dict[loss_name])
        return mean_dict

    def reset(self):
        self.loss_dict = {}

def split_by_key(result_file_name_arr, key='object_cat'):
    ret_dict = {}
    assert key in ['object_cat', 'hook_name', 'object_name']
    for result_file_name in result_file_name_arr:
        if key == 'object_cat':
            _, _, tmp_key, _ = decode_result_file_name(result_file_name)
        elif key == 'hook_name':
            tmp_key, _ = split_result_file_name(result_file_name)
        elif key == 'object_name':
            _, tmp_key = split_result_file_name(result_file_name)
        if not (tmp_key in ret_dict):
            ret_dict[tmp_key] = []
        ret_dict[tmp_key].append(result_file_name)

    ret_dict_len = {tmp:len(ret_dict[tmp]) for tmp in ret_dict}
    return ret_dict, ret_dict_len
class LossTrackerMore():
    def __init__(self):
        self.loss_dict = {}
        self.loss_dict_result_file_name = {}

    def add(self, loss_val, loss_name, result_file_name):
        if not loss_name in self.loss_dict:
            self.loss_dict[loss_name] = []
        if not loss_name in self.loss_dict_result_file_name:
            self.loss_dict_result_file_name[loss_name] = []
        self.loss_dict[loss_name].append(loss_val)
        self.loss_dict_result_file_name[loss_name].append(result_file_name)

    def add_dict(self, loss_dict, result_file_name):
        for loss_name in loss_dict:
            self.add(loss_dict[loss_name], loss_name, result_file_name)

    def per_object_cat(self):
        loss_dict_ret = {}
        for loss_name in self.loss_dict:
            object_name_arr = [split_result_file_name(tmp)[1] for tmp in self.loss_dict_result_file_name[loss_name]]
            object_cat_arr = [split_name(tmp)[0] for tmp in object_name_arr]

            loss_dict_per_obj_cat = {}
            for i, object_cat in enumerate(object_cat_arr):
                if not object_cat in loss_dict_per_obj_cat:
                    loss_dict_per_obj_cat[object_cat] = []
                loss_dict_per_obj_cat[object_cat].append(self.loss_dict[loss_name][i])
            loss_dict_ret[loss_name] = loss_dict_per_obj_cat
        return loss_dict_ret
    def get_result_file_name_per_object_cat(self):
        all_result_file_names = self.loss_dict_result_file_name[list(self.loss_dict.keys())[0]]
        return split_by_key(all_result_file_names, key='object_cat')[0]

    def get_failed_result_file_names(self):
        failed_dict = {}
        for loss_name in self.loss_dict_result_file_name:
            idx = np.where(np.array(self.loss_dict[loss_name]) == 0)[0]
            failed_result_file_name = [self.loss_dict_result_file_name[loss_name][tmp] for tmp in idx]
            failed_dict[loss_name] = {
                'idx': idx,
                'result_file_name': failed_result_file_name
            }
        return failed_dict
    def stat(self):
        mean_dict = {}
        for loss_name in self.loss_dict:
            print(loss_name, np.sum(self.loss_dict[loss_name]), len(self.loss_dict[loss_name]))
            mean_dict[loss_name] = np.mean(self.loss_dict[loss_name])

        loss_dict_obj_cat = self.per_object_cat()
        mean_dict_obj_cat = {}
        for loss_name in loss_dict_obj_cat:
            mean_dict_obj_cat[loss_name] = {
                tmp:np.mean(loss_dict_obj_cat[loss_name][tmp]) for tmp in loss_dict_obj_cat[loss_name]
            }
        return mean_dict, mean_dict_obj_cat

    def print(self):
        mean_dict, mean_dict_obj_cat = self.stat()
        result_file_name_per_object_cat = self.get_result_file_name_per_object_cat()

        print(loss_dict_to_str(mean_dict))

        for tmp in loss_dict_obj_cat_to_str(mean_dict_obj_cat, result_file_name_per_object_cat):
            print(tmp)


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Aug 12 14:39:35 2020

@author: trucdo
"""

#%% Define path, file names, and other parameters

# indicate name (exclude file extension) and path of fasta-formatted .txt file containing TFBS
# output will save to a genome_file_name_PWMmodel subfolder in this path
in_path = "path_to_directory"
tfbs_file_name = "Paeruginosa_fur_sites"

# indicate name of Genbank file (exclude file extension) of organism to scan for TFBS
genome_file_name = "name_of_Genbank_file_without_extension"

# indicate background frequency of each nucleotide for the organism under investigation
nucleotide_prob_dict = {"A": 0.165, "G": 0.335, "C": 0.335, "T": 0.165}

# indicate percent cutoff of top hits
percent_cutoff = 5

#%% Import modules

# to interact with OS-dependent functionality
import os

# to store transcription factor binding site as a Pandas dataframe and to work with arrays
import pandas as pd
#import numpy as np

# to handle Genbank and FASTA files with SeqIO subpackage from Bio package
from Bio import SeqIO

# for some mathematical calculations
import math

# to graph sequence logo
import logomaker as lm

# import matplotlib package for plotting data
import matplotlib.pyplot as plt

# to calculate and record process time
import time

# to handle regex search
import re

#%% Function definitions

def fasta_to_df(in_path, fasta_file_name):
# function to convert known TFBSs formatted as a fasta file into a pandas dataframe
# inputs: path to fasta file (str); name of fasta file (str)
# output: tuple of number of TFBSs and dataframe summarizing TFBSs
# in the dataframe, each row is a different TFBS and each column is a nucleotide within the site

    # navigate to folder containing fasta-formatted file of TFBSs
    os.chdir(in_path)

    # store each binding site as a dict of dict
    # key is name of site and value is an inner dict
    binding_site_dict = {}

    # open and parse the input FASTA file
    fasta_input_file = SeqIO.parse(fasta_file_name + ".txt", "fasta")

    # step through each binding site
    # user counter to keep track of number of TFBS

    no_tfbs = 0

    for binding_site in fasta_input_file:

        # inner dict in which key is the nucleotide position in the sequence and value is the nucleotide
        sequence_dict = {}

        # retrieve position and identity of the nucleotide at the position (i.e. parse the seq into a dict)
        # turn all nucleotides to uppercase
        for position in range(len(binding_site.seq)):
            sequence_dict[position+1] = binding_site.seq[position].upper()

        # assign parsed seq to dict with the key as TFBS name
        binding_site_dict[binding_site.description] = sequence_dict

        no_tfbs = no_tfbs + 1

    # close input fasta file handle
    # comment out beecause Della gives an error saying that fasta_input_file has no attribute close
    #fasta_input_file.close()

    # convert dict to dataframe and transpose so each binding site is a row and each column is a nucleotide in the seq
    binding_site_dataframe = pd.DataFrame(binding_site_dict).transpose()

    return (no_tfbs, binding_site_dataframe)

def df_to_pfm(tfbs_df):
# function to convert pandas dataframe summmarizing TFBSs to a position frequency matrix
# input: dataframe (output from fasta_to_df function)
# output: dataframe that summarizes the number and type of nucleotides at each position across all TFBSs

    # dict to store frequency of each nucleotide at each position across all TFBSs
    # key is position in the TFBS and value is inner dict
    frequency_dict = {}

    # iterate through the columns (i.e. position in seq) of the dataframe
    for position in tfbs_df.columns:

        # extract all nucleotides across all TFBS at current specific position in loop
        # then concatenate nucleotides in the list into a string
        nucl_list = tfbs_df[position].to_list()
        nucl_str = "".join(nucl_list)

        # inner dict wherein key is nucleotide and value is frequency of the nucleotide across all TFBS at a position
        nucleotide_dict = {}
        nucleotide_dict["A"] = nucl_str.count("A")
        nucleotide_dict["C"] = nucl_str.count("C")
        nucleotide_dict["G"] = nucl_str.count("G")
        nucleotide_dict["T"] = nucl_str.count("T")
        nucleotide_dict["total"] = len(nucl_str)

        frequency_dict[position] = nucleotide_dict

    # convert dict to dataframe
    pfm_df = pd.DataFrame(frequency_dict)

    return pfm_df

def pfm_to_pwm(tfbs_pfm, no_tfbs, nucl_prob_dict):
# function to convert position frequency matrix of TFBS to position weight matrix
# inputs: dataframe (index 1 output from df_to_pfm function); number of TFBSs (int; index 0 output of df_to_pfm)
# inputs: dict of background probability for each nucleotide (key is str; value is int)
# output: tuple in which
# index 0: dataframe that summarizes the log-scaled normalized frequency value of each nucleotide at each TFBS position
# index 1: information content (float) of calcluated position weight matrix

    # dict to store log-scaled pseudocount-corrected frequency at each position for all TFBSs
    # key is position in the TFBS and value is inner dict
    weight_dict = {}

    # to store information content of PWM
    info_content = 0

    # iterate through the columns (i.e. position in seq) of the dataframe
    for position in tfbs_pfm.columns:

        # extract the data series corresponding to the current column
        column_series = tfbs_pfm[position]

        # extract the total number of nucleotides at current position (equals no_tfbs in most cases)
        no_nucl = column_series["total"]

        # calculate pseudocount as square root of number of TFBSs used to build the model divided by 4 for the 4 nucl
        pseudocount = math.sqrt(no_tfbs) / 4

        # inner dict wherein key is nucleotide and value is weighted frequency of the nucleotide at current position
        nucleotide_dict = {}

        # step through each nucleotide
        for nucleotide in ["A", "C", "G", "T"]:

            # extract the raw count of current nucleotide in the current column
            nucleotide_count = column_series.loc[nucleotide]

            # calculate corrected probability for current nucleotide in current column
            cor_prob = (nucleotide_count + pseudocount) / (no_nucl + (4 * pseudocount))

            # calculate position weight matrix value of current nucleotide in current position
            # divide corrected nucleotide probability by expected background probability of current nucleotide
            # convert to log base 2
            nucleotide_dict[nucleotide] = math.log(cor_prob/nucl_prob_dict[nucleotide], 2)

            # update information content
            info_content = cor_prob * nucleotide_dict[nucleotide]

        weight_dict[position] = nucleotide_dict

    # convert dict to dataframe
    pwm_df = pd.DataFrame(weight_dict)

    return (pwm_df, info_content)

def find_consensus_seq(tfbs_pwm):
# function to define the highest-scoring sequence (i.e. consensus sequence) and max score from a position weight matrix
# input: position weight matrix dataframe (index 0 output of pfm_to_pwm function)
# output: tuple of consensus sequence (index 0) and its score (index 1)

    # initiate variables to store consensus sequence (str) and score (float)
    consensus_score = 0
    consensus_seq = ""

    # iterate through the columns (i.e. position in seq) of the dataframe
    for position in tfbs_pwm.columns:

        # extract the highest nucleotide PWM score at current position
        consensus_score = consensus_score + max(tfbs_pwm[position])

        # create a bool filter to identity best nucleotide at current position
        max_filter = tfbs_pwm[position] == max(tfbs_pwm[position])

        # apply bool filter to current column to extract identity of highest-scoring nucleotide
        # index method produces pandas index object, so convert to list object
        max_nucl_list = list(tfbs_pwm[position][max_filter].index)

        # in case there are multiple nucleotides that are equally high-scoring at current position
        # store both equally-scoring nucleotides separated by a pipe character
        # and put parentheses around the possible alternate nucleotides at this position
        # if max_nucl_list has only 1 item (i.e. only one high-scoring nucleotide), then no pipe character is added
        if len(max_nucl_list ) > 1:
            max_nucl = "(" + "|".join(max_nucl_list) + ")"
        else:
            max_nucl = "".join(max_nucl_list)

        consensus_seq = consensus_seq + max_nucl

    return (consensus_seq, consensus_score)

def seq_score(tfbs_pwm, dna_seq):
# function to score a DNA sequence based on its degree of match to a TFBS described by a position weight matrix
# inputs: position weight matrix dataframe (output of pfm_to_pwm function); DNA sequence of interest (str)
# output: dataframe summarizing weighted score for each nucleotide of the DNA sequence
# each column is a nucleotide and there are two rows
# the last column is the total score (sum of all nucleotide PWM scores), complete DNA and TFBS sequence
# the first row is the scores and the second row is the nucleotide/DNAsequence

    # check that the DNA sequence has the same length as the TFBSs and return error if different lengths

    if len(dna_seq) != len(tfbs_pwm.columns):
        print("seq_score function failed: DNA sequence length differs from TFBS length")
        return "seq_score function failed: DNA sequence length differs from TFBS length"
    else:

        # ensure DNA query sequence is all uppercase
        DNA_seq = dna_seq.upper()

        # dict to store PWM score of each nucleotide in the DNA sequence for all positions
        # key is position in the TFBS and value is inner dict of nucleotide scores/identity
        score_dict = {}

        # iterate through each nucleotide of the DNA sequence
        # keep track of the nucleotide position to query the correct PWM column
        column_no = 1
        # keep track of sum score of all nucleotides
        sum_score = 0;

        for nucleotide in DNA_seq:

            # inner dict wherein key is nucleotide "score" and "identity" and values follow
            inner_score_dict = {}

            # extract PWM value
            inner_score_dict["score"] = tfbs_pwm.loc[nucleotide,column_no]

            # extract current nucleotide
            inner_score_dict["identity"] = nucleotide

            # store inner dict in outer score_dict
            score_dict[column_no] = inner_score_dict

            # update sum_score
            sum_score = sum_score + inner_score_dict["score"]

            # increase column count by 1 to shift to next column
            column_no = column_no + 1

        # convert score dict to dataframe
        score_df = pd.DataFrame(score_dict, index = ["score", "identity"])

        # add a new final column to store sum of all nucleotide scores and complete DNA sequence
        score_df["sum"] = [sum_score, DNA_seq]

        return score_df

def gbk_to_fasta(in_path, genome_file_name):
# function to convert whole-genome Genbank seqence file for query organism to a FASTA file
# calculates both the forward and reverse-complement of genome sequence
# inputs: path to genome file (str); name of genome file (str)
# output: returns dict of seq_record objects; new fasta file is saved to in_path
# dict keys are names of chromosomes (e.g. chr_1)
# dict values are the seq_record object corresponding to the chromosome

    # navigate to folder containing genome sequence file
    os.chdir(in_path)

    # open Genbank file to read
    input_handle  = open(genome_file_name + ".gbk", "r")
    output_handle = open(genome_file_name + ".fasta", "w")

    # extract DNA sequence of all chromosomes and store in FASTA format
    # for each chromosome, also write out the reverse complement sequence
    # each FASTA entry is a different chromosome
    # store each seq_record (forward chromosome) as a key-value pair in dict

    seq_record_dict = {}

    for seq_record in SeqIO.parse(input_handle, "genbank"):

        # store seq_record object in dict
        seq_record_dict["chr_" + str(seq_record.name)] = seq_record

        # sequence as in Genbank record
        forward_seq = str(seq_record.seq)

        # reverse-complement of sequence
        reverse_seq = forward_seq.replace("A", "t")
        reverse_seq = reverse_seq.replace("C", "g")
        reverse_seq = reverse_seq.replace("G", "c")
        reverse_seq = reverse_seq.replace("T", "a")
        reverse_seq = reverse_seq.upper()
        reverse_seq = reverse_seq[::-1]

        output_handle.write(">chr_" + str(seq_record.name) + "\n" + forward_seq+ "\n")
        output_handle.write(">chr_" + str(seq_record.name) + "_reverse-complement\n" + reverse_seq+ "\n")

    output_handle.close()
    input_handle.close()

    return seq_record_dict

def scan_genome(in_path, genome_fasta_file, tfbs_pwm):
# function to scan whole-genome FASTA sequence file for candidate TFBSs using position weight matrix
# inputs: path to genome file (str); name of FASTA-formatted genome file (str; file output of gbk_to_fasta)
# inputs: position weight matrix dataframe (output of pfm_to_pwm function)
# output: dataframe summarizing sequence and genome position of all sequences with length equal to TFBS length
# dataframe also contains match score calculated for each sequence using PWM

    # navigate to folder containing genome sequence file
    os.chdir(in_path)

    # open FASTA file to read
    input_handle  = open(genome_fasta_file + ".fasta", "r")

    # length of TFBS
    tfbs_len = len(tfbs_pwm.columns)

    # outer dict to store extracted sequence
    # key is arbitrary number keeping track of order when sequence was extracted
    scan_dict = {}
    count = 1

    # scan each chromosome and extract sequence of the same length as TFBS
    for chromosome in SeqIO.parse(input_handle, "fasta"):

        # tack on the first X nucleotides of the chromosome to the end of the chromosome to account for circular chromosome
        # where X is the length of TFBS minus 1
        cir_chr = str(chromosome.seq).rstrip("\n") + str(chromosome.seq[0:(tfbs_len-1)]).rstrip("\n")
        cir_chr = cir_chr.upper()

        # use a sliding window of 1 base and extract sequence and its location on the chromosome
        # scan to the end of the chromosome
        for i in range(0, len(cir_chr) - tfbs_len + 1):

            # inner dict wherein keys: sequence, chromosome_no, location
            # all values are strings except for "location" key, whose value is an int
            inner_dict = {}

            inner_dict["sequence"] = cir_chr[i:i+tfbs_len]
            inner_dict["chromosome_no"] = str(chromosome.name)
            inner_dict["start_location"] = i+1

            # calculate match between the extracted sequence and TFBS using PWM
            # seq_score output is a dataframe of nucleotide score and identity
            # just extract the total PWM score for entire query sequence
            inner_dict["PWM_score"] = seq_score(tfbs_pwm, inner_dict["sequence"]).loc["score","sum"]

            # store inner_dict for current sequence
            scan_dict[count] = inner_dict

            # update count before proceeding to next sequence
            count = count + 1

    # convert dict to dataframe and transpose so each extracted sequence is a row and its info is a column
    extracted_fragments = pd.DataFrame(scan_dict).transpose()

    return extracted_fragments

def filter_cds(seq_record_dict, fragments_df):
# function to filter TFBS matches for those entirely within non-coding (i.e. regulatory) regions
# inputs: dict of SeqRecord objects (return output of gbk_to_fasta function)
# inputs: dataframe of fragments to assess (output of scan_genome before or after PWM filtering)
# output: filtered dataframe of all sequence fragments that are fully wthin non-coding regions
# dataframe also contains genome position and match score calculated for each sequence using PWM

    # iterate through all key-value pairs (i.e. all chromosomes and their SeqRecord objects)
    # add range of coordinates of all CDS to a set stored in a dict under chromosome name key
    # note that all coordinates correspond to forward-strand position
    cds_coor_dict = {}
    for chr_name, seq_record in seq_record_dict.items():
        # set to store coordinates of all CDS for current chromosome
        cds_coor_set = set()
        for seq_feature in seq_record.features:
            if seq_feature.type=="CDS":
                # account for the start and end location numbering following Python's zero-index scheme
                fstart = int(seq_feature.location.start) + 1
                fend = int(seq_feature.location.end)
                # add coordinates within the CDS range to set of all CDS coordinates in current chromosome
                cds_coor_set.update(list(range(fstart, fend + 1)))
        # store set of CDS coordinates for current chromosome in dict
        cds_coor_dict[chr_name] = cds_coor_set

    # iterate through all the sequence fragments (each row) stored in dataframe
    # check if the entire fragment is present within non-coding regions of the genome
    # add fully non-coding fragments to filtered dict

    noncds_dict = {}

    for index, row in fragments_df.iterrows():

        # retrieve name of current chromosome and set rc_flag to default False for forward strand
        chr_no = row["chromosome_no"]
        rc_flag = False

        # if current fragment is found on reverse strand
        # remove the reverse complement from chromosome name for later search
        # flag to convert reverse-complement coordinates to forward coordinates
        if re.search("_reverse-complement", chr_no):
            chr_no = chr_no.rstrip("_reverse-complement")
            rc_flag = True

        # retrieve length of current chromosome
        chr_len = len(seq_record_dict[chr_no].seq)

        # for forward strand
        # retrieve coordinates as is of all nucleotides that comprise the current sequence fragment as a list
        if rc_flag == False:
            start = row["start_location"]
            interval = len(row["sequence"])
            all_coor_list = list(range(start, start + interval))
        # for reverse-complement strand
        else:
            start = chr_len - row["start_location"] + 1
            interval = len(row["sequence"])
            all_coor_list = list(range(start, start - interval, -1))

        # iterate through all coordinates in list and check if nucleotide position is within CDS set for current chr
        # default is fragment fully within non-coding region
        # break search to exit loop if any nucleotide is within CDS set
        cds_flag = False

        # to account for extra bases at the end due to wrapping around to position 1 for circular chromosome
        # first convert nucleotides at the end to actual coordinates
        # the last X positions of a chromosome is equal to the first X positions of the same chr
        # where X = TFBS length - 1

        for coor in all_coor_list:
            if coor > chr_len:
                coor = coor - chr_len
            if coor in cds_coor_dict[chr_no]:
                cds_flag = True
                break

        # if entire sequence fragment is outside CDS, add to dict
        if cds_flag == False:
            noncds_dict[index] = row

    # convert dict to dataframe and transpose so each filtered sequence is a row and its info is a column
    noncds_fragments = pd.DataFrame(noncds_dict).transpose()

    return noncds_fragments

def find_cds_neighbor(seq_record_dict, fragments_df):
# function to find the adjacent CDS to the candidate TFBS
# inputs: dict of SeqRecord objects (return output of gbk_to_fasta function)
# inputs: dataframe of fragments to assess (output of scan_genome before or after filtering)
# output: dataframe of input sequence fragments with new column describing adjacent CDS

    # iterate through all key-value pairs (i.e. all chromosomes and their SeqRecord objects)
    # note that all coordinates correspond to forward-strand position

    # outer dict key is chromosome name and value is multi-layered dict
    chr_dict = {}

    for chr_name, seq_record in seq_record_dict.items():

        # middle dict separating start/end coordinates of forward and reverse-complement strands
        strand_dict = {}

        # in-most dicts wherein keys are start/end coordinates and values are tuple (CDS locus tags, annotated product)
        # separate start/end and forward/reverse strands into different dicts
        scoor_dict = {}
        ecoor_dict = {}
        rscoor_dict = {}
        recoor_dict = {}

        for seq_feature in seq_record.features:

            if seq_feature.type=="CDS":

                # extract start/end positions of CDS
                # account for Python vs. Genbank indexing by adding 1 to start position
                start_pos = int(seq_feature.location.start) + 1
                end_pos = int(seq_feature.location.end)

                # for forward CDS sequence
                if seq_feature.location.strand == 1:
                    scoor_dict[start_pos] = (seq_feature.qualifiers["locus_tag"][0], seq_feature.qualifiers["product"][0])
                    strand_dict["forward_start"] = scoor_dict
                    ecoor_dict[end_pos] = (seq_feature.qualifiers["locus_tag"][0], seq_feature.qualifiers["product"][0])
                    strand_dict["forward_end"] = ecoor_dict
                # for reverse-complement CDS sequence
                else:
                    rscoor_dict[start_pos] = (seq_feature.qualifiers["locus_tag"][0], seq_feature.qualifiers["product"][0])
                    strand_dict["reverse_start"] = rscoor_dict
                    recoor_dict[end_pos] = (seq_feature.qualifiers["locus_tag"][0], seq_feature.qualifiers["product"][0])
                    strand_dict["reverse_end"] = recoor_dict

        # store coordinates of current chromosome in multi-layered dict
        chr_dict[chr_name] = strand_dict

    # iterate through all candidate TFBS fragments and search for closest CDS on same strand
    cds_info_dict = {}

    for index, row in fragments_df.iterrows():

        # retrieve name and start position of current chromosome
        # set rc_flag to default False for forward strand
        chr_no = row["chromosome_no"]
        rc_flag = False
        fragment_loc = row["start_location"]

        # if current fragment is found on reverse strand
        # remove the reverse complement from chromosome name for later search
        # flag to convert reverse-complement coordinates to forward coordinates
        # convert start position on reverse-complement strand to corresponding forward position strand
        if re.search("_reverse-complement", chr_no):
            chr_no = chr_no.rstrip("_reverse-complement")
            rc_flag = True
            fragment_loc = len(seq_record_dict[chr_no].seq) - row["start_location"] + 1

        # to search on forward strand
        if rc_flag == False:
            strand = chr_dict[chr_no]
            coor = strand["forward_start"]
            all_coor = list(coor.keys())
            # add start position of current TFBS fragment to list and change to ascending order
            all_coor.append(fragment_loc)
            all_coor.sort()
            # find index of TFBS fragment in list and retrieve number in index + 1 position (i.e. closest forward CDS)
            closest_CDS_index = all_coor.index(fragment_loc) + 1
            # if TFBS is at the end of the list, then loop back to the first item in list
            if closest_CDS_index > len(all_coor) - 1:
                closest_CDS_index = 0
            # use retrieved position to obtain CDS locus tag and annnotation from multi-layered dict
            closest_CDS_pos = all_coor[closest_CDS_index]
            cds_name = coor[closest_CDS_pos][0]
            cds_annotation = coor[closest_CDS_pos][1]
        # to search on reverse strand
        else:
            strand = chr_dict[chr_no]
            coor = strand["reverse_end"]
            all_coor = list(coor.keys())
            # add start position of current TFBS fragment to list and change to ascending order
            all_coor.append(fragment_loc)
            all_coor.sort()
            # find index of TFBS fragment in list and retrieve number in index - 1 position (i.e. closest reverse CDS)
            # if TFBS is at the beginnning of the list, then loop forward to the last item in list (index of -1)
            closest_CDS_index = all_coor.index(fragment_loc) - 1
            # use retrieved position to obtain CDS locus tag from multi-layered dict
            closest_CDS_pos = all_coor[closest_CDS_index]
            cds_name = coor[closest_CDS_pos][0]
            cds_annotation = coor[closest_CDS_pos][1]

        # convert current row (i.e. current sequence fragment info) to dict
        # add two new key-value pairs for CDS_neighbor and CDS_annotation
        temp_dict = dict(row)
        temp_dict["CDS_neighbor"] = cds_name
        temp_dict["CDS_annotation"] = cds_annotation
        cds_info_dict[index] = temp_dict

    # convert dict to dataframe and transpose so each filtered sequence is a row and its info is a column
    cds_info = pd.DataFrame(cds_info_dict).transpose()

    return cds_info

#%% Main program

# record and print the process start time in second
start_time = time.time()
print("Start time: " + str(start_time) + " seconds")

# convert TFBSs into dataframe
no_tfbs, binding_site_df = fasta_to_df(in_path, tfbs_file_name)
# convert TFBS dataframe into position frequency matrix
tfbs_pfm = df_to_pfm(binding_site_df)
# convert position frequency matrix to position-specific scoring matrix and information content of PWM
tfbs_pwm, info_content = pfm_to_pwm(tfbs_pfm, no_tfbs, nucleotide_prob_dict)
# determine consensus sequence and max score based on position frequency matrix of input TFBSs
consensus_seq, consensus_score = find_consensus_seq(tfbs_pwm)

# convert whole-genome Genbank sequence file of query organism to FASTA format and save dict of SeqRecord objects
seq_record_dict = gbk_to_fasta(in_path,genome_file_name)

# scan entire genome and extract sequence fragments with length equal to TFBS length
genome_scan_df = scan_genome(in_path, genome_file_name, tfbs_pwm)

# filter genome for only those sequences with PWM scores > information content
ICfiltered_genome_df = genome_scan_df[genome_scan_df["PWM_score"] > info_content]

# filter IC-filtered genome further for TFBS matches within non-coding regions
CDS_IC_filtered_genome_df = filter_cds(seq_record_dict, ICfiltered_genome_df)

# also filter entire genome for TFBS matches within non-coding regions
CDS_filtered_genome_df = filter_cds(seq_record_dict, genome_scan_df)

# annotate CDS-IC filtered TFBS matches to find closest CDS
CDS_IC_neighbor_df = find_cds_neighbor(seq_record_dict, CDS_IC_filtered_genome_df)

# also annotate CDS filtered TFBS matches to find closest CDS
CDS_neighbor_df = find_cds_neighbor(seq_record_dict, CDS_filtered_genome_df)

# filter genome for only those sequences with PWM scores > user-defined cutoff
index_cutoff = math.floor(genome_scan_df.shape[0] * (100-percent_cutoff)/100)
ordered_df = genome_scan_df.sort_values("PWM_score", axis = 0, ascending = True)
ufiltered_genome_df = ordered_df.iloc[index_cutoff:,:]

# record and print the process end time in seconds
end_time = time.time()
print("End time: " + str(end_time) + " seconds")

# record and print the process time in min and second
elapsed_time = end_time-start_time
# if the process took less than a minute
if elapsed_time < 60:
    elapsed_min = 0
    elapsed_sec = elapsed_time
# if the process took more than a minute
else:
    elapsed_min = math.floor(elapsed_time/60)
    elapsed_sec = (elapsed_time/60 - elapsed_min) * 60

print("Elapsed time: " + str(elapsed_min) + " minutes and " + str(elapsed_sec) + " seconds")

#%% Data visualization

# make and navigate to new directory to store results
out_path = in_path + "/" + genome_file_name + "_PWMmodel/"
if os.path.exists(out_path):
    os.chdir(out_path)
else:
    os.mkdir(out_path)
    os.chdir(out_path)

# plot the sequence logo using the position frequency matrix of input TFBSs
# also save sequence logo to out_path
logo = lm.Logo(tfbs_pwm.transpose(), font_name = 'Arial Rounded MT Bold', fade_below=0.5, shade_below=0.5)
logo.ax.set_xlabel('Position',fontsize=14)
logo.ax.set_ylabel("Position-specific score", labelpad=-1,fontsize=14)
plt.savefig(tfbs_file_name + "_sequencelogo.png", bbox_inches = "tight")

# plot histogram to visualize distribution of PWM scores for all sequences
# save histogram to out_path
fig = plt.figure()
ax = fig.gca()
plt.title("All PWM scores from genome scan", fontsize=18.0)
plt.xlabel("PWM score", fontsize=18.0)
plt.ylabel("Number of sequences", fontsize=18.0)
plt.hist(genome_scan_df["PWM_score"], bins=20, alpha=0.6, edgecolor="black", linewidth=1.0)
# plot vertical line indicating threshold = information content and user-defined cutoff
plt.axvline(x = info_content, ymin = 0, color='seagreen', linewidth = 2)
plt.axvline(x = ufiltered_genome_df.iloc[0]["PWM_score"], ymin = 0, color='orange', linewidth = 2)
plt.xlim([-consensus_score,consensus_score])
plt.grid(False)
plt.tick_params(axis="both", direction="out", length=5.0, width=2.5, color="gray", labelsize=18.0)
ax.spines["bottom"].set_linewidth(1.5)
ax.spines["top"].set_linewidth(1.5)
ax.spines["left"].set_linewidth(1.5)
ax.spines["right"].set_linewidth(1.5)
plt.savefig(genome_file_name + "_predictedtfbs.png", bbox_inches = "tight")
plt.show()

#%% Write out results

tfbs_pfm.to_csv(tfbs_file_name + "_pfm.csv")
tfbs_pwm.to_csv(tfbs_file_name + "_pwm.csv")
genome_scan_df.to_csv(genome_file_name + "_allfragments.csv")
ICfiltered_genome_df.to_csv(genome_file_name + "_ICfilteredfragments.csv")
ufiltered_genome_df.to_csv(genome_file_name + "_Ufilteredfragments.csv")
CDS_IC_neighbor_df.to_csv(genome_file_name + "_CDS-IC-filteredfragments.csv")
CDS_neighbor_df.to_csv(genome_file_name + "_CDSfilteredfragments.csv")

out_file = open(genome_file_name + "_stats.txt", "w")
out_file.write("----------PWMmodel.py----------\n\n")
out_file.write("Input file(s):" + tfbs_file_name + ", " + genome_file_name + "\n")
out_file.write("Output file(s): " + tfbs_file_name + "_pfm.csv, " + tfbs_file_name + "_pwm.csv, ")
out_file.write(genome_file_name + "_allfragments.csv, " + genome_file_name + "_ICfilteredfragments.csv, ")
out_file.write(genome_file_name + "_Ufilteredfragments.csv, " + genome_file_name + "_CDS-IC-filteredfragments.csv, ")
out_file.write(genome_file_name + "_CDSfilteredfragments.csv, " + genome_file_name + "_stats.txt, " + genome_file_name + ".fasta, ")
out_file.write(tfbs_file_name + "_sequencelogo.png, " + genome_file_name + "_predictedtfbs.png\n")
out_file.write("Number of TFBSs used to build model: " + str(no_tfbs) + "\n")
out_file.write("Length of TFBSs: " + str(binding_site_df.shape[1]) + " bases\n")
out_file.write("Background probability of nucleotides: " + str(nucleotide_prob_dict) + "\n")
out_file.write("TFBS Consensus sequence: " + consensus_seq + "\n")
out_file.write("PWM score for perfect match to consensus sequence: " + str(consensus_score) + "\n")
out_file.write("Information content (IC): " + str(info_content) + "\n")
out_file.write("User-defined percent cutoff for top hits: " + str(percent_cutoff) + "\n")
out_file.write("Number of genome fragments assessed: " + str(genome_scan_df.shape[0]) + "\n")
out_file.write("Number of IC-filtered genome fragments: " + str(ICfiltered_genome_df.shape[0]) + "\n")
out_file.write("Number of CDS and IC-filtered genome fragments: " + str(CDS_IC_filtered_genome_df.shape[0]) + "\n")
out_file.write("Number of CDS-filtered genome fragments: " + str(CDS_filtered_genome_df.shape[0]) + "\n")
out_file.write("Number of user-filtered genome fragments " + str(ufiltered_genome_df.shape[0]) + "\n")
out_file.write("Elapsed time: " + str(elapsed_min) + " minutes and " + str(elapsed_sec) + " seconds\n")
out_file.close()

max_stars_repo_name: ajlinklab/PupB

max_issues_repo_name: ajlinklab/PupB

max_forks_repo_name: ajlinklab/PupB


text:
import cv2
import numpy as np
import os

path='./src/'
folders=os.listdir(path)
print(folders)
for folder in folders:
    folder_path=os.path.join(path,folder)
    contents=os.listdir(folder_path)
    save_path=os.path.join("./dataset/",folder)
    if not os.path.isdir(save_path):
        os.mkdir(save_path)
    for content in contents:
        print(content)
        image=cv2.imread(os.path.join(folder_path,content))
        cv2.imwrite(os.path.join(save_path, content),image)
        # cv2.waitKey(0)
        # print(save_path)
        print(content)
        print(content.find("."))
        fname=content[0:content.find(".")]
        print(fname)

        height,width,channel=image.shape
        start_row, start_col = int(0), int(0)
        end_row, end_col = int(height * .5), int(width)
        cropped_top = image[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'0.jpg'),cropped_top)

        image2=cv2.imread(os.path.join(save_path,fname+'0.jpg'))
        height,width,channel=image.shape
        start_row, start_col = int(0), int(0)
        end_row, end_col = int(height), int(width*.5)
        cropped_left = image2[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'2.jpg'),cropped_left)

        start_row, start_col = int(0), int(width*.5)
        end_row, end_col = int(height), int(width)
        cropped_right = image2[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'3.jpg'),cropped_right)

        height,width,channel=image.shape
        start_row, start_col = int(height * .5), int(0)
        end_row, end_col = int(height), int(width)
        cropped_bot = image[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'1.jpg'),cropped_bot)

        image2=cv2.imread(os.path.join(save_path,fname+'1.jpg'))
        height,width,channel=image.shape
        start_row, start_col = int(0), int(0)
        end_row, end_col = int(height), int(width*.5)
        cropped_left = image2[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'4.jpg'),cropped_left)

        start_row, start_col = int(0), int(width*.5)
        end_row, end_col = int(height), int(width)
        cropped_right = image2[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'5.jpg'),cropped_right)

        height,width,channel=image.shape
        start_row, start_col = int(0), int(0)
        end_row, end_col = int(height), int(width*.5)
        cropped_left = image[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'6.jpg'),cropped_left)

        start_row, start_col = int(0), int(width*.5)
        end_row, end_col = int(height), int(width)
        cropped_right = image[start_row:end_row , start_col:end_col]
        cv2.imwrite(os.path.join(save_path,fname+'7.jpg'),cropped_right)


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

import re
import pickle
import numpy as np
from gensim import models
from keras.preprocessing.text import Tokenizer
from keras.preprocessing.sequence import pad_sequences

# 读取数据
def read_data(file_path):
    with open(file_path, encoding='utf-8') as f:
        data = f.readlines()

    data = [re.split('\t', i) for i in data]
    q1 = [i[1] for i in data]

    q2 = [i[2] for i in data]
    q3 = [i[3] for i in data]

    q4 = [i[4] for i in data]
    label = [int(i[5]) for i in data]
    return q1, q2, q3, q4, label

train_q1, train_q2, train_q3, train_q4, train_label = read_data('../data/BQ_train_char_word.txt')
dev_q1, dev_q2, dev_q3, dev_q4, dev_label = read_data('../data/BQ_dev_char_word.txt')
test_q1, test_q2, test_q3, test_q4, test_label = read_data('../data/BQ_test_char_word.txt')
# 构造训练word2vec的语料库
corpus = train_q1 + train_q2 + train_q3 + train_q4 + test_q1 + test_q2 + test_q3 + test_q4 + dev_q1 + dev_q2 + dev_q3 + dev_q4
w2v_corpus = [i.split() for i in corpus]
#词表
word_set = set(' '.join(corpus).split())

MAX_SEQUENCE_LENGTH = 30  # sequence最大长度为30个词
EMBDIM = 300  # 词向量为300维

# 训练word2vec模型
w2v_model = models.Word2Vec(w2v_corpus, size=EMBDIM, window=5, min_count=1, sg=1, workers=4, seed=1234, iter=25)
w2v_model.save('w2v_model.pkl')

tokenizer = Tokenizer(num_words=len(word_set))
tokenizer.fit_on_texts(corpus)
#生成的词表长度，等于num_words=len(word_set)
L = len(tokenizer.word_index)

train_q1 = tokenizer.texts_to_sequences(train_q1)
train_q2 = tokenizer.texts_to_sequences(train_q2)

train_q3 = tokenizer.texts_to_sequences(train_q3)
train_q4 = tokenizer.texts_to_sequences(train_q4)

test_q1 = tokenizer.texts_to_sequences(test_q1)
test_q2 = tokenizer.texts_to_sequences(test_q2)

test_q3 = tokenizer.texts_to_sequences(test_q3)
test_q4 = tokenizer.texts_to_sequences(test_q4)

dev_q1 = tokenizer.texts_to_sequences(dev_q1)
dev_q2 = tokenizer.texts_to_sequences(dev_q2)

dev_q3 = tokenizer.texts_to_sequences(dev_q3)
dev_q4 = tokenizer.texts_to_sequences(dev_q4)

train_pad_q1 = pad_sequences(train_q1, maxlen=MAX_SEQUENCE_LENGTH)
train_pad_q2 = pad_sequences(train_q2, maxlen=MAX_SEQUENCE_LENGTH)
train_pad_q3 = pad_sequences(train_q3, maxlen=MAX_SEQUENCE_LENGTH)
train_pad_q4 = pad_sequences(train_q4, maxlen=MAX_SEQUENCE_LENGTH)

test_pad_q1 = pad_sequences(test_q1, maxlen=MAX_SEQUENCE_LENGTH)
test_pad_q2 = pad_sequences(test_q2, maxlen=MAX_SEQUENCE_LENGTH)
test_pad_q3 = pad_sequences(test_q3, maxlen=MAX_SEQUENCE_LENGTH)
test_pad_q4 = pad_sequences(test_q4, maxlen=MAX_SEQUENCE_LENGTH)

dev_pad_q1 = pad_sequences(dev_q1, maxlen=MAX_SEQUENCE_LENGTH)
dev_pad_q2 = pad_sequences(dev_q2, maxlen=MAX_SEQUENCE_LENGTH)
dev_pad_q3 = pad_sequences(dev_q3, maxlen=MAX_SEQUENCE_LENGTH)
dev_pad_q4 = pad_sequences(dev_q4, maxlen=MAX_SEQUENCE_LENGTH)

embedding_matrix = np.zeros([len(tokenizer.word_index) + 1, EMBDIM])

for word, idx in tokenizer.word_index.items():
    embedding_matrix[idx, :] = w2v_model.wv[word]

def save_pickle(fileobj, file_path):
    with open(file_path, 'wb') as f:
        pickle.dump(fileobj, f)

def load_pickle(file_path):
    with open(file_path, 'rb') as f:
        fileobj = pickle.load(f)
    return fileobj

model_data = {'train_q1': train_pad_q1, 'train_q2': train_pad_q2, 'train_q3': train_pad_q3, 'train_q4': train_pad_q4,  'train_label': train_label,
              'test_q1': test_pad_q1, 'test_q2': test_pad_q2, 'test_q3': test_pad_q3, 'test_q4': test_pad_q4, 'test_label': test_label,
              'dev_q1': dev_pad_q1, 'dev_q2': dev_pad_q2, 'dev_q3': dev_pad_q3, 'dev_q4': dev_pad_q4, 'dev_label': dev_label}

save_pickle(corpus, 'corpus.pkl')
save_pickle(model_data, 'model_data.pkl')
save_pickle(embedding_matrix, 'embedding_matrix.pkl')
save_pickle(tokenizer, 'tokenizer.pkl')

max_stars_repo_name: ZQpengyu/MAFFo

max_issues_repo_name: ZQpengyu/MAFFo

max_forks_repo_name: ZQpengyu/MAFFo