<jupyter_start><jupyter_text>Number of unique safegraph ids:<jupyter_code>philly_patterns_all['safegraph_place_id'].nunique()<jupyter_output><empty_output><jupyter_text>Get POIs with category info:<jupyter_code>categories = (
    philly_patterns_all.dropna(subset=["top_category"])
    .drop_duplicates(subset=["safegraph_place_id", "top_category"])
    .sort_values("date_range_start", ascending=False)
    .drop_duplicates(subset=["safegraph_place_id"])
)[["safegraph_place_id", "top_category"]]
categories.head()<jupyter_output><empty_output><jupyter_text>Number of POI with category info:<jupyter_code>categories['safegraph_place_id'].nunique()<jupyter_output><empty_output><jupyter_text>Merge in category info and drop those without a category:<jupyter_code>philly_patterns = philly_patterns_all.drop(labels=["top_category"], axis=1).merge(
    categories, on="safegraph_place_id"
)<jupyter_output><empty_output><jupyter_text>No more missing categories:<jupyter_code>philly_patterns['top_category'].isnull().sum()<jupyter_output><empty_output><jupyter_text>Trim to the columns we need:<jupyter_code>philly_patterns = philly_patterns[
    [
        "safegraph_place_id",
        "date_range_start",
        "date_range_end",
        "raw_visit_counts",
        "raw_visitor_counts",
        "poi_cbg",
        "top_category",
        'visitor_home_cbgs',
        'visits_by_day'
    ]
]<jupyter_output><empty_output><jupyter_text>Make datetime columns:<jupyter_code>philly_patterns["date_range_start"] = pd.to_datetime(
    philly_patterns["date_range_start"].astype(str).apply(lambda x: x[:10])
)

philly_patterns["date_range_end"] = pd.to_datetime(
    philly_patterns["date_range_end"].astype(str).apply(lambda x: x[:10])
)
# philly_patterns["poi_cbg"] = philly_patterns["poi_cbg"].astype(int)
census_path = DATA_DIR / 'raw' / 'open-census-data'
files = [file for file in census_path.glob('**/cbg_b02.csv')]
census = pd.read_csv(files[0])
census =  census[['census_block_group','B02001e1', 'B02001e2', 'B02001e3']]
census['White_pct'] = census['B02001e2']/census['B02001e1']
census['Black_pct'] = census['B02001e3']/census['B02001e1']
census = census.set_index('census_block_group')
census = census.fillna(0)
census.loc[421010176013]
def str_to_dict(x):
    result = x.replace('{','')
    result = result.replace('}', '')
    result = result.replace('"', '')
    if not result:
        return dict()
    result = result.split(',')
    result = [x.split(':') for x in result]
    result = {int(x[0]):int(x[1]) for x in result}
    return result
def estimate_visitors(x, race_col):
    result = 0
    for cbg in x:
        try:
            result += census.loc[cbg,race_col]*x[cbg]
        except:
            continue
    return result

def total_cbg_visitors(x):
    return sum(x[cbg] for cbg in x)
    
def compute_race_columns(x):
    d = str_to_dict(x)
    return (
        estimate_visitors(d, 'Black_pct'), estimate_visitors(d, 'White_pct'), 
        total_cbg_visitors(d)
    )
#philly_patterns['black_visitors'] = philly_patterns.visitor_home_cbgs.apply(lambda x: estimate_visitors(str_to_dict(x), 'Black_pct'))
#philly_patterns['white_visitors'] = philly_patterns.visitor_home_cbgs.apply(lambda x: estimate_visitors(str_to_dict(x), 'White_pct'))
#philly_patterns['cbg_visitors'] = philly_patterns.visitor_home_cbgs.apply(lambda x: total_cbg_visitors(str_to_dict(x)))
race_data = pd.DataFrame(
    philly_patterns.visitor_home_cbgs.apply(compute_race_columns).tolist(), 
    philly_patterns.index
)
race_data.columns = ['black_pct', 'white_pct', 'cbg_visitors']
race_data = philly_patterns[['safegraph_place_id']].join(race_data)
race_data = race_data.groupby('safegraph_place_id').agg('sum')
race_data['black_pct'] = race_data['black_pct']/race_data['cbg_visitors']
race_data['white_pct'] = race_data['white_pct']/race_data['cbg_visitors']
race_data = race_data.fillna(0)
race_data = race_data[['black_pct', 'white_pct']]<jupyter_output><empty_output><jupyter_text>Histogram Features<jupyter_code>def str_to_list(x):
    result = x.replace('[','')
    result = result.replace(']', '')
    result = result.replace('"', '')
    if not result:
        return []
    result = result.split(',')
    result = [int(x) for x in result]
    return result
def make_bucket_counts(x, n_buckets):
    #output (
    # elements in x = 0, # elements in x = 1, # elements in x b/w 2 and 10, 
    # elements in x b/w 11 and 100, # elements in x bigger than 100
    #)
    d = {i:0 for i in range(n_buckets)}
    for n in x: 
        if n == 0: 
            d[0]+=1 
        else:
            i = min(int(np.ceil(np.log10(n))) + 1, n_buckets - 1)
            d[i] += 1 
    return d
col_names = ['0', '1', '2-10', '11-100', '>100']
col_names = [col+'_visits' for col in col_names]
histogram_features = pd.DataFrame(
    philly_patterns.visits_by_day.apply(
        lambda x: make_bucket_counts(str_to_list(x),5)
    ).tolist(), philly_patterns.index
)
histogram_features = philly_patterns[['safegraph_place_id']].join(histogram_features)
histogram_features = histogram_features.groupby('safegraph_place_id').agg('sum')
histogram_features.columns = col_names
histogram_features.head()
# 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-2020-12"


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)
    .dropna()
)
def get_date_range(x): 
    mask = (date_ranges['date_range_start'] <= x) & (date_ranges['date_range_end'] > x)
    if np.any(mask):
        return date_ranges.loc[mask,'date_range_start'].values[0]
    else:
        return np.nan
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"]
)
philly_patterns.head()<jupyter_output><empty_output><jupyter_text>Tests<jupyter_code>philly_patterns.top_category.value_counts().head(n=20)<jupyter_output><empty_output><jupyter_text>Let's check restaurants:<jupyter_code>restaurants = philly_patterns.query(
    "top_category == 'Restaurants and Other Eating Places'"
)
X = restaurants.groupby("date_range_start")["visitor_count_norm"].mean()

# Normalize by Jan 2019 (could also do Jan 2020)
(X / X.iloc[0]).plot()
X = restaurants.groupby("date_range_start")["visit_count_norm"].mean()

# Normalize by Jan 2019 (could also do Jan 2020)
(X / X.iloc[0]).plot()<jupyter_output><empty_output><jupyter_text>Let's check grocery stores:<jupyter_code>grocery_stores = philly_patterns.query("top_category == 'Grocery Stores'")
X = grocery_stores.groupby("date_range_start")["visitor_count_norm"].mean()

# Normalize by Jan 2019 (could also do Jan 2020)
(X / X.iloc[0]).plot()<jupyter_output><empty_output><jupyter_text>**Note: "visit_count_norm" seems weird!**<jupyter_code>X = grocery_stores.groupby("date_range_start")["visit_count_norm"].mean()

# Normalize by Jan 2019 (could also do Jan 2020)
(X / X.iloc[0]).plot()<jupyter_output><empty_output><jupyter_text>Analysis<jupyter_code>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 pivot(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)
    pivoted = df.pivot_table(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
    return pivoted

def Stdrshp(df, pivot_cols, auxiliaries):
    df = pivot(df, pivot_cols)
    for aux in auxiliaries:
         df = df.join(aux)
    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, n_init=50).fit(df) # NH NOTE: USE HIGHER N_INIT TO AVOID RANDOMNESS
    df['cluster'] = kmeans.labels_
    return (df, kmeans.inertia_)

def run_Kmeans(df, pivot_cols, auxiliaries, k = 5):
    normalized = Stdrshp(cut_outliers(df, pivot_cols), pivot_cols, auxiliaries)
    return make_clusters(normalized, k)
pivot_cols = ["visitor_count_norm"] ## NH NOTE: try using visitor_count_norm only
auxiliaries = [race_data, histogram_features]  #'B19013e1','income_missing']
import matplotlib.pyplot as plt
scaled = Stdrshp(philly_patterns, pivot_cols, auxiliaries)
n_clusters = [k for k in range(3,40)]
inertias = []
for k in n_clusters:
    print(k)
    df, inertia = make_clusters(scaled, k = k)
    inertias.append(inertia)
    
plt.plot(n_clusters,inertias)
scaled.head()
from kneed import KneeLocator
# Initialize the knee algorithm
kn = KneeLocator(n_clusters, inertias, curve='convex', direction='decreasing')
OPTIMAL_CLUSTERS = kn.knee
OPTIMAL_CLUSTERS
philly_patterns_clustered, pandemic_inertia = run_Kmeans(
    philly_patterns, pivot_cols, auxiliaries, k=OPTIMAL_CLUSTERS
)
# prepandemic, prepandemic_inertia = run_Kmeans(
#     prepandemic, pivot_cols, constant_cols, k=OPTIMAL_CLUSTERS
# )
# pandemic, pandemic_inertia = run_Kmeans(
#     philly_patterns, pivot_cols, constant_cols, k=OPTIMAL_CLUSTERS
# )
combined = philly_patterns_clustered.copy()
combined[['cluster']].to_csv(DATA_DIR / 'processed' / 'sg_clusters.csv', index = True)
# # 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, 'cluster')
#worst case scenario
combined['newcol'] = 1
get_inertia(combined, cols, 'newcol')
combined.groupby('cluster').size()
# combined.groupby('clusterpre').size()<jupyter_output><empty_output><jupyter_text>**NH**: Plot the mean value over time by cluster<jupyter_code>fig, ax = plt.subplots()
for cluster_num in combined['cluster'].unique():
    cluster = combined.query(f"cluster == {cluster_num}")
    cluster[cluster.columns[:-6]].T.mean(axis=1).plot(ax=ax, legend=False)<jupyter_output><empty_output><jupyter_text>**NH:** seems like it sort of worked in terms of differentiating both by overall foot traffic and the shape of the traffic over time<jupyter_code>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.cluster, combined.top_category_simple, dropna=False)
len(combined)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1413850255818981, 'fraction_numerical': 0.018087482885349858, 'mean_word_length': 5.662506000960153, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 66, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '633983', 'license': 'MIT', 'n_tokens_mistral': 5359, 'n_tokens_neox': 5054, 'n_words': 854, 'path': 'notebooks/NH-Places Clustering First Try.ipynb', 'repo_name': 'hkronenb/HK_DATS_Practicum', 'size': '13877'}