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

# +
import numpy as np
import json
from scipy import optimize
from os import listdir
import re
import os
from os.path import isfile, join

import pandas as pd
from joblib import Parallel, delayed, dump
# %matplotlib inline
import matplotlib.pyplot as plt
from skimage import measure
import matplotlib
import seaborn as sns
import csv
from scipy.stats import stats
# -

NJOBS = 35


# +
#sns.set(style="white")
#sns.set_context("paper")

def figsize(scale,ratio):
    fig_width_pt = 468.0                      # Get this from LaTeX using \the\textwidth
    inches_per_pt = 1.0/72.27                       # Convert pt to inch
    golden_mean = (np.sqrt(5.0)-1.0)/2.0            # Aesthetic ratio (you could change this)
    fig_width = fig_width_pt*inches_per_pt*scale    # width in inches
    if(ratio == "golden"):                          # Golden ratio ...
        fig_height = fig_width*golden_mean          # height in inches
    else:                                           # ... or other ratio
        fig_height = fig_width*ratio
        
    fig_size = [fig_width,fig_height]
    return fig_size

pgf_with_latex = {'backend': 'pdf',
                                            'axes.labelsize': 8,
                                            'xtick.labelsize': 8,
                                            'ytick.labelsize': 8,
                                            'legend.fontsize': 8,
                                            'lines.markersize': 3,
                                            'font.size': 8,
                                            'font.family': u'sans-serif',
                                            'font.sans-serif': ['Arial'],
                                            'pdf.fonttype': 42,
                                              'ps.fonttype': 42,
                                            'text.usetex': False}
matplotlib.rcParams.update(pgf_with_latex)

# +
from collections import defaultdict
import gzip
from tqdm import tqdm_notebook
from itertools import chain

def collect_distance_file(prefix, directoryname, tileid):
    with gzip.open(directoryname + '{}{}.json.gz'.format(prefix, tileid), 'rt') as f:
        d = json.load(f)
    return d['distances'][0]

def collect_distances(directoryname, df_classes, distance_name='energy', modelname=''):
    dours = {}
    drybski = {}
    prefix_name = '{}_'.format(distance_name)
    if modelname:
        prefix_name = '{}_{}_'.format(distance_name, modelname)
        
    prefix_name_rybski = '{}_{}_'.format(distance_name, 'rybski')
      
    for c in [1,2,3,4]:
        tiles = df_classes[df_classes['class'] == c]['tileid'].values
        dours[c] = [r for r in Parallel(n_jobs=NJOBS)(delayed(collect_distance_file)(prefix_name, directoryname, t) for t in tqdm_notebook(tiles)  if isfile(directoryname + '{}{}.json.gz'.format(prefix_name, t)))]
        
        drybski[c] = [r for r in Parallel(n_jobs=NJOBS)(delayed(collect_distance_file)(prefix_name_rybski, directoryname, t) for t in tqdm_notebook(tiles) if isfile(directoryname + '{}{}.json.gz'.format(prefix_name, t)))]
        
    return dours, drybski


def create_dataframe_distances(dours, drybski):
    data = [[ 10, 'tom', 'c']]
    df = pd.DataFrame(data, columns = ['dist', 'Class', 'model']) 

    for c in [1,2,3,4]:
        array_distances = np.hstack((np.array(dours[c]), np.array(drybski[c]))).reshape(-1, 1)


        repeated = np.repeat(c, len(array_distances)).reshape(-1, 1)
        repeated_model = np.hstack((np.repeat('multi', len(dours[c])), np.repeat('rybski', len(dours[c])))).reshape(-1, 1)

        df = pd.concat((df, pd.DataFrame(np.hstack((array_distances, repeated, repeated_model)), columns = ['dist', 'Class', 'model']) ))

    df = df[df['Class']  != 'tom']
    df['dist'] = df['dist'].astype('float32')
    df = df.sort_values(['model', 'Class'])
    return df

def compute_kolmogorov_classes(df):
    for c in ['1', '2', '3', '4']:
        data1 = df[(df['Class'] == c) & (df['model'] == 'multi')]['dist'].values
        data2 = df[(df['Class'] == c) & (df['model'] == 'rybski')]['dist'].values
        d, p = stats.ks_2samp(data1, data2)
        print("KL class {}: {} (p-value {})".format(c, d, p))
    


# -

df_classes = pd.read_csv('../data/generated_files/quantiles_classes.csv', dtype={'tileid': str})
df_classes.head()

df = pd.read_csv('../data/generated_files/summary_tiles_05x05.csv', dtype={'tileid': 'str'})
df['perc_urban'] = df['urban_area_km2'] / df['original_km2']
df['perc_constraint'] = 1 - (df['tile_km2'] / df['original_km2'])

df = pd.merge(df[['tileid', 'perc_urban', 'perc_constraint']], df_classes, on='tileid')

# +
rootdir = "../data/generated_files/simulations/distances/"

dours, drybski = collect_distances(rootdir, df_classes, distance_name='energy')
# -

df = create_dataframe_distances(dours, drybski)
df.head()

print("Model multi")
data1 = df[(df['model'] == 'multi')]['dist'].values
data2 = df[(df['model'] == 'rybski')]['dist'].values
d, p = stats.ks_2samp(data1, data2)
print("KL all: {} (p-value {})".format(d, p))

compute_kolmogorov_classes(df)

# +
colors = np.array([[109, 110, 113],
                   [230, 231, 232]])/255

f, axs = plt.subplots(1,1, figsize=figsize(1, 0.8))
g = sns.boxplot(x="Class", y="dist", hue='model', data=df, linewidth=0.5, showfliers=False,
           palette=sns.color_palette(colors))

leg = g.axes.get_legend()
leg.set_title('Model')
renames = {'multi': 'Multi parameter', 'rybski': 'Single parameter'}
for t in leg.texts: 
    old = t.get_text()
    t.set_text(renames[old])
    
    
original = np.median(df[df['model'] == 'rybski']['dist'].values)
decrease = original - np.median(df[df['model'] == 'multi']['dist'].values)

print("Median diff:", decrease/original*100)

plt.savefig('../figures/manuscript/boxplot_distances.pdf', dpi=150, bbox_inches='tight', pad_inches=0.05)
# -

doursKL, drybskiKL = collect_distances(rootdir, df_classes, distance_name='kl')
dfKL = create_dataframe_distances(doursKL, drybskiKL)
dfKL.head()

doursJS, drybskiJS = collect_distances(rootdir, df_classes, distance_name='js')
dfJS = create_dataframe_distances(doursJS, drybskiJS)
dfJS.head()


# +
def edit_distance_legend(g):
    leg = g.axes.get_legend()
    leg.set_title('Model')
    renames = {'multi': 'Multi parameter', 'rybski': 'Single parameter'}
    for t in leg.texts: 
        old = t.get_text()
        t.set_text(renames[old])

colors = np.array([[109, 110, 113],
                   [230, 231, 232]])/255

f, axs = plt.subplots(1,3, figsize=figsize(1.5, 0.4))
g = sns.boxplot(x="Class", y="dist", hue='model', data=df, linewidth=0.5, showfliers=False,
           palette=sns.color_palette(colors), ax=axs[0])
edit_distance_legend(g)
    
g = sns.boxplot(x="Class", y="dist", hue='model', data=dfKL, linewidth=0.5, showfliers=False,
           palette=sns.color_palette(colors), ax=axs[1])
edit_distance_legend(g)

g = sns.boxplot(x="Class", y="dist", hue='model', data=dfJS, linewidth=0.5, showfliers=False,
           palette=sns.color_palette(colors), ax=axs[2])
edit_distance_legend(g)

#axs[0].set_ylim((0.2, 2.2))

axs[0].set_ylabel("Distance")
axs[1].set_ylabel("Distance")
axs[2].set_ylabel("Distance")

axs[0].set_title("Energy distance")
axs[1].set_title("Kullback Leibler divergence")
axs[2].set_title("Jensen–Shannon divergence")

f.tight_layout()
plt.savefig('../figures/manuscript/appendix_boxplot_distances.pdf', dpi=150, bbox_inches='tight', pad_inches=0.05)
# -

dmultiprob, drybski = collect_distances(rootdir, df_classes, distance_name='energy', modelname='multi')

df_multiprob = create_dataframe_distances(dmultiprob, drybski)
df_multiprob.head()

print("Model multiprob")
data1 = df_multiprob[(df_multiprob['model'] == 'multi')]['dist'].values
data2 = df_multiprob[(df_multiprob['model'] == 'rybski')]['dist'].values
d, p = stats.ks_2samp(data1, data2)
print("KL all: {} (p-value {})".format(d, p))

compute_kolmogorov_classes(df_multiprob)

# ## Classes real vs simulated

df_classes = pd.read_csv('../data/generated_files/quantiles_classes.csv', dtype={'tileid': str})
print(len(df_classes))
df_classes.head()

# +
from functools import reduce

dfs = []
for m in ['1000', 'rybski', 'marco']:#
    df_temp = pd.read_csv('../data/generated_files/quantiles_classes_{}_all.csv'.format(m), dtype={'tileid': 'str'})
    df_temp = df_temp.rename(columns={'class': 'class_{}'.format(m), 'quantile': 'quantile_{}'.format(m)})
    df_temp['class_{}'.format(m)] = df_temp['class_{}'.format(m)].astype('int')
    dfs.append(df_temp)
    
df_merged_classes = reduce(lambda left,right: pd.merge(left,right, on='tileid'), dfs)
df_merged_classes = pd.merge(df_merged_classes, df_classes.rename(columns={'class': 'class_truth', 'quantile': 'quantile_truth'}), on='tileid')
print(len(df_merged_classes))
df_merged_classes.head()
# -

df_merged_classes.to_csv('../data/generated_files/quantiles_classes_merged_all.csv', index=False)

df_melted_classes = df_merged_classes.melt(id_vars=['tileid'], value_vars=['class_1000', 'class_rybski', 'class_truth'], #, 'class_marco'
                      var_name='model', value_name='class')
df_melted_classes.loc[:, 'model'] = df_melted_classes['model'].str.replace('class_', '')
df_melted_classes.head()

df_count_classes = df_melted_classes.groupby(['model', 'class'], as_index=False).count()
df_count_classes = df_count_classes.rename(columns={'tileid': 'num'})
df_count_classes.head()

# +
from  matplotlib.ticker import PercentFormatter

g = sns.catplot(x="class", y='num', hue="model", 
                data=df_count_classes[df_count_classes.model.isin({'rybski', '1000'})].sort_values('model', ascending=False),
                height=5, kind="bar", palette="muted")
#plt.axhline(y=500, linestyle='--', color='gray')
plt.ylabel("Number of tiles")

plt.savefig('../figures/manuscript/boxplot_classes.pdf', dpi=150, bbox_inches='tight', pad_inches=0.05)
# -

len(pd.merge(df_merged_classes, df_classes, on='tileid'))

len(df_merged_classes), len(df_classes)

# +
from sklearn.metrics import f1_score
                         
def f1_score_classes(true_classes, sim_classes, class_num):
    n = len(sim_classes)
    y_true = np.zeros(n)
    y_pred = np.zeros(n)
    y_true[true_classes == class_num] = 1
    y_pred[sim_classes == class_num] = 1
    #y_pred[sim_classes != class_num] = 0
    
    return f1_score(y_true, y_pred)

results = []

df_merged_classes2 = df_merged_classes.copy()
#print(len(df_merged_classes2))
#df_merged_classes2 = df_merged_classes2[~df_merged_classes2.tileid.isin(df_summ[df_summ.perc_constraint > 0.6].tileid.values)]
#print(len(df_merged_classes2))

for m in ['1000', 'rybski', 'marco']: #
    
    print("MODEL", m)
    for class_num in [1,2,3,4]:
        df_sub = df_merged_classes2[(df_merged_classes2['class_truth'] == class_num)]
        true_classes = df_sub['class_truth'].values
        sim_classes = df_sub['class_{}'.format(m)].values
        print(class_num, f1_score_classes(true_classes, sim_classes, class_num))
        results.append((m, class_num, f1_score_classes(true_classes, sim_classes, class_num)))
    print("")

# -

for m in ['1000', 'rybski', 'marco']: #
    
    print("MODEL", m)
    df_sub = df_merged_classes2.copy()
    true_classes = df_sub['class_truth'].values
    sim_classes = df_sub['class_{}'.format(m)].values
    print(class_num, f1_score(true_classes, sim_classes, average='weighted'))
    print("")


df = pd.DataFrame(data=results, columns=['model', 'class', 'f1'])
df.head()

# +
original = np.median(df[df['model'] == 'rybski']['f1'].values)
decrease = original - np.median(df[df['model'] == '1000']['f1'].values)

print("Median diff:", decrease/original*100)

# +
flatui = [np.array([151, 152, 155])/255., np.array([230,231,232])/255.]
palette = sns.color_palette(flatui)

g = sns.catplot(x="class", y='f1', hue="model", data=df[df.model.isin({'1000', 'rybski'})],
                height=5, kind="bar", palette=palette)
plt.ylabel("F1-score")

plt.savefig('../figures/manuscript/boxplot_classes.pdf', dpi=150, bbox_inches='tight', pad_inches=0.05)

# +
g = sns.catplot(x="class", y='f1', hue="model", data=df,
                height=3, kind="bar", palette=palette, aspect=1.5)
plt.ylabel("F1-score")

plt.savefig('../figures/manuscript/multi_boxplot_classes.pdf', dpi=150, bbox_inches='tight', pad_inches=0.05)
# -



StarcoderdataJupyterScriptsDedupFiltered