# -*- coding: utf-8 -*-
"""
Created on Mon Sep 14 17:41:37 2020

@author: Lucas & Caro
"""



import streamlit as st
import base64

from pathlib import Path
from sklearn.metrics import roc_curve, precision_recall_curve, average_precision_score
from sklearn.model_selection import train_test_split
from rdkit.ML.Scoring import Scoring
from sklearn.metrics import confusion_matrix
import math
import pandas as pd
import statistics
# import plotly as py
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import numpy as np

from sklearn.metrics import accuracy_score
from sklearn.metrics import balanced_accuracy_score
from sklearn.metrics import f1_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import matthews_corrcoef
#---------------------------------#
# Page layout
## Page expands to full width
st.set_page_config(page_title='LIDEB Tools - Metrics', page_icon="📏", layout='wide')

######
# Function to put a picture as header   
def img_to_bytes(img_path):
    img_bytes = Path(img_path).read_bytes()
    encoded = base64.b64encode(img_bytes).decode()
    return encoded

from PIL import Image
image = Image.open('cropped-header.png')
st.image(image)
#####
st.write("&nbsp[![Website](https://img.shields.io/badge/website-LIDeB-blue)](https://lideb.biol.unlp.edu.ar)&nbsp[![Twitter Follow](https://img.shields.io/twitter/follow/LIDeB_UNLP?style=social)](https://twitter.com/intent/follow?screen_name=LIDeB_UNLP)")
st.markdown("###### Made in 🐍 and [![this is an image link](https://i.imgur.com/iIOA6kU.png)](https://www.streamlit.io/)&nbsp with :heart: by [Lucas Alberca](https://twitter.com/capigol) and [Caro Bellera](https://twitter.com/carobellera)")
#####
st.subheader(":pushpin:" "About Us")
st.markdown("We are a drug discovery team with an interest in the development of publicly available open-source customizable cheminformatics tools to be used in computer-assisted drug discovery. We belong to the Laboratory of Bioactive Research and Development (LIDeB) of the National University of La Plata (UNLP), Argentina. Our research group is focused on computer-guided drug repurposing and rational discovery of new drug candidates to treat epilepsy and neglected tropical diseases.")


# Introduction
#---------------------------------#

st.write("""
# LIDeB Tools - Metrics
**Metrics to evaluate the performance of the classificatory models **

Web App to evaluate the perfomance of classificatory models by calculation of AUCROC, BEDROC, Accuracy, F-measure, PR and EF.

Additionaly, using the generated **PPV surfaces** you can select an adequate score threshold for prospective virtual screenings.

The tool uses the following packages [RDKIT](https://www.rdkit.org/docs/index.html), [Scikit-learn](https://scikit-learn.org/stable/), [Plotly](https://plotly.com/python/)
    
""")


# image = Image.open('clustering_workflow.png')
# st.image(image, caption='Clustering Workflow')

st.sidebar.header('Upload your File')

uploaded_file_1 = st.sidebar.file_uploader("Upload a TXT file with the scores", type=["txt"])
st.sidebar.markdown("""
[Example TXT input file](https://raw.githubusercontent.com/cbellera/metrics/main/example_file_metrics.txt)
""")

st.sidebar.header('SCORE THRESHOLD')    
threshold_cm = st.sidebar.number_input('Threshold to split between actives and inactives', value=0.5)

st.sidebar.header('BEDROC')    
alpha_bedroc = st.sidebar.slider('Alpha BEDROC', 5, 100, 20, 5)

st.sidebar.header('ENRICHMENT FACTOR')    
fraccion_enriquesimiento = st.sidebar.slider('Fraction enrichment', 0.0, 0.95, 0.01 , 0.01)

st.sidebar.header('STANDARD DEVIATION')    
repeticiones = st.sidebar.slider('Number of iterations', 100, 1000, 100 , 50)
particion = st.sidebar.slider('Fraction of evaluation', 0.5, 0.95, 0.85 , 0.05)




st.sidebar.title(":speech_balloon: Contact Us")
st.sidebar.info(
"""
If you are looking to contact us, please
[:e-mail:](mailto:lideb@biol.unlp.edu.ar) or [Twitter](https://twitter.com/LIDeB_UNLP)
""")



def metrics_calculation(uploaded_file_1):
    
    # Leemos el archivo cargado
    scores1 = uploaded_file_1
    # Lo ordenamos
    scores1.sort_values(by=["score"],axis = 0, ascending = False,inplace = True, na_position ='last')
    
    y_modelo = list(scores1["score"])
    y_real = list(scores1["class"])
    
    clases_predichas = []
    for value in y_modelo:
        if value > threshold_cm:
            clases_predichas.append(1)
        else:
            clases_predichas.append(0)
            
    clases_predichas_serie = pd.Series(clases_predichas)
    # result = scores1["class"].eq(clases_predichas_serie)
    # proporcion_bc = round(sum(result) / len(result),4) # Numero de compuestos bien clasificados / numero de compuestos del dataset
    
    st.markdown(":point_down: **Here you can see the perfomance of the model for the selected threshold**", unsafe_allow_html=True)
    
    metric1 = round(accuracy_score(y_real, clases_predichas_serie),4)
    metric2 = round(balanced_accuracy_score(y_real, clases_predichas_serie),4)
    metric4 = round(f1_score(y_real, clases_predichas_serie),4)
    metric5 = round(precision_score(y_real, clases_predichas_serie),4)
    metric6 = round(recall_score(y_real, clases_predichas_serie),4)
    metric7 = round(matthews_corrcoef(y_real, clases_predichas_serie),4)
    
    results = [metric1,metric2,metric4,metric5,metric6,metric7]
    metric_name = ["Accuracy","Balanced accuracy","F-measure","Precision","Recall", "MCC"]
    table_metrics = pd.Series(results,index=metric_name,name="Value")
    st.write(table_metrics)
   
    st.markdown(":point_down: **Here you can find the confusion matrix for the selected threshold**", unsafe_allow_html=True)
    matrix = pd.DataFrame(confusion_matrix(y_real, clases_predichas))
    # matrix=matrix.rename(index={0: 'Predicted Positive', 1:"Predicted Negative"},columns={0:"Real Positive",1:"Real Negative"})
    matrix=matrix.rename(index={0: 'Real Negative', 1:"Real Positive"},columns={0:"Predicted Negative",1:"Predicted Positive"})

    st.write(matrix)
    
    # Hacemos los calculos para X repeticiones
    i=0   
    aucs=[]
    bedrocs=[]
    efs=[]
    prs=[]

    # aupr_max= 1
    
    while i <= repeticiones:    # Para hacer el cross validation y poder sacar la desviacion del auc
    
        scores_modelo, X_test, clases_real, y_test = train_test_split(y_modelo, y_real, test_size= 1 - particion, random_state=i,stratify = y_real)
            
        array = tuple(zip(scores_modelo,clases_real))
        array_1 = pd.DataFrame(array)
        array_1.sort_values(0, axis = 0, ascending = False,inplace = True, na_position ='last')
        array_metricas = tuple(zip(array_1[0],array_1[1]))
        
        auc = Scoring.CalcAUC(scores= array_metricas, col= 1)   
        aucs.append(auc)
        bedroc = Scoring.CalcBEDROC(scores= array_metricas, col= 1,alpha=alpha_bedroc)
        bedrocs.append(bedroc)
        enrichment = Scoring.CalcEnrichment(scores= array_metricas, col= 1,fractions=[fraccion_enriquesimiento])
        efs.append(enrichment[0])
        auc_PR = average_precision_score(clases_real, scores_modelo)                        # Calculo el AUPR
        prs.append(auc_PR)
        
        i = i + 1
        
    
    # # Esto ese para el EFmax #
    n_total = len(y_real)           # TOT(tot)
    fraccion_ef = math.ceil(n_total*fraccion_enriquesimiento) # devuelve el numero entero mayor (TOT(1%))
    n_activos = sum(y_real)     # ACT(TOT)
    if fraccion_ef > n_activos:
        efmax = (n_activos/fraccion_ef)/(n_activos/n_total)
    else:
        efmax = 1/(n_activos/n_total)
    
    # Sacamos las medias y las desviaciones
    auc_ok= round(statistics.mean(aucs),4)
    sd_auc = round(statistics.stdev(aucs),4)
    
    bedroc_ok= round(statistics.mean(bedrocs),4)
    sd_bedroc = round(statistics.stdev(bedrocs), 4)
    
    efs1 = [x/efmax for x in efs]
    ef_ok= round(statistics.mean(efs1),4)
    if ef_ok > 1:
        ef_ok = 1.000
    sd_ef = round(statistics.stdev(efs1),4)
    
    pr_ok= round(statistics.mean(prs),4)
    sd_pr = round(statistics.stdev(prs),4)
    
    resultado = [auc_ok,bedroc_ok,ef_ok, pr_ok]
    desviaciones = [sd_auc,sd_bedroc, sd_ef, sd_pr]
    
    final = pd.DataFrame([resultado,desviaciones]).T
    
    final=final.rename(index={0: 'AUC', 1:"BEDROC_" + str(alpha_bedroc), 2: 'EFrel_' + str(fraccion_enriquesimiento), 3:'AUC_PR'},columns={0:"Mean",1:"SD"})
    
    st.markdown(":point_down: **Here you can see the calculated metrics**", unsafe_allow_html=True)
    st.write(final)
    # st.markdown(":point_down: **Here you can download this metrics**", unsafe_allow_html=True)
    st.markdown(filedownload_1(final), unsafe_allow_html=True)

    
    return scores1, auc_ok, sd_auc,pr_ok, sd_pr, y_modelo, y_real

def roc_aupr_curves(auc_ok, sd_auc, pr_ok, sd_pr , y_modelo, y_real):

    title1= "ROC CURVE " + " " + "(" + "AUCROC = " + str(format(auc_ok, '.4f')) + " ± " + str(format(sd_auc, '.4f')) + ")"
    title2="PR CURVE " + " " + "(" + "AUPR = " + str(format(pr_ok, '.4f')) + " ± " + str(format(sd_pr, '.4f')) + ")"
    fig = make_subplots(rows=1, cols=2, subplot_titles=(title1,title2))
    
    
    # # CURVA ROC
    fpr, tpr, thresholds = roc_curve(y_score= y_modelo,y_true= y_real)
    fig.add_trace(go.Scatter(x=fpr, y=tpr, fill='tozeroy', showlegend=False),row=1, col=1)
    fig.add_shape(type='line', line=dict(dash='dash'), x0=0, x1=1, y0=0, y1=1,row=1, col=1)
    fig.update_xaxes(title_text="False Positive Rate (FPR)", row=1, col=1)
    fig.update_yaxes(title_text="True Positive Rate (TPR)", row=1, col=1)

    
    # # CURVA PRECISION RECALL
    precision, recall, thresholds = precision_recall_curve(y_true= y_real, probas_pred = y_modelo)
    fig.add_trace(go.Scatter(x=recall, y=precision, fill='tozeroy', showlegend=False),row=1, col=2)
    fig.add_shape(type='line', line=dict(dash='dash'),x0=0, x1=1, y0=1, y1=0,row=1, col=2)
    
    fig.update_layout(margin = dict(t=60,r=20,b=20,l=20), autosize = True)
    fig.update_xaxes(title_text="Recall", row=1, col=2)
    fig.update_yaxes(title_text="Precision", row=1, col=2)
    st.plotly_chart(fig)


def ppv_plot(scores1):  
# importamos los paquetes necesarios:
    y_modelo = list(scores1["score"])
    y_real = list(scores1["class"])
    
    # Determino la sensibilidad, especificidad y los valores de corte
    fpr, tpr, thresholds = roc_curve(y_score= y_modelo, y_true= y_real,drop_intermediate=False)  
    se = tpr
    sp = 1 - fpr
    se_sp = se/sp 
     
    tabla_final = pd.DataFrame(list(zip(thresholds,fpr,tpr,se,sp,se_sp)), 
                   columns =['Cutoff', 'fpr',"tpr","Se","Sp","Se/Sp"]) 
    
    # Calculos para agregar el valor de PPV a dos Ya distintos a la tabla
    
    prev1 = [0.001,0.005,0.01]
    
    se1=tabla_final["Se"].tolist()
    sp1=tabla_final["Sp"].tolist()
    
    tamanio = list(range(1, len(se1), 1))
    
    ppv_ok1 =[]
    
    for prev in prev1:
        ppv1a=[]
        for i in tamanio:
            ppv1 = (se1[i]*prev)/(se1[i]*prev + (1 - sp1[i])*(1-prev))
            ppv1a.append(ppv1)
        ppv_ok1.append(ppv1a)    
      
    # ppv_para_dos_ya = pd.DataFrame(ppv_ok1).T
    
    tabla_final_1 = pd.DataFrame(list(zip(thresholds,se,sp,se_sp,ppv_ok1[0],ppv_ok1[1],ppv_ok1[2])), 
                    columns =['Cutoff', "Se","Sp","Se/Sp","PPV (Ya = 0.001)","PPV (Ya = 0.005)","PPV (Ya = 0.01)"]) 
    
    
    # preparacion de datos para el grafico de  ppv
    
    prevalencia = pd.DataFrame(np.arange(0.00000000001, 0.012, 0.001))
    ya=prevalencia[0].tolist()
    
    se1=tabla_final["Se"].tolist()
    sp1=tabla_final["Sp"].tolist()
    
    tamanio = list(range(1, len(se1), 1))
    
    ppv_ok =[]
    
    for x in ya:
        ppv1=[]
        for i in tamanio:
            ppv = (se1[i]*x)/(se1[i]*x + (1 - sp1[i])*(1-x))
            ppv1.append(ppv)
        ppv_ok.append(ppv1)    
      
    ppv_tabla = pd.DataFrame(ppv_ok).T
    
    
    # GRAFICO PPV
    
    fig = go.Figure(go.Surface(opacity=1, 
        contours = {
            "x": {"show": False, "start": 0, "end": 0.01, "size": 0.04, "color":"white"},
            "y": {"show": False, "start": 0, "end": 2, "size": 0.04, "color":"white"},
            "z": {"show": False, "start": 0, "end": 1, "size": 0.2,"color":"white"}
        },
        x = ya,
        y = se_sp,
        z = ppv_tabla))
    
    fig.update_layout(
          scene = dict(
            xaxis = dict(title="Ya", nticks=10, range=[0,0.011],gridwidth = 10,linewidth = 10,
                         backgroundcolor="rgb(200, 200, 230)", gridcolor="white", showbackground=True, zerolinecolor="white"),
            yaxis = dict(title="Se/Sp",nticks=10, range=[0,2],gridwidth = 10,linewidth = 10,
                         backgroundcolor="rgb(230, 200,230)",
                            gridcolor="white",
                            showbackground=True,
                            zerolinecolor="white"),
            zaxis = dict(title="PPV",nticks=10, range=[0,1],gridwidth = 10,linewidth = 10,
                         backgroundcolor="rgb(230, 230,200)",
                            gridcolor="white",
                            showbackground=True,
                            zerolinecolor="white"),),
        margin=dict(r=20, l=10, b=10, t=10),
        font=dict(
            size=14,
            color="black")
        )
    st.plotly_chart(fig)
    
    st.markdown(":point_down: **Here you can see the table with cutoffs, Se, Sp and PPVs**", unsafe_allow_html=True)
    st.write(tabla_final_1)
    return tabla_final_1

#%%
### Exporting files ###

def filedownload(df):
    csv = df.to_csv(index=False,header=True)
    b64 = base64.b64encode(csv.encode()).decode()  # strings <-> bytes conversions
    href = f'<a href="data:file/csv;base64,{b64}" download="score_thresholds.csv">Download CSV File with the previous table</a>'
    return href

def filedownload_1(df):
    csv = df.to_csv(index=True,header=True)
    b64 = base64.b64encode(csv.encode()).decode()  # strings <-> bytes conversions
    href = f'<a href="data:file/csv;base64,{b64}" download="model_metrics.csv">Download CSV File with the metrics</a>'
    return href





if uploaded_file_1 is not None:
    scores = pd.read_csv(uploaded_file_1, sep='\t', delimiter=None, header='infer', names=None)

    scores1, auc_ok, sd_auc, pr_ok, sd_pr , y_modelo, y_real = metrics_calculation(scores)
    st.write("")
    st.markdown(":point_down: **Here you can see the ROC and AUPR curves**", unsafe_allow_html=True)
    roc_aupr_curves(auc_ok, sd_auc, pr_ok, sd_pr , y_modelo, y_real)
        
    st.markdown(":point_down: **Here you can see the PPV surface**", unsafe_allow_html=True)
    tabla_final_1 = ppv_plot(scores1)
    
    st.markdown(filedownload(tabla_final_1), unsafe_allow_html=True)
    
else:
    st.info('Awaiting for TXT file to be uploaded.')
    if st.button('Press to use Example File'):
        st.write("A file with 1132 predictions (only 22 real actives) have been loaded as example")
        scores = pd.read_csv("example_file_metrics.txt",sep="\t",header="infer")
        scores1, auc_ok, sd_auc, pr_ok, sd_pr , y_modelo, y_real = metrics_calculation(scores)
        st.write("")
        st.markdown(":point_down: **Here you can see the ROC and AUPR curves**", unsafe_allow_html=True)
        roc_aupr_curves(auc_ok, sd_auc, pr_ok, sd_pr , y_modelo, y_real)
            
        st.markdown(":point_down: **Here you can see the PPV surface**", unsafe_allow_html=True)
        tabla_final_1 = ppv_plot(scores1)
        
        st.markdown(filedownload(tabla_final_1), unsafe_allow_html=True)

   
    

    
