# ---
# 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
# ---

# +
import pandas as pd, numpy as np, seaborn as sns, geopandas as gpd, matplotlib.pyplot as plt
from scipy.stats import f_oneway
import pingouin as pin
import matplotlib

from shapely.geometry import Point, mapping
# %matplotlib inline

import warnings
warnings.simplefilter(action="ignore")

pd.set_option("display.precision", 2)
pd.options.display.float_format = '{:20.2f}'.format
pd.set_option('display.float_format', lambda x: '%.3f' % x)
pd.options.mode.chained_assignment = None

import cityImage as ci, ABManalysis as af

#libraries for clustering
from sklearn.cluster import KMeans, AgglomerativeClustering
from sklearn.metrics import silhouette_samples, silhouette_score
# -

# ### Preparing and cleaning the routes

# ### Coordinate System of the case study area for cartographic visualisations

city_name = 'Muenster'
epsg = 25832
crs = 'EPSG:'+str(epsg)

# # Questionnaire
# ## 1. Pre-processing

raw = pd.read_csv("Input/empiricalABM/Muenster_responses_matrix.csv") # raw responses
raw['startdate'] = raw.apply(lambda row: row['startdate'][11:], axis = 1)
raw['datestamp'] = raw.apply(lambda row: row['datestamp'][11:], axis = 1)
raw['duration'] = raw.apply(lambda row: af.compute_duration(row['startdate'], row['datestamp']), axis = 1)    
raw['PD2'] = raw['PD2'].fillna(raw['PD2'].mean())
raw['age'] = raw['PD2'].astype(int)
print('Number of total participants:', str(len(raw)))

# ### 1.1 Disregard subjects who took less than 20 minutes to complete the study

# cleaning
raw = raw[raw['duration']>= 20].copy()
print('Number of participants whose record was not disregarded ', str(len(raw)))
limit = {'A1':18, 'A2':25, 'A3':33, 'A4':41, 'A5':49, 'A6':57, 'A7':65, 'A8':73, 'A9':150}
raw['limit'] = raw[~raw.FQ1.isnull()].apply(lambda row: limit[row['FQ1']], axis = 1)
# print("age issue", len(raw[~raw.FQ1.isnull()][['PD2', 'FQ1', 'limit']][raw.age > raw.limit]) + 
#       len(raw[~raw.FQ1.isnull()][['PD2', 'FQ1', 'limit']][raw.age < raw.limit-7]))

raw.rename(columns={"PD1": "sex"}, inplace = True)
values = ['A1', 'A2', 'A3', 'A4']
new_values = ["male", "female", "non-binary", "prefer not to"]
for n, value in enumerate(values):
    raw['sex'].replace(value, new_values[n], inplace = True)


# ### 1.1 Disregarding subjects who always chose the same response when asked how to proceed along the route

# +
def check(index):
    response = list(raw[to_check].loc[index][raw[to_check].loc[index].notna()])
    percentageA1 = response.count("A1")/len(response)
    percentageA2 = response.count("A2")/len(response)
    if ((percentageA1 > 0.65) | (percentageA2 > 0.65)) & (raw.loc[index].duration < 30):
        return True
    return False

# checking videos
video_columns = [col for col in raw if col.startswith('VD')]
to_remove = ['VD000a','VD000','VD099[SQ001]', 'VD100a', 'VD100', 'VD199[SQ001]','VD200a','VD200', 'VD299[SQ001]']
to_check = [item for item in video_columns if item not in to_remove]
raw['allSame'] = raw.apply(lambda row: check(row.name), axis = 1)
# -

# ## 2. Demographic information

# ### 2.1 General information

# +
nr = 'Number'
pr = 'participants'

N = len(raw)
print(nr+' of female '+pr +': '+str(len(raw[raw.sex == 'female']))+', '+str(round(len(raw[raw.sex == 'female'])/N*100,1))+'%')
print(nr+' of male '+pr +': '+str( len(raw[raw.sex == 'male']))+', '+str(round(len(raw[raw.sex == 'male'])/N*100,1))+'%')
print(nr+' of non-binary '+pr +': '+str(len(raw[raw.sex == 'non-binary']))+', '+str(round(len(raw[raw.sex == 'non-binary'])/N*100, 1))+'%')
print(nr+' of '+pr+' who preferred not to declare their gender: '+str(len(raw[raw.sex == 'prefer not to']))+', '+
      str(round(len(raw[raw.sex == 'prefer not to'])/N*100, 1))+'%')
print()
print("Particpants' mean age:", round(raw['age'].mean(),3))
print("Particpants' std age:", round(raw['age'].std(), 3))
print()
print('Mean duration:', str(round(raw.duration.mean(), 2))+" minutes")
# -

# ### 2.2 Age categories

# +
Ga = len(raw[raw.age < 18])
Gb = len(raw[(raw.age >= 18) & (raw.age <= 25)])
Gc = len(raw[(raw.age >= 26) & (raw.age <= 33)])
Gd = len(raw[(raw.age >= 34) & (raw.age <= 41)])
Ge = len(raw[(raw.age >= 42) & (raw.age <= 49)])
Gf = len(raw[(raw.age >= 50) & (raw.age <= 57)])
Gg = len(raw[(raw.age >= 58) & (raw.age <= 65)])
Gh = len(raw[(raw.age >= 66) & (raw.age <= 73)])
Gi = len(raw[(raw.age >= 74)])

print("Number of participants per age group")
print()
print('< 18: '+str(Ga)+', '+str(round(Ga/N*100, 1))+ '%')
print('18 - 25: '+str(Gb)+', '+str(round(Gb/N*100, 1))+ '%')
print('26 - 33: '+str(Gc)+', '+str(round(Gc/N*100, 1))+ '%')
print('34 - 41: '+str(Gd)+', '+str(round(Gd/N*100, 1))+ '%')
print('42 - 49: '+str(Ge)+', '+str(round(Ge/N*100, 1))+ '%')
print('50 - 57: '+str(Gf)+', '+str(round(Gf/N*100, 1))+ '%')
print('58 - 65: '+str(Gg)+', '+str(round(Gg/N*100, 1))+ '%')
print('66 - 73: '+str(Gh)+', '+str(round(Gh/N*100, 1))+ '%')
print('> 74:'+str(Gi)+', '+str(round(Gi/N*100, 1))+ '%')
# -

# ### 2.3 Links or relationship to the case-study area

# +
Ga = len(raw[~raw['PD3[SQ001]'].isnull()])
Gb = len(raw[~raw['PD3[SQ002]'].isnull()])
Gc = len(raw[~raw['PD3[SQ003]'].isnull()])
Gd = len(raw[~raw['PD3[SQ004]'].isnull()])
Ge = len(raw[~raw['PD3[SQ005]'].isnull()])

nr = "Number of participants who "

print(nr+'are tourists: '+str(Ga)+', '+str(round(Ga/N*100,1))+ '%')
print(nr+'work here: '+str(Gb)+', '+str(round(Gb/N*100,1))+ '%')
print(nr+'live her: '+str(Gc)+', '+str(round(Gc/N*100,1))+ '%')
print(nr+'study here: '+str(Gd)+', '+str(round(Gd/N*100,1))+ '%')
print(nr+'are occasional visitors: '+str(Ge)+', '+str(round(Ge/N*100,1))+ '%')
# -

# ### 2.4 Reasons for walking

# +
columns = ['WB2[SQ001]', 'WB2[SQ002]','WB2[SQ003]','WB2[SQ004]','WB2[SQ005]']
labels = ['Commuting to/from work', 'Commuting to/from school/college/university', 'For social activities', 
          'For exercise or free-time activities', 'For other daily errands and commitments']
          
for column in columns: 
    raw[column].replace('A1', 1, inplace = True)
    raw[column].replace('A2', 2, inplace = True)
    raw[column].replace('A3', 3, inplace = True)
    raw[column].replace('A4', 4, inplace = True)
    raw[column].replace('A5', 5, inplace = True)

total_prob = 0.0
for n, column in enumerate(columns):
    prob = round(raw[column].mean(),3)
    total_prob += prob
    print(labels[n]+':', str(round(raw[column].mean(),3))) # value on the Likert scale
    print("percentage: ", str(round(raw[column].mean()*100/5,3))+'%')
    if n < len(columns)-1: 
        print()
# -

# #### *Reasons for Walking: Transformation into probabilities*

for n, column in enumerate(columns):
    actual_prob = round(raw[column].mean()/total_prob,2)
    print(str(labels[n]) + ': '+str(actual_prob))

# ## 3. Familiriaty, spatial knowledge and preferences
# ### 3.1 Time spent living in the case-study area (*familiarity*)

raw[['PD4[SQ001]', 'PD4[SQ002]']] = raw[['PD4[SQ001]', 'PD4[SQ002]']].fillna(0.0)
raw['PD4[SQ002]'].replace(1984.0, 2021-1984, inplace = True)
raw['familiarity'] = raw['PD4[SQ001]']/12 + raw['PD4[SQ002]']
raw['familiarity'].plot.kde()

# ### 3.2 Self-reported spatial knowledge

columns = ['WB4[SQ001]', 'WB4[SQ002]','WB4[SQ003]','WB4[SQ004]']
for column in columns: 
    raw[column].replace('A1', 1, inplace = True)
    raw[column].replace('A2', 2, inplace = True)
    raw[column].replace('A3', 3, inplace = True)
    raw[column].replace('A4', 4, inplace = True)
    raw[column].replace('A5', 5, inplace = True)
raw['knowledge'] = (raw['WB4[SQ001]'] + raw['WB4[SQ002]'] + raw['WB4[SQ003]'] +raw['WB4[SQ004]'])/4
raw['knowledge'].plot.kde()

columns = ['WB4[SQ001]', 'WB4[SQ002]','WB4[SQ003]','WB4[SQ004]', 'knowledge', 'familiarity']
raw[columns].corr(method='pearson')

# ### 3.3 Preference for and aversion to barriers

# +
columns = ['SP4[SQ002]', 'SP4[SQ006]','SP4[SQ001]','SP4[SQ005]']
for column in columns: 
    raw[column].replace('A1', 0.00, inplace = True)
    raw[column].replace('A2', 0.25, inplace = True)
    raw[column].replace('A3', 0.50, inplace = True)
    raw[column].replace('A4', 0.75, inplace = True)
    raw[column].replace('A5', 1.00, inplace = True)

raw['preferenceNatural'] = (raw['SP4[SQ002]'] + raw['SP4[SQ006]'])/2 ## preference for Natural Barriers
raw['aversionSevering'] = raw['SP4[SQ001]'] ## Aversion to Severing Barriers
# -

# ## 4. Video Tasks: route choice behaviour variables
# #### *Preliminary columns cleaning*

# +
video_columns = [col for col in raw if col.startswith('VD')]
to_remove = ['VD000a','VD000','VD099[SQ001]', 'VD100a', 'VD100', 'VD199[SQ001]','VD200a','VD200', 'VD299[SQ001]']

columns = ['id'] + video_columns
columns = [item for item in columns if item not in to_remove] # remove not necessary video columns
responses = raw[columns].copy() # only response to the video taks

# +
# cleaning and rename columns and prepare a legible dataframe 
for column in responses.columns:
    if column == 'id':
        continue
    responses.rename({column: column[2:]}, axis=1, inplace = True)
    
for column in responses.columns:
    if column == 'id':
        continue
    if ('a' in column) | ('b' in column): 
        continue
    new_column = '-'.join(column[i:i+3] for i in range(0, len(column), 3))
    responses.rename({column: new_column}, axis=1, inplace = True)
    
responses.rename({'150a151152': '150a-151-152'}, axis=1, inplace = True)
responses.rename({'150b126': '150b-126'}, axis=1, inplace = True)
# -

# #### *Loading the routes used in the survey* 

# +
survey_routes = gpd.read_file("Outputs/Routes_sections.shp").to_crs(crs)
summary = pd.DataFrame(columns = ['id']+list(survey_routes.video.unique()))
del summary[None]

summary = summary.reindex(sorted(summary.columns), axis=1)

for row in responses.itertuples():
    sectors = responses.loc[row.Index].notna().dot(responses.columns+',').rstrip(',')
    sectors = sectors.replace('-', ',')
    sectors_list = sectors.strip('').split(',')
    for sector in sectors_list:
        summary.at[row.Index, sector] = 1
        summary.at[row.Index, 'id'] = responses.loc[row.Index].id
        
summary.drop('114', inplace = True, axis = 1)
summary.fillna(0, inplace = True)
for column in summary.columns: 
    summary.rename(columns={column: str(column)}, inplace = True)
# -

summary.head()

# ### 4.1 Obtaining the general statistics

survey_routes['routeChoice'] = survey_routes.apply(lambda row: row['routeChoic'].replace(" ",""), axis = 1)
survey_routes['routeChoice'] = survey_routes.apply(lambda row: row['routeChoice'].strip('][').split(','), axis = 1)
survey_routes.drop('routeChoic', inplace = True, axis = 1)

# #### *Computing for each subject how much they resorted to certain urban elements or road costs*

# +
route_variables = ['usingElements', 'noElements','onlyDistance', 'onlyAngular', 'distanceHeuristic', 'angularHeuristic', 
                   'regions', 'routeMarks', 'barriers','distantLandmarks', 'preferenceNatural', 'aversionSevering']
other_variables = ["length", "minimisation_length","combined_length"]

video0 = [col for col in summary if col.startswith('0')]
video1 = [col for col in summary if col.startswith('1')]
video2 = [col for col in summary if col.startswith('2')]
videos = [video0, video1, video2]
routes_stats = af.set_routes_stats(summary, survey_routes, videos)
# -

raw.index = raw['id']
routes_stats.index.name = None
routes_stats['preferenceNatural'] = raw['preferenceNatural']
routes_stats['aversionSevering'] = raw['aversionSevering']
routes_stats['knowledge'] = raw['knowledge']
routes_stats.head(10)

## appending variables regarding the preference for natural/severing barriers and spatial knowledge
raw.index = raw['id']
routes_stats.index.name = None
routes_stats['preferenceNatural'] = raw['preferenceNatural']
routes_stats['aversionSevering'] = raw['aversionSevering']
routes_stats['knowledge'] = raw['knowledge']
routes_stats.head(10)

# #### *Obtaining the input matrix for cluster analysis*

# +
input_matrix = routes_stats[route_variables].copy()
input_matrix.replace(0.0, 0.01, inplace = True)
cluster_variables =  ['onlyDistance', 'onlyAngular', 'distanceHeuristic', 'angularHeuristic', 'regions', 
                      'routeMarks', 'barriers','distantLandmarks']

for column in cluster_variables:
    if column in ['onlyDistance', 'onlyAngular']:
        input_matrix[column] = input_matrix[column] * input_matrix['noElements'] 
    else:
        input_matrix[column] = input_matrix[column] * input_matrix['usingElements']

input_matrix = input_matrix.astype(float)
# -

# ### 4.2 Visualising overall importance of each urban element (probabilities)

# +
sns.set()
sns.set_color_codes()
fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize=(25, 12)) 

## obtaining dataframe for boxplot visualisation
tab = pd.DataFrame(columns = {'variable', 'value'}) 
labels = ['using\nUE', 'only\nGMH', 'distance\nGMH', 'angular\nGMH', 'distance\nLMH', 'angular\nLMH', 'regions', 
          'on-route\nmarks', 'barriers', 'distant\nlandmarks', 'natural\nbarriers', 'severing\nbarriers']

tab = pd.DataFrame(columns = {'variable', 'value'}) 
index = 0
for subject in input_matrix.index:
    for n, variable in enumerate(route_variables):
        tab.at[index, 'variable'] = labels[n]
        tab.at[index, 'value'] = input_matrix.loc[subject][variable]
        index += 1

palette = ['coral']*2+['gold']*8+['mediumseagreen']*2
ax = sns.boxplot(x="variable", y="value", data=tab, palette = palette)
ax.set_ylim(0.0, 1.01)
ax.set_ylabel('', fontsize = 27, labelpad = 30, fontfamily = 'Times New Roman')
ax.set_xlabel('', labelpad = 30, fontfamily = 'Times New Roman')

for tick in ax.get_yticklabels(): 
    tick.set_fontname('Times New Roman')
for tick in ax.get_xticklabels(): 
    tick.set_fontname('Times New Roman')
    
ax.tick_params(axis='both', labelsize= 27, pad = 20)
# -

fig.savefig("Outputs/Figures/empiricalABM/f4.pdf", bbox_inches='tight')

input_matrix[cluster_variables].mean()

input_matrix[cluster_variables].corr(method='pearson')

# ## 5. Cluster Analysis
# ### 5.1 Variables transformation

X = input_matrix[cluster_variables].copy()
fig = plt.figure(figsize = (20, 10))
for n, column in enumerate(list(X.columns)):
    ax = fig.add_subplot(2,4,n+1)
    ax = X[column].plot.kde()
    ax.set_title(column)

# #### *Logarithmic transformation and standardisation*

# +
X_log = af.log_transf(input_matrix[cluster_variables])
X_log_stand = X_log.copy()
    
for column in X_log.columns:
    X_log_stand[column] = af.standardise_column(X_log, column)
# -

# ### 5.2 Obtaining and comparing different clustered structures obtained with the k-means algorithm

# +
list_scores = []

def pipe_kmeans(X_toFit):

    X_tmp = X_toFit.copy()
    for n_clusters in range(2, 10):
        clusterer = KMeans(n_clusters = n_clusters, n_init = 2000).fit(X_toFit)
        labels = clusterer.labels_

        score = round(silhouette_score(X_tmp, labels, metric= 'sqeuclidean'), 3)
        param_score = {'algorithm': 'k-means', 'n_clusters' : n_clusters, 'score' : score,
                       'clusterer' : clusterer}
        list_scores.append(param_score)

pipe_kmeans(X_log_stand)
clustering = pd.DataFrame(list_scores)
clustering.sort_values(by = 'score', ascending = False)


# +
def pipe_kmeans_VCR(X):

    VRCs = []
    X_tmp = X.copy()
    X_toFit = X.copy()
    
    for n_clusters in range(2, 11):
        
        clusterer = KMeans(n_clusters = n_clusters, n_init = 2000).fit(X_toFit)
        X_tmp['cluster'] = clusterer.labels_

        f = 0.0
        for var in cluster_variables:
            arrays = []
            for cluster in list(X_tmp['cluster'].unique()):
                arrays.append(X_tmp[var+'_log'][X_tmp.cluster == cluster])
            
            if n_clusters == 2:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1])).statistic
            if n_clusters == 3:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2])).statistic
            if n_clusters == 4:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), 
                                     np.array(arrays[3])).statistic
            if n_clusters == 5:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), np.array(arrays[3]), 
                          np.array(arrays[4])).statistic
            if n_clusters == 6:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), np.array(arrays[3]), 
                          np.array(arrays[4]),np.array(arrays[5])).statistic
            if n_clusters == 7:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), np.array(arrays[3]), 
                          np.array(arrays[4]),np.array(arrays[5]),np.array(arrays[6])).statistic
            if n_clusters == 8:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), np.array(arrays[3]), 
                          np.array(arrays[4]), np.array(arrays[5]), np.array(arrays[6]), np.array(arrays[7])).statistic
            if n_clusters == 9:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), np.array(arrays[3]), 
                          np.array(arrays[4]),np.array(arrays[5]), np.array(arrays[6]), np.array(arrays[7]), np.array(arrays[8])
                                    ).statistic
            if n_clusters == 10:
                statistic = f_oneway(np.array(arrays[0]), np.array(arrays[1]), np.array(arrays[2]), np.array(arrays[3]), 
                          np.array(arrays[4]), np.array(arrays[5]), np.array(arrays[6]), np.array(arrays[7]), 
                                     np.array(arrays[8]), np.array(arrays[9])).statistic
            
            f += statistic
        VRCs.append(f)   
    return VRCs

VRCs = pipe_kmeans_VCR(X_log_stand)
# -

# ### 5.3 Choosing the best structure

for n, VRC in enumerate(VRCs):
    if (n == 0) | (n == len(VRCs)-1):
        continue
    index = (VRCs[n+1]-VRC)-(VRC-VRCs[n-1])
    print("structure with", n+2, "clusters: Omega is", round(index, 3), "VRC is", round(VRC,3))

chosen = clustering.loc[4]
labels_cluster = np.array(chosen.clusterer.labels_.copy())
labels_cluster += 1

input_with_cluster = input_matrix.copy()
input_with_cluster['cluster'] = labels_cluster

# ### 5.4 Examining the chosen partition

input_with_cluster['knowledge'] = raw['knowledge']
cluster_stats = input_with_cluster.groupby("cluster").mean()
cluster_stats

# +
sns.set()
sns.set_color_codes()

colors = ['royalblue', 'gold', 'orchid', 'tomato', 'lightgreen', 'sandybrown', 'sienna', 'darkgreen',
         'c', 'black']

clusters_stats = input_with_cluster.groupby(by = 'cluster').mean()
figsize = (22.5, (15/2*3))
fig = plt.figure(figsize = figsize)
plot_variables = ['onlyDistance',  'regions', 'routeMarks','angularHeuristic', 'onlyAngular', 'barriers', 
                  'distantLandmarks', 'distanceHeuristic']
cluster_v = np.array(plot_variables)
labels = ['distance GMH', 'regions','on-route marks','angular LMH','angular GMH','barriers','distant landmarks',
          'distance LMH'] 
plt.rcParams['font.family'] = 'Times New Roman'
to_plot = list(clusters_stats.index)

for cluster in to_plot:

    if cluster != 'population':
        tmp = clusters_stats.loc[cluster, cluster_v].values

    angles = np.linspace(0, 2*np.pi, len(labels), endpoint=False)
    
    # close the plot
    tmp = np.concatenate((tmp,[tmp[0]]))
    angles = np.concatenate((angles,[angles[0]]))
    
    if cluster != 'population':
        ax = fig.add_subplot(3, 3, cluster, polar = True)
        color = colors[cluster-1]
    
    ax.plot(angles, tmp, '-', color = color, linewidth=2)
    ax.fill(angles, tmp, color = color, alpha=0.25)
    ax.set_thetagrids((angles * 180/np.pi)[0:len(plot_variables)], labels, fontsize = 15)
    ax.yaxis.set_ticks([0.10, 0.20, 0.30, 0.40, 0.50])
    ax.tick_params(axis='y', labelsize= 15)
    ax.tick_params(axis='x', labelsize= 15, pad = 5)
    ax.set_rlabel_position(150)
    
    if cluster != 'population':
        nr = len(input_with_cluster[input_with_cluster.cluster == cluster])
        title = "Cluster "+str(cluster) + " (N = "+str(nr)+")"      
    ax.set_title(title, va = 'bottom', fontsize = 20, pad = 50, fontfamily = 'Times New Roman')

fig.subplots_adjust(wspace=0.20, hspace = 0.60)
# -

fig.savefig("Outputs/Figures/empiricalABM/f5.pdf", bbox_inches='tight')

# #### *Spatial skills - Do they explain the variation in the route choice behaviour variables?*

clusters = list(set(labels_cluster))
for cluster in clusters:
    knowledge = input_with_cluster[input_with_cluster.cluster == cluster]['knowledge']
    print("cluster", cluster, "spatial skills/knowledge mean:", round(knowledge.mean(), 3), "std:", round(knowledge.std(), 3))

# +
tab = pd.DataFrame(columns = {'cluster', 'variable', 'value'}) 
labels = ['distance min.', 'angular min.', 'distance LMH', 'angular LMH', 'regions', 
          'on-route marks', 'barriers', 'distant landmarks', 'knowledge']

tab = pd.DataFrame(columns = {'cluster', 'variable', 'value'}) 
index = 0
for subject in input_with_cluster.index:
    for n, variable in enumerate(cluster_variables+['knowledge']):
        tab.at[index, 'variable'] = labels[n]
        tab.at[index, 'cluster'] = "Cluster "+str(int(input_with_cluster.loc[subject]['cluster']))
        tab.at[index, 'value'] = input_with_cluster.loc[subject][variable]
        index += 1

tab['value'] = tab.value.astype(float)
tab.sort_values('cluster', inplace = True)
# -

## using the t-test to verify whether differences in the spatial skills/knowledge are significant
pin.pairwise_ttests(data = tab[tab.variable == 'knowledge'], dv= 'value', between = 'cluster')

# #### *Demographic characteristics*

raw['cluster'] = labels_cluster
print('Age')
cluster_deom_stats = raw.groupby("cluster")['age'].mean()
print(cluster_deom_stats)

print('Gender')
cluster_deom_stats = raw.groupby(["cluster", 'sex'])['sex'].count()
for cluster in range(1, max(labels_cluster)+1):
    print("cluster", cluster)
    for gender in list(cluster_deom_stats.loc[cluster].index):
        print(gender, round(cluster_deom_stats.loc[cluster][gender]/cluster_deom_stats.loc[cluster].sum() *100, 1))
    print()

# ## 6. Final input data for building agent typologies in the ABM

# +
columns = []
for column in input_with_cluster.columns:
    if column in ['knowledge', 'cluster']:
        continue
    columns.append(column+"_mean")
    columns.append(column+"_std")
    
groups = ['group'+str(cluster) for cluster in input_with_cluster['cluster'].unique()]
indexes = groups + ['population', 'nullGroup']
clusters_gdf = pd.DataFrame(index = indexes, columns = columns)
variables = route_variables


clusters = input_with_cluster['cluster'].unique()
for index in indexes:
    for variable in variables:
        if index not in ['population', 'nullGroup']:
            cluster = int(index[5:])
            clusters_gdf.at[index, variable+'_mean'] = input_with_cluster[input_with_cluster.cluster == cluster][variable].mean()
            clusters_gdf.at[index, variable+'_std'] = input_with_cluster[input_with_cluster.cluster == cluster][variable].std()
        else:
            clusters_gdf.at[index, variable+'_mean'] = input_with_cluster[variable].mean()
            clusters_gdf.at[index, variable+'_std'] = input_with_cluster[variable].std()
    
    if index not in ['population', 'nullGroup']:
        cluster = int(index[5:])
        clusters_gdf.at[index, 'portion'] = len(input_with_cluster[input_with_cluster.cluster == cluster])/len(input_with_cluster)
    else:
        clusters_gdf.at[index, 'portion'] = 1.00
# -

clusters_gdf

## exporting the files 
clusters_gdf.to_csv('Outputs/empiricalABM/'+city_name+'_clusters.csv') ## this is imported into the ABM
routes_stats.to_csv('Outputs/empiricalABM/'+city_name+'_routes_stats.csv')

StarcoderdataJupyterScriptsDedupFiltered