# ---
# 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 pandas as pd
from safegraph_py_functions import safegraph_py_functions as sgpy
import os
from dotenv import load_dotenv, find_dotenv
from loguru import logger

# +
# find .env automagically by walking up directories until it's found
dotenv_path = find_dotenv()

# load up the entries as environment variables
load_dotenv(dotenv_path)
os.chdir(os.environ.get("ROOT_DIR"))
from src import DATA_DIR
raw_data_dir = DATA_DIR / 'raw'
# -

# if import is 0 it reads the data from the existing file
# otherwise reads in the raw data an makes a unified dataset
IMPORT = 0 


# Read in all patterns files in the monthly-patterns folder
def get_files():
    patterns_path = raw_data_dir / "monthly-patterns"
    files = []
    for f in patterns_path.glob("**/*.csv.gz"):
        files.append(f)
    return files


def filter_to_philly(df):
    # zip codes are read as integers rather than strings so we add leading zeros.
    # this is not strictly necessary since Philadelphia zipcodes don't have leading zeros.
    
    # Philadelphia selection
    # HK: adding leading zeros because some zipcodes in MA are 0191X.
    df['postal_code'] = df['postal_code'].apply(lambda x: ('00000'+str(x))[-5:])
    in_philly = df['postal_code'].astype(str).str.startswith("191")
    df = df.loc[in_philly]
    df = df[['safegraph_place_id','date_range_start','postal_code', 
             'raw_visit_counts', 'raw_visitor_counts']]
    
    return df


if IMPORT == 1:
    philly_patterns = []
    files = get_files()
    for i, f in enumerate(files):
        print(f)
        philly_patterns.append(filter_to_philly(pd.read_csv(f)))
    philly_patterns = pd.concat(philly_patterns)
    philly_patterns.to_csv(
    DATA_DIR / "processed" / "kmeans_patterns.csv.tar.gz", index=False
    )
else:
    philly_patterns = pd.read_csv(DATA_DIR / "processed" / "philly_patterns.csv.tar.gz", low_memory = False)

philly_patterns.columns

philly_patterns= philly_patterns[['safegraph_place_id','date_range_start', 'date_range_end',
             'raw_visit_counts', 'raw_visitor_counts', 'poi_cbg','top_category']]
philly_patterns['poi_cbg'] = philly_patterns['poi_cbg'].astype(int)
philly_patterns['date_range_start'] = pd.to_datetime(philly_patterns['date_range_start'])
philly_patterns['date_range_end'] = pd.to_datetime(philly_patterns['date_range_end'].apply(lambda x: x[:10]))

census_path = DATA_DIR / 'raw' / 'open-census-data'
files = [file for file in census_path.glob('**/cbg_b19.csv')]
census = pd.read_csv(files[0])
census =  census[['census_block_group','B19013e1']]

philly_patterns = philly_patterns.merge(census, left_on = 'poi_cbg', right_on = 'census_block_group', how = 'left')
philly_patterns.drop(columns = 'census_block_group', inplace = True)

philly_patterns['income_missing'] = philly_patterns['B19013e1'].isna().astype(int)
philly_patterns['B19013e1'].fillna(value=philly_patterns['B19013e1'].mean(), inplace=True)


# +
# for files with information disaggregated at the state level, keep only the country-wide info
def keep_total_level(norm_stats):
    if 'region' in norm_stats.columns:
        if len(norm_stats[norm_stats['region'] == 'ALL_STATES']) == 0:
            raise ValueError('no region named "ALL_STATES"')
        norm_stats = norm_stats[norm_stats['region'] == 'ALL_STATES']
        norm_stats = norm_stats.drop(columns = ['region'])
    return norm_stats


patterns_path = raw_data_dir / "monthly-patterns"

    
norm_files = [f for f in patterns_path.glob("**/normalization_stats.csv")]
#read in normalization data
norm_stats = pd.concat([keep_total_level(pd.read_csv(file)) for file in norm_files])
norm_stats['year'] = norm_stats['year'].astype(int)
norm_stats['month'] = norm_stats['month'].astype(int)
norm_stats['day'] = norm_stats['day'].astype(int)
# HK: I only downloaded patterns data from 2019 onwards due to memory constraints
norm_stats = norm_stats[norm_stats['year'] >= 2019]
norm_stats.reset_index(inplace = True, drop = True) 
# -

norm_stats.head()

norm_stats['date'] = pd.to_datetime((norm_stats.year*10000+norm_stats.month*100+norm_stats.day).apply(str),format='%Y%m%d')

date_ranges = philly_patterns[['date_range_start', 'date_range_end']].drop_duplicates().sort_values('date_range_start').reset_index(drop = True)

date_ranges


def get_date_range(x): 
    mask = (date_ranges['date_range_start'] <= x) & (date_ranges['date_range_end'] > x)
    return date_ranges.loc[mask,'date_range_start'].values[0]


norm_stats['date_range_start'] = norm_stats.date.apply(get_date_range)

norm_stats = norm_stats[['date_range_start', 'total_visits', 'total_devices_seen']].groupby('date_range_start').agg('sum').reset_index()

philly_patterns = philly_patterns.merge(norm_stats, on = 'date_range_start')
philly_patterns['visit_count_norm'] = 1000000*philly_patterns['raw_visit_counts']/philly_patterns['total_visits']
philly_patterns['visitor_count_norm'] = 1000000*philly_patterns['raw_visitor_counts']/philly_patterns['total_devices_seen']

check = philly_patterns.groupby(['safegraph_place_id']).nunique()

check.loc[check['B19013e1']>1]

philly_patterns[philly_patterns['safegraph_place_id'] == 'sg:002e331cebec461d932e13cee0e01ca3']

mask = philly_patterns['date_range_start'] < '2020-03-01' 
prepandemic = philly_patterns.loc[mask,].copy()
pandemic = philly_patterns.loc[~mask,].copy()

prepandemic.date_range_start.unique()

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler


# +
def cut_outliers(df, cols):
    for col in cols:
        ub = 2*df[col].quantile(0.99)
        mask = df[col] > ub
        df.loc[mask,col] = ub
    return df

def make_dummies(df, cols):
    return pd.get_dummies(df)
    

def Stdrshp(df, pivot_cols, constant_cols): 
    constant = df.loc[
        df.groupby('safegraph_place_id').cumcount() == 0,
        ['safegraph_place_id']+constant_cols
    ]
    constant.set_index('safegraph_place_id', inplace = True)
    #constant = make_dummies(constant, pivot_cols)
    pivoted = df.pivot(index='safegraph_place_id',columns='date_range_start')[pivot_cols]
    new_cols = [('{1} {0}'.format(*tup)) for tup in pivoted.columns]
    pivoted.columns = new_cols
    df = pivoted.join(constant)
    df = df.fillna(0)
    scaled = StandardScaler().fit_transform(df)
    scaled_features_df = pd.DataFrame(scaled, index=df.index, columns=df.columns)
    return scaled_features_df

def make_clusters(df, k = 5):
    kmeans = KMeans(n_clusters=k).fit(df)
    df['cluster'] = kmeans.labels_
    return (df, kmeans.inertia_)

def run_Kmeans(df, pivot_cols, constant_cols, k = 5):
    normalized = Stdrshp(cut_outliers(df, pivot_cols), pivot_cols, constant_cols)
    return make_clusters(normalized, k)


# -

pivot_cols = ['visit_count_norm', 'visitor_count_norm']
constant_cols = ['B19013e1','income_missing']

scaled = Stdrshp(cut_outliers(prepandemic, pivot_cols), pivot_cols, constant_cols)
n_clusters = [k for k in range(3,13)]
inertias = []
for k in n_clusters:
    df, inertia = make_clusters(scaled, k = k)
    inertias.append(inertia)

scaled.head()

import matplotlib.pyplot as plt

plt.plot(n_clusters,inertias)

# +
scaled = Stdrshp(prepandemic, pivot_cols, constant_cols)
n_clusters = [k for k in range(3,13)]
inertias = []
for k in n_clusters:
    df, inertia = make_clusters(scaled, k = k)
    inertias.append(inertia)
    
plt.plot(n_clusters,inertias)
# -

prepandemic, prepandemic_inertia = run_Kmeans(prepandemic, pivot_cols, constant_cols, k = 5)
pandemic, pandemic_inertia = run_Kmeans(pandemic, pivot_cols, constant_cols, k = 5)

pandemic_inertia

# We want to compare pandemic_inertia to the inertia we would have in 2020 if we used the
# 2019 clustering
combined = pandemic.join(prepandemic[['cluster']], how='left',rsuffix='pre')
combined[['clusterpre']] = combined[['clusterpre']].fillna(value=-1)
combined['clusterpre'] = combined['clusterpre'].astype(int)

cols = combined.columns
cols = cols[:-2]
n = len(combined.columns)


# +
def get_inertia(df, cols, cluster_col):
    mask = df[cluster_col] > -1
    n = len(df.index)
    filtered = df[mask]
    m = len(filtered.index)
    se = ((df[cols] - df.groupby(cluster_col)[cols].transform('mean'))**2).to_numpy().sum()
    return se*(n/m)

get_inertia(combined, cols, 'clusterpre')
# -

#worst case scenario
combined['newcol'] = 1
get_inertia(combined, cols, 'newcol')

combined.groupby('cluster').size()

combined.groupby('clusterpre').size()

big_cats = philly_patterns.groupby('top_category').size().sort_values(ascending = False)[:5]

big_cats = big_cats.index


def category_groups(x):
    if x in big_cats: 
        return x 
    else: 
        return 'other'


philly_patterns['top_category_simple'] = philly_patterns.top_category.apply(category_groups)

philly_patterns['top_category_simple'].unique()

mask = philly_patterns['date_range_start'] < '2020-03-01' 
cats = philly_patterns.loc[~mask,].copy()
cats = cats.loc[
        cats.groupby('safegraph_place_id').cumcount() == 0,
        ['safegraph_place_id', 'top_category_simple']
    ]

combined.drop(columns = 'top_category_simple', inplace = True)

combined = combined.join(cats.set_index('safegraph_place_id'), how = 'left')

pd.crosstab(combined.clusterpre, combined.top_category_simple, dropna=False)



StarcoderdataJupyterScriptsDedupFiltered