# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# # Experiment 02: Study performance stability over time
#
# Study how well models trained on images from early dates perform at test time on images from later dates.  This is meant to investigate how stable model performance is over time, as news rooms' image publishing pipelines (possibly) evolve.  
#
# For each source, sort the images chronologically by the news article date, then split the images into a training subset (with images from early dates), and a test set (with images from later dates.)  
#
# Then train models using the images from early dates, and then test the models on images from late dates.
#
# Only include QM features since they greatly outperformed CL features.  Only Study Naive Bayes model here, so the focus is on the effect of time, not the effect of the model (and since NB was a top performing model.)

# %matplotlib widget
# %load_ext autoreload
# %autoreload 2

# +
import os
import sys
import subprocess
import random
import pickle
import numpy as np
import pandas as pd

from sklearn import preprocessing
from sklearn.model_selection import train_test_split

import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
# from tqdm.autonotebook import tqdm
from tqdm.notebook import tqdm

import uncertainties

from image_compression_attribution.common.code.models import quant_matrices, compr_levels
from image_compression_attribution.common.code.summarize_quant_matrices import summarize_compression_features

pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 500)
pd.set_option('display.max_colwidth', 500)

from sklearn.metrics import make_scorer, roc_curve
from scipy.optimize import brentq
from scipy.interpolate import interp1d

#WARNING: this method does not seem to work well when there are large gaps
#in the ROC curve. Hence, only use this if you have interpolated between
#ROC curve data points to fill in the roc curve on a grid with small intervals.
#https://github.com/scikit-learn/scikit-learn/issues/15247#issuecomment-542138349
def calculate_eer(fpr, tpr):
    '''
    Returns the equal error rate for a binary classifier output.
    '''
    eer = brentq(lambda x : 1. - x - interp1d(fpr, tpr)(x), 0., 1.)
    return eer
  
#---------------------------------------------------------------
#Code to combine mean value and uncertainty estimate into
#one formatted string, like 3.14 +/- .02 becomes "3.14(2)"

import string

class ShorthandFormatter(string.Formatter):
  """https://pythonhosted.org/uncertainties/user_guide.html"""

  def format_field(self, value, format_spec):
    if isinstance(value, uncertainties.UFloat):
      return value.format(format_spec+'S')  # Shorthand option added
    # Special formatting for other types can be added here (floats, etc.)
    else:
      # Usual formatting:
      return super(ShorthandFormatter, self).format_field(
        value, format_spec)

def uncertainty_format_arrays(mean_vals, uncertainty_vals):
  frmtr_uncertainty = ShorthandFormatter()
  
  vals_formatted = []
  for mean, uncert in zip(mean_vals, uncertainty_vals):
    number = uncertainties.ufloat(mean, uncert)
    str_formatted = frmtr_uncertainty.format("{0:.1u}", number)
    vals_formatted.append(str_formatted)
  return vals_formatted


# -

RND_SEED=1234
np.random.seed(RND_SEED)

# +
SUMMARY_FILE = "/app/dataset/data.csv"

RESULTS_FOLDER = "results/exp_02"
os.makedirs(RESULTS_FOLDER, exist_ok=True)

df = pd.read_csv(SUMMARY_FILE)
#We'll work with timestamps, so need to convert to a datetime for ease of use
df['timestamp'] = pd.to_datetime(df['timestamp'], utc=True)

#Drop non-image files, e.g. html files returned
#due to download errors

df, df_dropped = df[ df['mime'].str.startswith('image') ].reset_index(drop=True), \
  df[ ~df['mime'].str.startswith('image') ].reset_index(drop=True)

sources = sorted(list(df['source'].unique()))


# +
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_curve, auc
from sklearn.ensemble import IsolationForest
#Guide to LabelEncoder:  
#https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.LabelEncoder.html


#create numerical class labels for quantization categorical names (suitable for 
#use as  ML training feature vector)
le_qs = preprocessing.LabelEncoder()
le_qs.fit(df['q_name'])
df['q_name_class'] = le_qs.transform(df['q_name'])

sources = sorted(list(df['source'].unique()))
le_sources = preprocessing.LabelEncoder()
le_sources.fit(sources)
df['source_class'] = le_sources.transform(df['source'])
# -

df





df.groupby('source')['timestamp'].min()

df.groupby('source')['timestamp'].max()

# +
timespan = df.groupby('source')['timestamp'].max() - df.groupby('source')['timestamp'].min()
timespan_list = timespan.tolist()

timespan_list2 = [(x.days + x.seconds/86400.0)/365.2425 for x in timespan_list]
timespan_years = pd.Series(timespan_list2, index=sources)

print("mean timespan = {:.1f} years".format(timespan_years.mean()))
with open(os.path.join(RESULTS_FOLDER,"timespan_years_mean.txt"),"w") as file1:
  file1.write("{:.1f}".format(timespan_years.mean()))

print("min timespan = {:.2f} years".format(timespan_years.min()))
with open(os.path.join(RESULTS_FOLDER,"timespan_years_min.txt"),"w") as file1:
  file1.write("{:.2f}".format(timespan_years.min()))

print("max timespan = {:.1f} years".format(timespan_years.max()))
with open(os.path.join(RESULTS_FOLDER,"timespan_years_max.txt"),"w") as file1:
  file1.write("{:.1f}".format(timespan_years.max()))

timespan_years
# -

# ## Comment:
# We see that for most sources we have images across large time spans, which is desireable for these experiments; the average time span is over 10 years, and the timespans range from 0.89 years to 21.4 years









# ## Method:
# Note: date ranges of different sources may not overlap. So instead of picking 1 cutoff date to be shared for all sources, we will instead have a separate cutoff date for each source, to split each source into two parts. Put another way, sort articles from each source by date, and within a source, split the articles into two parts: before and after the date cutoff.
# 1. All articles (from each source) from before the cutoff date form the train set --  the first 60% of articles.
# 1. Then form the test set:
#     1. Select all articles from a source after the cutoff date -- last 40% of articles.
#     1. Randomly sample an equal number of articles from the remaining sources, also after each's cutoff date.
#     1. Combine them to form the test set.
# 1. Since the composition of the test set varies, repeat this 5x to quantify uncertainty.



# First do some precomputation: For each source, sort the articles by date, then split articles from each source into an early portion and later portion. Early portions can be used for training, later portions can be used for testing.

# +
PERCENT_TEST = 0.40

df_articles = df[['articleHash', 'timestamp', 'source', 'source_class']].drop_duplicates()

articles_predate = {}
articles_postdate = {}

for source in sources:

  #get all articles from the source, sorted by date
  df_articles_from_source = df_articles[df_articles['source']==source].sort_values(by="timestamp")

  num_test_articles_from_source = int(PERCENT_TEST*len(df_articles_from_source))
  num_train_articles_from_source = len(df_articles_from_source) - num_test_articles_from_source

  df_art_from_source_predate = df_articles_from_source.iloc[0:num_train_articles_from_source,:]
  df_art_from_source_postdate = df_articles_from_source.iloc[num_train_articles_from_source:,:]

  articles_predate[source] = df_art_from_source_predate
  articles_postdate[source] = df_art_from_source_postdate
# -



# +
#Prepare Train and Test Split.

all_q_name_vals = sorted(df['q_name'].unique())

#Sample from articles (so we can keep images from articles grouped together)
df_articles = df[['articleHash', 'timestamp', 'source', 'source_class']].drop_duplicates()

NUM_TRIALS = 5
results_per_trial_qm = {}

for trial in tqdm(range(NUM_TRIALS)):
  
  numsamples_balanced_testset=[]
  AUCs_qm = []
  results_qm={}

  for source in sources:

    remaining_sources = [s for s in sources if s != source]

    #-----------------------------------
    #Form train/test data split. Train set first:

    #All articles (from every source) from before their cutoff date form the train set
    df_train_articles = None
    for src in sources:
      if df_train_articles is None:
        df_train_articles = articles_predate[src]
      else:
        df_train_articles = pd.concat([df_train_articles, articles_predate[src] ])

    #----------------------
    #Test set:

    #All articles from a source after its cutoff date contributes to test set:
    df_articles_test_from_source = articles_postdate[source]
    num_test_articles_from_source = len(df_articles_test_from_source)

    #-------

    #collect all articles not from remaining sources from after their cutoff dates
    df_articles_postdate_not_from_source = None
    for src in remaining_sources:
      if df_articles_postdate_not_from_source is None:
        df_articles_postdate_not_from_source = articles_postdate[src]
      else:
        df_articles_postdate_not_from_source = pd.concat([df_articles_postdate_not_from_source, articles_postdate[src] ])

    #Randomly sample an equal number of articles from the remaining sources, after their cutoff dates.
    num_test_articles_not_from_source = num_test_articles_from_source
    df_articles_test_not_from_source = df_articles_postdate_not_from_source.sample(num_test_articles_not_from_source)

    #------
    #combine to build the test set
    df_test_articles = pd.concat([df_articles_test_from_source, df_articles_test_not_from_source])

    #----------------------
    #Get all images articles in train/test splits:
    df_train = df[ df['articleHash'].isin(df_train_articles['articleHash']) ].reset_index()
    df_test = df[ df['articleHash'].isin(df_test_articles['articleHash']) ].reset_index()

    #Set ground truth label: 1 if image misattributed, else 0
    df_test['is_misattributed'] = np.array(df_test['source']!=source, dtype=int)

    #-----------------------------------
    #Fit models

    #quantization matrices
    qm_model = quant_matrices.attribution_quant_matrices()
    qm_model.fit(df_train[['source', 'q_name']], compr_category_names=all_q_name_vals)

    #-----------------------------------
    #prediction on test set
    claimed_source_list = [source]*len(df_test)

    LLRs_isfake_qm, probs_fromsource_qm, probs_notfromsource_qm, \
      unrecognized_sources_qm = qm_model.predict(df_test['q_name'], claimed_source_list)
    df_test['LLR_qm'] = LLRs_isfake_qm

    #Determine if prediction is wrong
    misclassified_qm = (df_test['is_misattributed'] - .5) * LLRs_isfake_qm < 0
    df_test['misclassified_qm'] = misclassified_qm

    #-----------------------------------
    #Use hypothesis test score to compute ROC curve for this source:

    numsamples_balanced_testset.append(len(df_test))

    fpr, tpr, thresholds = roc_curve(df_test['is_misattributed'], df_test['LLR_qm'], pos_label=1)
    roc_auc = auc(fpr, tpr)
    AUCs_qm.append(roc_auc)

    results_qm[source] = {'source': source, 'fpr': fpr, 'tpr':tpr, 
                       'auc':roc_auc, 'numsamples':len(df_test),
                       'scores_isfake': df_test['LLR_qm'],
                       'label_isfake': df_test['is_misattributed'],
                       'df_test':df_test}
    results_per_trial_qm[trial] = results_qm
# -



# ## Summarize results

# +
FPR_THRESHOLD = 0.005  # compute TPR @ this FPR = 0.5%

AUCs_mean_cl = []
AUCs_std_cl = []
tpr_at_fpr_mean_cl = []

AUCs_mean_qm = []
AUCs_std_qm = []
tpr_at_fpr_mean_qm = []


#quantization matrices qm
for source in sources:
  AUCs_per_trial = []
  tpr_per_trial = []
  fpr_per_trial = []
  tprs_at_fpr_threshold = []
  for trial in range(NUM_TRIALS):
    AUCs_per_trial.append(results_per_trial_qm[trial][source]['auc'])
    fpr = results_per_trial_qm[trial][source]['fpr']
    tpr = results_per_trial_qm[trial][source]['tpr']
    fpr_per_trial.append(fpr)
    tpr_per_trial.append(tpr)
    tprs_at_fpr_threshold.append( np.interp(FPR_THRESHOLD, fpr, tpr) )
  AUCs_mean_qm.append(np.mean(AUCs_per_trial))
  AUCs_std_qm.append(np.std(AUCs_per_trial))
  tpr_at_fpr_mean_qm.append(np.mean(tprs_at_fpr_threshold))
  
df_summary = pd.DataFrame({'source':sources, 'test_size':numsamples_balanced_testset,
  'AUC_mean_qm':AUCs_mean_qm, 'AUC_std_qm':AUCs_std_qm, 
  'tpr_at_fpr_mean_qm':tpr_at_fpr_mean_qm,                         
  } )

df_summary['AUC_formatted_qm'] = uncertainty_format_arrays(df_summary['AUC_mean_qm'], df_summary['AUC_std_qm'])

df_summary

# -



# # Plot multiple ROC curves on one graph with uncertainty bands

# +
EERs_mean_qm = []

#------------
#New-EER
EERs_all_qm = []
EERs_std_qm = []
#------------

plt.figure(figsize=(6,5))
plt.plot([0, 1], [0, 1], color="black", linestyle="--")
plt.plot(np.linspace(0,1,100), 1-np.linspace(0,1,100), color="red", linestyle="--")
interp_fpr = np.linspace(0, 1, 1000)
for source in sources[0:15]:
  #New-EER
  EERs_per_src_qm = []
  
  interp_tprs = []
  
  #interpolate between fpr,tpr datapoints to compute tpr at regular fpr intervals
  for trial in range(NUM_TRIALS):
    fpr = results_per_trial_qm[trial][source]['fpr']
    tpr = results_per_trial_qm[trial][source]['tpr']
    interp_tpr = np.interp(interp_fpr, fpr, tpr)
    interp_tpr[0] = 0.0
    interp_tprs.append(interp_tpr)
  #------------
    #New-EER
    EERs_per_src_qm.append(calculate_eer(interp_fpr, interp_tpr)) #get EERs across all trials for this source
  EERs_std_qm.append( np.std(EERs_per_src_qm) ) #gives a std of EER for each source, across all 5 trials
  EERs_all_qm.append(EERs_per_src_qm) #all data: first index gives src, second index gives trial
  #------------

  mean_tpr = np.mean(interp_tprs, axis=0)
  mean_tpr[-1] = 1.0
  
  EERs_mean_qm.append(calculate_eer(interp_fpr, mean_tpr))

  std_tpr = np.std(interp_tprs, axis=0)
  tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
  tprs_lower = np.maximum(mean_tpr - std_tpr, 0)

  plt.plot(interp_fpr, mean_tpr,  linestyle="-", label=source)

  plt.fill_between(interp_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2)

  auc_mean = float(df_summary.loc[df_summary['source']==source, 'AUC_mean_qm'])
  auc_std = float(df_summary.loc[ df_summary['source']==source, 'AUC_std_qm'])
  tpr_at_fpr_mean = float(df_summary.loc[ df_summary['source']==source, 'tpr_at_fpr_mean_qm'])
  numsamples = int(df_summary.loc[ df_summary['source']==source, 'test_size'])
plt.xlabel("FPR")
plt.ylabel("TPR")
plt.legend()
plt.title("Verification: Time shift (part 1)")
plt.tight_layout()
plt.show()
#uncomment to save:
plt.savefig(os.path.join(RESULTS_FOLDER,"roc_curves_all_curves1.pdf"), bbox_inches='tight')

plt.figure(figsize=(6,5))
plt.plot([0, 1], [0, 1], color="black", linestyle="--")
plt.plot(np.linspace(0,1,100), 1-np.linspace(0,1,100), color="red", linestyle="--")
interp_fpr = np.linspace(0, 1, 1000)
for source in sources[15:]:
  
  #New-EER
  EERs_per_src_qm = []
  
  interp_tprs = []

  #interpolate between fpr,tpr datapoints to compute tpr at regular fpr intervals
  for trial in range(NUM_TRIALS):
    fpr = results_per_trial_qm[trial][source]['fpr']
    tpr = results_per_trial_qm[trial][source]['tpr']
    interp_tpr = np.interp(interp_fpr, fpr, tpr)
    interp_tpr[0] = 0.0
    interp_tprs.append(interp_tpr)
  #------------
    #New-EER
    EERs_per_src_qm.append(calculate_eer(interp_fpr, interp_tpr)) #get EERs across all trials for this source
  EERs_std_qm.append( np.std(EERs_per_src_qm) ) #gives a std of EER for each source, across all 5 trials
  EERs_all_qm.append(EERs_per_src_qm) #all data: first index gives src, second index gives trial
  #------------

  mean_tpr = np.mean(interp_tprs, axis=0)
  mean_tpr[-1] = 1.0
  
  EERs_mean_qm.append(calculate_eer(interp_fpr, mean_tpr))

  std_tpr = np.std(interp_tprs, axis=0)
  tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
  tprs_lower = np.maximum(mean_tpr - std_tpr, 0)

  plt.plot(interp_fpr, mean_tpr,  linestyle="-", label=source)

  plt.fill_between(interp_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2)

  auc_mean = float(df_summary.loc[df_summary['source']==source, 'AUC_mean_qm'])
  auc_std = float(df_summary.loc[ df_summary['source']==source, 'AUC_std_qm'])
  tpr_at_fpr_mean = float(df_summary.loc[ df_summary['source']==source, 'tpr_at_fpr_mean_qm'])
  numsamples = int(df_summary.loc[ df_summary['source']==source, 'test_size'])
plt.xlabel("FPR")
plt.ylabel("TPR")
plt.legend()
plt.title("Verification: Time shift (part 2)")
plt.tight_layout()
plt.show()
#uncomment to save:
plt.savefig(os.path.join(RESULTS_FOLDER,"roc_curves_all_curves2.pdf"), bbox_inches='tight')


df_summary['EER_mean_qm'] = EERs_mean_qm

#New-EER
df_summary['EER_std_qm'] = EERs_std_qm
# -



# +

df_latex = df_summary[['source', 'test_size', 'AUC_formatted_qm',  'tpr_at_fpr_mean_qm',  'EER_mean_qm']]
df_latex.columns=['source', 'test size', 'AUC', 'tpr@fpr', 'EER']
df_latex
# -

# 3 sig figs use '%.3g'; 3 digits use '%.3f'
latex_table = df_latex.to_latex(index=False, float_format='%.3g')
with open(os.path.join(RESULTS_FOLDER,"table1.tex"),"w") as file1:
  file1.write(latex_table)









# +
df_metricplot =  df_summary.sort_values(by='AUC_mean_qm', ascending=False).reset_index(drop=True)
sources_metricplot = list(df_metricplot['source'])
AUCs_metricplot = list(df_metricplot['AUC_mean_qm'])

# plt.figure(figsize=(6,3.5))
plt.figure(figsize=(6,2.6))
x_vals = [i for i,_ in enumerate(sources_metricplot)]

# plt.plot(x_vals, df_metricplot['EER_mean_qm'], linestyle='--', marker='.', label="EER", color="tab:blue")
plt.errorbar(x_vals, df_metricplot['EER_mean_qm'], yerr=df_metricplot['EER_std_qm'], fmt=".", linestyle="--",
             label="EER", color="tab:blue", mfc="tab:blue", mec='tab:blue', ecolor="tab:blue", capsize=2)


plt.errorbar(x_vals, AUCs_metricplot, yerr=df_metricplot['AUC_std_qm'], fmt=".", linestyle="--", 
             label="AUC", color="tab:orange", mfc="tab:orange", mec='tab:orange', ecolor="tab:orange", capsize=2)


plt.xticks(x_vals, sources_metricplot, rotation=90)

handles, labels = plt.gca().get_legend_handles_labels()
handles = [handles[1], handles[0]]
labels = [labels[1], labels[0]]
plt.legend(handles, labels, loc="center left")

plt.title("Verification metrics: Time generalization")

# plt.tight_layout()
plt.yticks(np.arange(0.0, 1.2, 0.2))
plt.show()
#uncomment to save:
plt.savefig(os.path.join(RESULTS_FOLDER,"verification_metrics_plot.pdf"), bbox_inches='tight')
# -



df_summary[['source', 'AUC_mean_qm']]













# ### Save results used to make the figure, so I can add the curves to plots made in exp 01 code (to combine figs from experiments 1 and 2)

df_metricplot.to_csv(os.path.join(RESULTS_FOLDER,"exp_02_metrics_plot_data.csv"), index=False)



StarcoderdataJupyterScriptsDedupFiltered