PFP-WGAN/conditional_new_nn_hirearchical_seq.py

#!/usr/bin/env python

"""
parts of codes are based on the work in https://github.com/bio-ontology-research-group/deepgo. 

"""
from __future__ import division

import logging
import sys
import time
from collections import deque
from multiprocessing import Pool

import click as ck
import numpy as np
import pandas as pd
import tensorflow as tf
from keras import backend as K
from keras.callbacks import EarlyStopping, ModelCheckpoint
from keras.layers import (
    Dense, Input, SpatialDropout1D, Conv1D, MaxPooling1D,
    Flatten, Concatenate, Add, Maximum, Embedding, BatchNormalization, Activation, Dropout)
from keras.losses import binary_crossentropy
from keras.models import Sequential, Model, load_model
from keras.preprocessing import sequence
from scipy.spatial import distance
from sklearn.metrics import log_loss
from sklearn.metrics import roc_curve, auc, matthews_corrcoef
from keras.layers import Lambda
from sklearn.metrics import precision_recall_curve

from utils import (
    get_gene_ontology,
    get_go_set,
    get_anchestors,
    get_parents,
    DataGenerator,
    FUNC_DICT,
    get_height,
    get_ipro)
from conditional_wgan_wrapper_post import WGAN_wrapper, wasserstein_loss, generator_recunstruction_loss_new

config = tf.ConfigProto()
config.gpu_options.allow_growth = True
sess = tf.Session(config=config)
K.set_session(sess)

logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.INFO)
sys.setrecursionlimit(100000)

DATA_ROOT = 'data/swiss/'
MAXLEN = 1000
REPLEN = 256
ind = 0


@ck.command()
@ck.option(
    '--function',
    default='bp',
    help='Ontology id (mf, bp, cc)')
@ck.option(
    '--device',
    default='gpu:0',
    help='GPU or CPU device id')
@ck.option(
    '--org',
    default= None,
    help='Organism id for filtering test set')
@ck.option('--train',default = True, is_flag=True)
@ck.option('--param', default=0, help='Param index 0-7')
def main(function, device, org, train, param):
    global FUNCTION
    FUNCTION = function
    global GO_ID
    GO_ID = FUNC_DICT[FUNCTION]
    global go
    go = get_gene_ontology('go.obo')
    global ORG
    ORG = org
    func_df = pd.read_pickle(DATA_ROOT + FUNCTION + '.pkl')
    global functions
    functions = func_df['functions'].values
    global func_set
    func_set = set(functions)
    global all_functions
    all_functions = get_go_set(go, GO_ID)
    logging.info('Functions: %s %d' % (FUNCTION, len(functions)))
    if ORG is not None:
        logging.info('Organism %s' % ORG)
    global go_indexes
    go_indexes = dict()
    for ind, go_id in enumerate(functions):
        go_indexes[go_id] = ind
    global node_names
    node_names = set()
    with tf.device('/' + device):
        params = {
            'fc_output': 1024,
            'learning_rate': 0.001,
            'embedding_dims': 128,
            'embedding_dropout': 0.2,
            'nb_conv': 1,
            'nb_dense': 1,
            'filter_length': 128,
            'nb_filter': 32,
            'pool_length': 64,
            'stride': 32        
        }
        model(params, is_train=train)

def load_data():

    df = pd.read_pickle(DATA_ROOT + 'train' + '-' + FUNCTION + '.pkl')
    n = len(df)
    index = df.index.values
    valid_n = int(n * 0.8)
    train_df = df.loc[index[:valid_n]]
    valid_df = df.loc[index[valid_n:]]
    test_df = pd.read_pickle(DATA_ROOT + 'test' + '-' + FUNCTION + '.pkl')
    print( test_df['orgs'] )
    if ORG is not None:
        logging.info('Unfiltered test size: %d' % len(test_df))
        test_df = test_df[test_df['orgs'] == ORG]
        logging.info('Filtered test size: %d' % len(test_df))

    # Filter by type
    # org_df = pd.read_pickle('data/prokaryotes.pkl')
    # orgs = org_df['orgs']
    # test_df = test_df[test_df['orgs'].isin(orgs)]

    def reshape(values):
        values = np.hstack(values).reshape(
            len(values), len(values[0]))
        return values

    def normalize_minmax(values):
        mn = np.min(values)
        mx = np.max(values)
        if mx - mn != 0.0:
            return (values - mn) / (mx - mn)
        return values - mn

    def get_values(data_frame):
        print(data_frame['labels'].values.shape)
        labels = reshape(data_frame['labels'].values)
        ngrams = sequence.pad_sequences(
            data_frame['ngrams'].values, maxlen=MAXLEN)
        ngrams = reshape(ngrams)
        rep = reshape(data_frame['embeddings'].values)
        data = ngrams
        return data, labels

    train = get_values(train_df)
    valid = get_values(valid_df)
    test = get_values(test_df)

    return train, valid, test, train_df, valid_df, test_df


def get_feature_model(params):
    embedding_dims = params['embedding_dims']
    max_features = 8001
    model = Sequential()
    model.add(Embedding(
        max_features,
        embedding_dims,
        input_length=MAXLEN))
    model.add(SpatialDropout1D(0.4))
    for i in range(params['nb_conv']):
        model.add(Conv1D(
            activation="relu",
            padding="valid",
            strides=1,
            filters=params['nb_filter'],
            kernel_size=params['filter_length']))
    model.add(MaxPooling1D(strides=params['stride'], pool_size=params['pool_length']))
    model.add(Flatten())
    return model

def merge_outputs(outputs, name):
    if len(outputs) == 1:
        return outputs[0]
    ## return merge(outputs, mode='concat', name=name, concat_axis=1)
    return Concatenate(axis=1, name=name)(outputs)


def merge_nets(nets, name):
    if len(nets) == 1:
        return nets[0]
    ## return merge(nets, mode='sum', name=name)
    return Add(name=name)(nets)

def get_node_name(go_id, unique=False):
    name = go_id.split(':')[1]
    if not unique:
        return name
    if name not in node_names:
        node_names.add(name)
        return name
    i = 1
    while (name + '_' + str(i)) in node_names:
        i += 1
    name = name + '_' + str(i)
    node_names.add(name)
    return name


def get_layers(inputs):
    q = deque()
    layers = {}
    name = get_node_name(GO_ID)
    layers[GO_ID] = {'net': inputs}
    for node_id in go[GO_ID]['children']:
        if node_id in func_set:
            q.append((node_id, inputs))
    while len(q) > 0:
        node_id, net = q.popleft()
        parent_nets = [inputs]
        # for p_id in get_parents(go, node_id):
        #     if p_id in func_set:
        #         parent_nets.append(layers[p_id]['net'])
        # if len(parent_nets) > 1:
        #     name = get_node_name(node_id) + '_parents'
        #     net = merge(
        #         parent_nets, mode='concat', concat_axis=1, name=name)
        name = get_node_name(node_id)
        net, output = get_function_node(name, inputs)
        if node_id not in layers:
            layers[node_id] = {'net': net, 'output': output}
            for n_id in go[node_id]['children']:
                if n_id in func_set and n_id not in layers:
                    ok = True
                    for p_id in get_parents(go, n_id):
                        if p_id in func_set and p_id not in layers:
                            ok = False
                    if ok:
                        q.append((n_id, net))

    for node_id in functions:
        childs = set(go[node_id]['children']).intersection(func_set)
        if len(childs) > 0:
            outputs = [layers[node_id]['output']]
            for ch_id in childs:
                outputs.append(layers[ch_id]['output'])
            name = get_node_name(node_id) + '_max'
            ## layers[node_id]['output'] = merge(
            ##     outputs, mode='max', name=name)
            layers[node_id]['output'] = Maximum(name=name)(outputs)
    return layers

def get_function_node(name, inputs):
    output_name = name + '_out'
    # net = Dense(256, name=name, activation='relu')(inputs)
    output = Dense(1, name=output_name, activation='sigmoid')(inputs)
    return output, output




def get_generator(params, n_classes):
    inputs = Input(shape=(MAXLEN,), dtype='int32', name='input1')
    feature_model = get_feature_model(params)(inputs)
    net = Dense(300, activation='relu')(feature_model)
    net = BatchNormalization()(net)
    layers = get_layers(net)
    output_models = []
    for i in range(len(functions)):
        output_models.append(layers[functions[i]]['output'])
    net = Concatenate(axis=1)(output_models)
    output = Dense(n_classes, activation='sigmoid')(net)
    model = Model(inputs=inputs, outputs=output)

    return model




def get_discriminator(params, n_classes, dropout_rate=0.5):
    inputs = Input(shape=(n_classes, ))
    inputs2 = Input(shape =(MAXLEN,), dtype ='int32', name='d_input2')
    x2 = Embedding(8001,128, input_length=MAXLEN)(inputs2)
    x2 = Conv1D(filters =1 , kernel_size= 1, padding = 'valid', activation ='relu', strides=1)(x2)
    x2 = Lambda(lambda x: K.squeeze(x, 2))(x2)

    #for i in range(params['nb_conv']):
    #    x2 = Conv1D ( activation="relu", padding="valid", strides=1, filters=params['nb_filter'],kernel_size=params['filter_length'])(x2)
    

    #x2 =MaxPooling1D(strides=params['stride'], pool_size=params['pool_length'])(x2)
    #x2 =  Flatten()(x2)
    
    size = 40
    x = inputs
    x = Dropout(dropout_rate)(x)
    x = Dense(size)(x)
    x = BatchNormalization()(x)
    x = Activation('relu')(x)


    size = 40
    x2 = Dropout(dropout_rate)(x2)
    x2 = Dense(size)(x2)
    x2 = BatchNormalization()(x2)
    x2 = Activation('relu')(x2)


    x = Concatenate(axis =1 , name = 'merged2')([x, x2])
    layer_sizes = [80, 40,30]
    for size in layer_sizes:
        x = Dropout(dropout_rate)(x)
        x = Dense(size)(x)
        x = BatchNormalization()(x)
        x = Activation('relu')(x)


    outputs = Dense(1)(x)
    model = Model(inputs = [inputs ,inputs2], outputs=outputs, name='Discriminator')

    return model



def get_model(params,nb_classes, batch_size, GRADIENT_PENALTY_WEIGHT=10):
    generator = get_generator(params, nb_classes)
    discriminator = get_discriminator(params, nb_classes)

    generator_model, discriminator_model = \
        WGAN_wrapper(generator=generator,
                     discriminator=discriminator,
                     generator_input_shape=(MAXLEN,),
                     discriminator_input_shape=(nb_classes,),
                     discriminator_input_shape2 = (MAXLEN, ),
                     batch_size=batch_size,
                     gradient_penalty_weight=GRADIENT_PENALTY_WEIGHT)

    logging.info('Compilation finished')
    return generator_model, discriminator_model


def train_wgan(generator_model, discriminator_model, batch_size, epochs,
               x_train, y_train, x_val, y_val, generator_model_path, discriminator_model_path,
               TRAINING_RATIO=10, N_WARM_UP=0):
    BATCH_SIZE = batch_size
    N_EPOCH = epochs
    
   
    positive_y = np.ones((batch_size, 1), dtype=np.float32)
    zero_y = positive_y * 0
    negative_y = -positive_y
    positive_full_y = np.ones((BATCH_SIZE * TRAINING_RATIO, 1), dtype=np.float32)
    dummy_y = np.zeros((BATCH_SIZE, 1), dtype=np.float32)
    
    positive_full_enable_train = np.ones((len(x_train), 1), dtype = np.float32 )
    positive_full_enable_val = np.ones((len(x_val), 1), dtype =np.float32 )
    #positive_enable_train = np.ones((1, batch_size),dtype = np.float32 )
    #positive_full_train_enable = np.ones((1,BATCH_SIZE * TRAINING_RATIO ), dtype=np.float32 )
    best_validation_loss = None

    for epoch in range(N_EPOCH):
        # np.random.shuffle(X_train)
        print("Epoch: ", epoch)
        print("Number of batches: ", int(y_train.shape[0] // BATCH_SIZE))
        discriminator_loss = []
        generator_loss = []
        minibatches_size = BATCH_SIZE * TRAINING_RATIO

        shuffled_indexes = np.random.permutation(x_train.shape[0])
        shuffled_indexes_2 = np.random.permutation(x_train.shape[0])

        for i in range(int(y_train.shape[0] // (BATCH_SIZE * TRAINING_RATIO))):
            batch_indexes = shuffled_indexes[i * minibatches_size:(i + 1) * minibatches_size]
            batch_indexes_2 = shuffled_indexes_2[i * minibatches_size:(i + 1) * minibatches_size]
            x = x_train[batch_indexes]
            y = y_train[batch_indexes]
            y_2 = y_train[batch_indexes_2]
            x_2 = x_train[batch_indexes_2]
            if epoch < N_WARM_UP:
                for j in range(TRAINING_RATIO):
                    x_batch = x[j * BATCH_SIZE:(j + 1) * BATCH_SIZE]
                    y_batch = y[j * BATCH_SIZE:(j + 1) * BATCH_SIZE]

                    generator_loss.append(generator_model.train_on_batch([x_batch, positive_y], [y_batch, zero_y]))
            else:
                for j in range(TRAINING_RATIO):
                    x_batch = x[j * BATCH_SIZE:(j + 1) * BATCH_SIZE]
                    y_batch_2 = y_2[j * BATCH_SIZE:(j + 1) * BATCH_SIZE]
                    x_batch_2 = x_2[j * BATCH_SIZE:(j + 1) * BATCH_SIZE]
                    # noise = np.random.rand(BATCH_SIZE, 100).astype(np.float32)
                    noise = x_batch
                    #print(sum(y_batch_2))
                    discriminator_loss.append(discriminator_model.train_on_batch(
                        [y_batch_2,  noise, x_batch_2 ],
                        [positive_y, negative_y, dummy_y]))

                generator_loss.append(generator_model.train_on_batch([x,positive_full_y], [y, positive_full_y]))
        # Still needs some code to display losses from the generator and discriminator, progress bars, etc.
        predicted_y_train, _ = generator_model.predict([x_train , positive_full_enable_train], batch_size=BATCH_SIZE)
        predicted_y_val, _ = generator_model.predict([ x_val , positive_full_enable_val ], batch_size=BATCH_SIZE)
        
        #print(sum(sum(positive_full_enable_train)))
        #print(predicted_y_train)
        train_loss = log_loss(y_train, predicted_y_train)
        val_loss = log_loss(y_val, predicted_y_val)

        print("train loss: {:.4f}, validation loss: {:.4f}, discriminator loss: {:.4f}".format(
            train_loss, val_loss,
            (np.sum(np.asarray(discriminator_loss)) if discriminator_loss else -1) / x_train.shape[0]))

        if best_validation_loss is None or best_validation_loss > val_loss:
            print('\nEpoch %05d: improved from %0.5f,'
                  ' saving model to %s and %s'
                  % (epoch + 1, val_loss, generator_model_path, discriminator_model_path))

            best_validation_loss = val_loss
            generator_model.save(generator_model_path, overwrite=True)
            discriminator_model.save(discriminator_model_path, overwrite=True)


def model(params, batch_size=20, nb_epoch=40, is_train=True):
    # set parameters:
    nb_classes = len(functions)
    start_time = time.time()
    logging.info("Loading Data")
    

    train, val, test, train_df, valid_df, test_df = load_data()
    train_df = pd.concat([train_df, valid_df])
    test_gos = test_df['gos'].values
    train_data, train_labels = train
    val_data, val_labels = val
    test_data, test_labels = test


    logging.info("Data loaded in %d sec" % (time.time() - start_time))
    logging.info("Training data size: %d" % len(train_data))
    logging.info("Validation data size: %d" % len(val_data))
    logging.info("Test data size: %d" % len(test_data))
    generator_model_path = DATA_ROOT + 'models/new_model_seq_' + FUNCTION + '.h5'
    discriminator_model_path = DATA_ROOT + 'models/new_model_disc_seq_' + FUNCTION + '.h5'


    logging.info('Starting training the model')

    train_generator = DataGenerator(batch_size, nb_classes)
    train_generator.fit(train_data, train_labels)
    valid_generator = DataGenerator(batch_size, nb_classes)
    valid_generator.fit(val_data, val_labels)
    test_generator = DataGenerator(batch_size, nb_classes)
    test_generator.fit(test_data, test_labels)

    if is_train:
        generator_model, discriminator_model = get_model(params, nb_classes, batch_size)
        train_wgan(generator_model, discriminator_model, batch_size=batch_size, epochs=nb_epoch,
                   x_train=train_data, y_train=train_labels, x_val=val_data, y_val=val_labels,
                   generator_model_path=generator_model_path,
                   discriminator_model_path=discriminator_model_path)

    logging.info('Loading best model')
    model = load_model(generator_model_path,
                       custom_objects={'generator_recunstruction_loss_new': generator_recunstruction_loss_new,
                                       'wasserstein_loss': wasserstein_loss})

    logging.info('Predicting')
    preds = model.predict_generator(test_generator, steps=len(test_data) / batch_size)[0]

    logging.info('Computing performance')
    f, p, r, t, preds_max = compute_performance(preds, test_labels) #, test_gos)
    roc_auc = compute_roc(preds, test_labels)
    mcc = compute_mcc(preds_max, test_labels)
    aupr , _ = compute_aupr(preds, test_labels)
    m_pr_max, m_rc_max, m_f1_max, M_pr_max, M_rc_max, M_f1_max = micro_macro_function_centric_f1(preds.T, test_labels.T)

    logging.info('Protein centric macro Th, PR, RC, F1: \t %f %f %f %f' % (t, p, r, f))
    logging.info('ROC AUC: \t %f ' % (roc_auc, ))
    logging.info('MCC: \t %f ' % (mcc, ))
    logging.info('AUPR: \t %f ' % (aupr, ))
    logging.info('Function centric macro PR, RC, F1: \t %f %f %f'  % (M_pr_max, M_rc_max, M_f1_max) )
    logging.info('Function centric micro PR, RC, F1: \t %f %f %f'  % (m_pr_max, m_rc_max, m_f1_max) )

    function_centric_performance(functions, preds.T, test_labels.T, train_labels.T)



def load_prot_ipro():
    proteins = list()
    ipros = list()
    with open(DATA_ROOT + 'swissprot_ipro.tab') as f:
        for line in f:
            it = line.strip().split('\t')
            if len(it) != 3:
                continue
            prot = it[1]
            iprs = it[2].split(';')
            proteins.append(prot)
            ipros.append(iprs)
    return pd.DataFrame({'proteins': proteins, 'ipros': ipros})


def performanc_by_interpro():
    pred_df = pd.read_pickle(DATA_ROOT + 'test-' + FUNCTION + '-preds.pkl')
    ipro_df = load_prot_ipro()
    df = pred_df.merge(ipro_df, on='proteins', how='left')
    ipro = get_ipro()

    def reshape(values):
        values = np.hstack(values).reshape(
            len(values), len(values[0]))
        return values

    for ipro_id in ipro:
        if len(ipro[ipro_id]['parents']) > 0:
            continue
        labels = list()
        predictions = list()
        gos = list()
        for i, row in df.iterrows():
            if not isinstance(row['ipros'], list):
                continue
            if ipro_id in row['ipros']:
                labels.append(row['labels'])
                predictions.append(row['predictions'])
                gos.append(row['gos'])
        pr = 0
        rc = 0
        total = 0
        p_total = 0
        for i in range(len(labels)):
            tp = np.sum(labels[i] * predictions[i])
            fp = np.sum(predictions[i]) - tp
            fn = np.sum(labels[i]) - tp
            all_gos = set()
            for go_id in gos[i]:
                if go_id in all_functions:
                    all_gos |= get_anchestors(go, go_id)
            all_gos.discard(GO_ID)
            all_gos -= func_set
            fn += len(all_gos)
            if tp == 0 and fp == 0 and fn == 0:
                continue
            total += 1
            if tp != 0:
                p_total += 1
                precision = tp / (1.0 * (tp + fp))
                recall = tp / (1.0 * (tp + fn))
                pr += precision
                rc += recall
        if total > 0 and p_total > 0:
            rc /= total
            pr /= p_total
            if pr + rc > 0:
                f = 2 * pr * rc / (pr + rc)
                logging.info('%s\t%d\t%f\t%f\t%f' % (
                    ipro_id, len(labels), f, pr, rc))


def function_centric_performance(functions, preds, labels, labels_train):
    results = []
    preds = np.round(preds, 2)
    for i in range(preds.shape[0]):
        f_max = 0
        p_max = 0
        r_max = 0
        for t in range(1, 100):
            threshold = t / 100.0
            predictions = (preds[i, :] > threshold).astype(np.int32)
            tp = np.sum(predictions * labels[i, :])
            fp = np.sum(predictions) - tp
            fn = np.sum(labels[i, :]) - tp
            if tp > 0:
                precision = tp / (1.0 * (tp + fp))
                recall = tp / (1.0 * (tp + fn))
                f = 2 * precision * recall / (precision + recall)
            else:
                if fp == 0 and fn == 0:
                    precision = 1
                    recall = 1
                    f = 1
                else:
                    precision = 0
                    recall = 0
                    f = 0

            if f_max < f:
                f_max = f
                p_max = precision
                r_max = recall
        num_prots_train = np.sum(labels_train[i, :])
        height = get_height(go, functions[i])
        results.append([functions[i], num_prots_train, height, f_max, p_max, r_max])
    results = pd.DataFrame(results)
    results.to_csv('Con_GodGanSeq_results_' + FUNCTION + '.txt', sep='\t', index=False)



def function_centric_performance_backup(functions, preds, labels, labels_train):
    results = []
    preds = np.round(preds, 2)
    for i in range(len(functions)):
        f_max = 0
        p_max = 0
        r_max = 0
        x = list()
        y = list()
        total = 0
        for t in range(1, 100):
            threshold = t / 100.0
            predictions = (preds[i, :] > threshold).astype(np.int32)
            tp = np.sum(predictions * labels[i, :])
            fp = np.sum(predictions) - tp
            fn = np.sum(labels[i, :]) - tp
            if tp >0:    
                sn = tp / (1.0 * np.sum(labels[i, :]))
                sp = np.sum((predictions ^ 1) * (labels[i, :] ^ 1))
                sp /= 1.0 * np.sum(labels[i, :] ^ 1)
                fpr = 1 - sp
                x.append(fpr)
                y.append(sn)
                precision = tp / (1.0 * (tp + fp))
                recall = tp / (1.0 * (tp + fn))
                f = 2 * precision * recall / (precision + recall)
                total +=1
                if f_max < f:
                    f_max = f
                    p_max = precision
                    r_max = recall
        num_prots = np.sum(labels[i, :])
        num_prots_train = np.sum(labels_train[i,:])
        if total >1 :
            roc_auc = auc(x, y)
        else:
            roc_auc =0
        height = get_height(go , functions[i])
        results.append([functions[i], f_max, p_max, r_max, num_prots, num_prots_train, height,roc_auc])
    results = pd.DataFrame(results)    
    #results.to_csv('new_results.txt' , sep='\t' , index = False)
    results.to_csv('Con_GodGanSeq_results_'+FUNCTION +'.txt', sep='\t', index=False)
    #results = np.array(results)
    #p_mean = (np.sum(results[:,2])) / len(functions)
    #r_mean =  (np.sum(results[:,3])) / len(functions)
    #f_mean = (2*p_mean*r_mean)/(p_mean+r_mean)
    #roc_auc_mean = (np.sum(results[:,7])) / len(functions)
    #print('Function centric performance (macro) ' '%f %f %f %f' % (f_mean, p_mean, r_mean, roc_auc_mean))



def micro_macro_function_centric_f1_backup(preds, labels):
    preds = np.round(preds, 2)
    m_f1_max = 0
    M_f1_max = 0
    for t in range(1, 100):
        threshold = t / 100.0
        predictions = (preds > threshold).astype(np.int32)
        m_tp = 0
        m_fp = 0
        m_fn = 0
        M_pr = 0
        M_rc = 0
        total = 0
        p_total = 0
        for i in range(len(preds)):
            tp = np.sum(predictions[i, :] * labels[i, :])
            fp = np.sum(predictions[i, :]) - tp
            fn = np.sum(labels[i, :]) - tp
            if tp == 0 and fp == 0 and fn == 0:
                continue
            total += 1
            if tp > 0:
                pr = tp / (1.0 * (tp + fp))
                rc = tp / (1.0 * (tp + fn))
                m_tp += tp
                m_fp += fp
                m_fn += fn
                M_pr += pr
                M_rc += rc
                p_total += 1

        if p_total == 0:
            continue
        if total > 0:
            m_tp /= total
            m_fn /= total
            m_fp /= total
            m_pr = m_tp / (1.0 * (m_tp + m_fp))
            m_rc = m_tp / (1.0 * (m_tp + m_fn))
            M_pr /= p_total
            M_rc /= total
            m_f1 = 2 * m_pr * m_rc / (m_pr + m_rc)
            M_f1 = 2 * M_pr * M_rc / (M_pr + M_rc)

            if m_f1 > m_f1_max:
                m_f1_max = m_f1
                m_pr_max = m_pr
                m_rc_max = m_rc

            if M_f1 > M_f1_max:
                M_f1_max = M_f1
                M_pr_max = M_pr
                M_rc_max = M_rc

    return m_pr_max, m_rc_max, m_f1_max, M_pr_max, M_rc_max, M_f1_max





def micro_macro_function_centric_f1(preds, labels):
    preds = np.round(preds, 2)
    m_f1_max = 0
    M_f1_max = 0
    for t in range(1, 200):
        threshold = t / 200.0
        predictions = (preds > threshold).astype(np.int32)
        m_tp = 0
        m_fp = 0
        m_fn = 0
        M_pr = 0
        M_rc = 0
        for i in range(preds.shape[0]):
            tp = np.sum(predictions[i, :] * labels[i, :])
            fp = np.sum(predictions[i, :]) - tp
            fn = np.sum(labels[i, :]) - tp
            m_tp += tp
            m_fp += fp
            m_fn += fn
            if tp > 0:
                pr = 1.0 * tp / (1.0 * (tp + fp))
                rc = 1.0 * tp / (1.0 * (tp + fn))
            else:
                if fp == 0 and fn == 0:
                    pr = 1
                    rc = 1
                else:
                    pr = 0
                    rc = 0
            M_pr += pr
            M_rc += rc

        if m_tp > 0:
            m_pr = 1.0 * m_tp / (1.0 * (m_tp + m_fp))
            m_rc = 1.0 * m_tp / (1.0 * (m_tp + m_fn))
            m_f1 = 2.0 * m_pr * m_rc / (m_pr + m_rc)
        else:
            if m_fp == 0 and m_fn == 0:
                m_pr = 1
                m_rc = 1
                m_f1 = 1
            else:
                m_pr = 0
                m_rc = 0
                m_f1 = 0

        M_pr /= preds.shape[0]
        M_rc /= preds.shape[0]
        if M_pr == 0 and M_rc == 0:
            M_f1 = 0
        else:
            M_f1 = 2.0 * M_pr * M_rc / (M_pr + M_rc)
        if m_f1 > m_f1_max:
            m_f1_max = m_f1
            m_pr_max = m_pr
            m_rc_max = m_rc

        if M_f1 > M_f1_max:
            M_f1_max = M_f1
            M_pr_max = M_pr
            M_rc_max = M_rc

    return m_pr_max, m_rc_max, m_f1_max, M_pr_max, M_rc_max, M_f1_max


def compute_roc(preds, labels):
    # Compute ROC curve and ROC area for each class
    fpr, tpr, _ = roc_curve(labels.flatten(), preds.flatten())
    roc_auc = auc(fpr, tpr)
    return roc_auc


def compute_aupr(preds, labels):
    # Compute ROC curve and ROC area for each class
    pr, rc, threshold =precision_recall_curve(labels.flatten(), preds.flatten())
    pr_auc = auc(rc, pr)
    #pr, rc, threshold =precision_recall_curve(labels.flatten(), preds.flatten(),average ='macro' )
    M_pr_auc = 0 
    return pr_auc, M_pr_auc



def compute_mcc(preds, labels):
    # Compute ROC curve and ROC area for each class
    mcc = matthews_corrcoef(labels.flatten(), preds.flatten())
    return mcc


def compute_performance(preds, labels): #, gos):
    preds = np.round(preds, 2)
    f_max = 0
    p_max = 0
    r_max = 0
    t_max = 0
    for t in range(1, 100):
        threshold = t / 100.0
        predictions = (preds > threshold).astype(np.int32)
        total = 0
        f = 0.0
        p = 0.0
        r = 0.0
        p_total = 0
        for i in range(labels.shape[0]):
            tp = np.sum(predictions[i, :] * labels[i, :])
            fp = np.sum(predictions[i, :]) - tp
            fn = np.sum(labels[i, :]) - tp
            all_gos = set()
            #for go_id in gos[i]:
            #    if go_id in all_functions:
            #        all_gos |= get_anchestors(go, go_id)
            #all_gos.discard(GO_ID)
            #all_gos -= func_set
            #fn += len(all_gos)
            if tp == 0 and fp == 0 and fn == 0:
                continue
            total += 1
            if tp != 0:
                p_total += 1
                precision = tp / (1.0 * (tp + fp))
                recall = tp / (1.0 * (tp + fn))
                p += precision
                r += recall
        if p_total == 0:
            continue
        r /= total
        p /= p_total
        if p + r > 0:
            f = 2 * p * r / (p + r)
            if f_max < f:
                f_max = f
                p_max = p
                r_max = r
                t_max = threshold
                predictions_max = predictions
    return f_max, p_max, r_max, t_max, predictions_max


def get_gos(pred):
    mdist = 1.0
    mgos = None
    for i in range(len(labels_gos)):
        labels, gos = labels_gos[i]
        dist = distance.cosine(pred, labels)
        if mdist > dist:
            mdist = dist
            mgos = gos
    return mgos


def compute_similarity_performance(train_df, test_df, preds):
    logging.info("Computing similarity performance")
    logging.info("Training data size %d" % len(train_df))
    train_labels = train_df['labels'].values
    train_gos = train_df['gos'].values
    global labels_gos
    labels_gos = zip(train_labels, train_gos)
    p = Pool(64)
    pred_gos = p.map(get_gos, preds)
    total = 0
    p = 0.0
    r = 0.0
    f = 0.0
    test_gos = test_df['gos'].values
    for gos, tgos in zip(pred_gos, test_gos):
        preds = set()
        test = set()
        for go_id in gos:
            if go_id in all_functions:
                preds |= get_anchestors(go, go_id)
        for go_id in tgos:
            if go_id in all_functions:
                test |= get_anchestors(go, go_id)
        tp = len(preds.intersection(test))
        fp = len(preds - test)
        fn = len(test - preds)
        if tp == 0 and fp == 0 and fn == 0:
            continue
        total += 1
        if tp != 0:
            precision = tp / (1.0 * (tp + fp))
            recall = tp / (1.0 * (tp + fn))
            p += precision
            r += recall
            f += 2 * precision * recall / (precision + recall)
    return f / total, p / total, r / total


def print_report(report, go_id):
    with open(DATA_ROOT + 'reports.txt', 'a') as f:
        f.write('Classification report for ' + go_id + '\n')
        f.write(report + '\n')


if __name__ == '__main__':
    main()