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

# # Analyzing connections in Baltimore shooting data
#
# This notebook takes the data standardized in the first step and runs various calculations on it — primarily related to cases connected by a common suspect or victim.

import pandas as pd
import networkx as nx
import geopandas as gpd
from shapely.geometry import Point, LineString
from pyproj import Proj, transform
import matplotlib.pyplot as plt


# This function lets us write basic findings to a logfile:

def log(obj, description):
    print(description + ":")
    with open("../outputs/log.txt", "a") as f:
        f.write(f"{description}:\n\n{obj}\n\n---\n")
    return obj


# ## Load persons named in Baltimore database

# +
persons = pd.read_csv("../inputs/baltimore-shootings-persons.csv")

persons.head().T


# -

# # Calculate rolling case-counts and closure rates

def get_annual_rolling_rates(df):
    return (
        df
        [[
            "CASE_DATE",
            "CASE_NUMBER",
            "CASE_STATUS",
        ]]
        .drop_duplicates()
        .assign(
            is_closed = lambda df: df["CASE_STATUS"] == "CLOSED",
            CASE_DATE = lambda df: pd.to_datetime(df["CASE_DATE"])
        )
        .groupby([
            pd.Grouper(
                key = "CASE_DATE",
                freq = "D",
            ),
        ])
        .pipe(lambda grp: pd.DataFrame({
            "incidents_total": grp.size(),
            "incidents_closed": grp["is_closed"].sum(),
        }))
        .fillna(0)
        .astype(int)
        .rolling(365)
        .sum()
        .dropna()
        .assign(
            closure_rate =  lambda df: df["incidents_closed"] / df["incidents_total"]
        )
        .reset_index()
    )


# +
rolling_annual_case_counts = (        
    persons
    .loc[lambda df: df["VICT_TYPE"] == "NON FATAL SHOOTING"]
    .pipe(get_annual_rolling_rates)
    [[ "CASE_DATE", "closure_rate" ]]
    .merge(
        persons
        .loc[lambda df: df["VICT_TYPE"] == "HOMICIDE SHOOTING"]
        .pipe(get_annual_rolling_rates)
        [[ "CASE_DATE", "closure_rate" ]],
        on = [ "CASE_DATE" ],
        how = "left",
        suffixes = [ "_nonfatal", "_fatal" ]
    )
    .merge(
        persons
        .pipe(get_annual_rolling_rates)
        [[ "CASE_DATE", "incidents_total" ]],
        how = "left",        
        on = [ "CASE_DATE" ],
        suffixes = [ "_nonfatal", "_fatal" ]
    )
    .set_index("CASE_DATE")
    .loc["2013":"2016"]
    .astype({ "incidents_total": int })
)

rolling_annual_case_counts.to_csv(
    "../outputs/baltimore-rolling-case-counts-and-closure-rate.csv"
)

rolling_annual_case_counts.head()

# +
fig, ax = plt.subplots(
    figsize = (10, 6),
)

(
    rolling_annual_case_counts
    [[
        "closure_rate_nonfatal",
        "closure_rate_fatal"
    ]]
    .plot(
        ax = ax,
        lw = 2,
        alpha = 0.75
    )
)

ax.figure.set_facecolor("white")

ax.yaxis.grid()
ax.set_yticks(pd.np.arange(0, 1.1, 0.1))
ax.set_yticklabels([ "{:.0f}%".format(y * 100)
    for y in ax.get_yticks() ])

ax.set_xlabel(None)
ax.set_ylabel("Closure Rate")
ax.set_title("Closure Rate — 12-Month Average")

ax.set_ylim(0, 1)
ax.axvline("2015-04-19", color = "red")
pass

# +
fig, ax = plt.subplots(
    figsize = (10, 3),
)

(
    rolling_annual_case_counts
    ["incidents_total"]
    .plot(
        lw = 3,
        color = "gray"
    )
)

ax.axvline("2015-04-19", color = "red")

ax.set_ylim(0, None)
ax.yaxis.grid(True)

ax.figure.set_facecolor("white")

ax.set_xlabel(None)
ax.set_xticks(ax.get_xticks())

ax.legend(["Total Incidents"])
pass
# -

log((    
    rolling_annual_case_counts
    .loc["2015-04-19"]
    ["incidents_total"]
), "Number of shooting cases in 12-month period ending 2015-04-19")

log((    
    rolling_annual_case_counts
    .loc["2016-04-19"]
    ["incidents_total"]    
), "Number of shooting cases in 12-month period ending 2016-04-19")

persons["CASE_STATUS"].value_counts()

# Ensure that there's only one case status per incident
assert (
    persons
    .groupby([
        "CASE_NUMBER",
    ])
    ["CASE_STATUS"]
    .nunique()
    .max()
) == 1

log(round((
    persons
    .loc[lambda df: df["CASE_DATE"] >= "2015-04-20"]
    .loc[lambda df: df["CASE_DATE"] < "2016"]
    [[
        "CASE_NUMBER",
        "CASE_STATUS"
    ]]
    .drop_duplicates()
    .pipe(lambda df: df["CASE_STATUS"] == "CLOSED")
    .mean()
), 3), "Closure rate for cases from April 20, 2015 through end of 2015")

log(round(
    persons
    .loc[lambda df: df["CASE_DATE"] >= "2015"]
    .loc[lambda df: df["CASE_DATE"] < "2017"]
    [[
        "CASE_NUMBER",
        "CASE_STATUS"
    ]]
    .drop_duplicates()
    .pipe(lambda df: df["CASE_STATUS"] == "OPEN")
    .sum()
), "Cases from 2015 and 2016 marked 'OPEN', as of mid-2017")

# ## Calculate basic statistics on cases, suspects, and victims

log((    
    persons["CASE_NUMBER"].nunique()
), "Number of distinct case numbers")

log((    
    persons["CASE_DATE"].min()
), "Earliest case date")

log((    
    persons["CASE_DATE"].max()
), "Latest case date")

log(round((
    persons
    .loc[lambda df: df["DESCRIPTION"] == "SUSPECT"]
    .loc[lambda df: df["name_is_vague"] == 0]
    ["CASE_NUMBER"]
    .nunique() 
    / 
    persons["CASE_NUMBER"].nunique()
), 3), "Proportion of cases with a named suspect")

(
    persons
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    ["VICT_TYPE"]
    .value_counts(normalize = True)
    .round(3)
    .pipe(log, "Proportion of victims by victim type")
)

(
    persons
    .loc[lambda df: df["DESCRIPTION"] == "SUSPECT"]
    ["case_vict_type_generic"]
    .value_counts(normalize = True)
    .round(3)
    .pipe(log, "Proportion of suspects by victim type")
)

(
    persons
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    ["RACE"]
    .value_counts(normalize = True)
    .round(3)
    .pipe(log, "Proportion of victims by race")
)

# ## Connect cases via shared victims/suspects

# +
# Here, we remove "persons of interest", and also deduplicate the handful
# of suspects who are listed more than once for the same case

person_victims_and_suspects = (   
    persons
    .loc[lambda df: df["DESCRIPTION"] != "PERSON OF INTEREST"]
    [[
        "CASE_NUMBER",
        "CASE_DATE",
        "DESCRIPTION",
        "VICT_TYPE",
        "person_uid"
    ]]
    .drop_duplicates()
)

person_victims_and_suspects.head()
# -

# By joining the persons data to itself, we can identify instances where two cases share a common person — someone who is either a suspect or victim in both cases (although not necessarily the same role in both):

# +
case_connections = (
    person_victims_and_suspects
    .drop(columns = [ "VICT_TYPE", "CASE_DATE" ])
    .pipe(lambda df: (
        df
        .merge(
            df,
            how = "left",
            on = [ "person_uid" ]
        )
    ))
    # Here, we limit the results to those where the matching case
    # is higher in the sort-order than the starting case,
    # to remove duplicates and self-matches
    .loc[lambda df: df["CASE_NUMBER_x"] < df["CASE_NUMBER_y"]]
)

case_connections.head()
# -

log((    
    case_connections
    ["person_uid"]
    .nunique()
), "Number of people connected to multiple cases")

(
    pd.concat([
        case_connections["CASE_NUMBER_x"],
        case_connections["CASE_NUMBER_y"],
    ])
    .to_frame("CASE_NUMBER")
    .nunique()
    .pipe(log, "Number of shootings (unique_case_numbers) connected by a common victim/suspect")
)

# +
connected_cases = (
    case_connections
    [[
        "CASE_NUMBER_x",
        "CASE_NUMBER_y"
    ]]
    .drop_duplicates()
)

connected_cases.head()
# -

# Here we create a network graph (via `networkx`) to identify groups of cases connected to one another — directly or indirectly:

# +
graph = nx.Graph()

graph.add_edges_from(
    connected_cases
    [[
        "CASE_NUMBER_x",
        "CASE_NUMBER_y"
    ]]
    .apply(tuple, axis = 1)
    .values
)

# +
groups = pd.concat([
    pd.DataFrame([ { "CASE_NUMBER": x, "group": i }  for x in subset ])
    for i, subset in enumerate(nx.connected_components(graph)) 
])

groups.head()
# -

log((
    groups["group"].nunique()
), "Number of groups of connected shootings")

# The table below counts the proportion of grouped cases of each size — 76% of groups are just two cases; 14% are groups of three cases, 5% are groups of four cases, et cetera:

(
    groups
    ["group"]
    .value_counts()
    .value_counts(normalize = True)
    .sort_index()
    .to_frame("proportion")
)

# ## Calculate statistics about these groups

# Most of the groups involve both a fatal and non-fatal shooting:

(
    person_victims_and_suspects
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    .merge(
        groups,
        on = "CASE_NUMBER",
        how = "left",
        validate = "m:1"
    )
    .dropna(subset = [ "group" ])
    .groupby([
        "group"
    ])
    ["VICT_TYPE"]
    .unique()
    .apply(sorted)
    .apply(tuple)
    .value_counts()
    .pipe(log, "Connected-shooting groups, by distinct victim types")
)

# If you count only the highest-level offense for each case — i.e., only homicide for cases that involved both fatal and non-fatal shootings — the numbers are similar, but slightly different:

# +
highest_level_victim_by_case = (
    persons
    .dropna(subset = [ "VICT_TYPE" ])
    .groupby("CASE_NUMBER")
    ["VICT_TYPE"]
    .min() # Here, `min` selects "HOMICIDE SHOOTING" over "NON FATAL SHOOTING"
    .reset_index()
)

highest_level_victim_by_case.head()
# -

(
    highest_level_victim_by_case
    .merge(
        person_victims_and_suspects
        .loc[lambda df: df["DESCRIPTION"] == "VICTIM"],
        how = "left",
        on = "CASE_NUMBER",
        suffixes = [ "_highest", "_individual" ]
    )
    .groupby([
        "VICT_TYPE_highest",
        "VICT_TYPE_individual"
    ])
    .size()
    .to_frame("count")
    .reset_index()
)

(
    groups
    .merge(
        highest_level_victim_by_case,
        how = "left",
        on = "CASE_NUMBER",
        validate = "1:1",
    )
    .groupby("group")
    ["VICT_TYPE"]
    .unique()
    .apply(sorted)
    .apply(tuple)
    .value_counts()
    .pipe(log, "Connected-shooting groups, by distinct *highest-level* victim types")    
)

# ## Calculate statistics about victims and suspects named in multiple cases

# Here, we count the number of distinct cases in which someone was listed as a victim, and the same for as a suspect:

# +
connected_person_counts = (
    person_victims_and_suspects
    
    # Select only people who serve as a connection between cases
    .loc[lambda df: df["person_uid"].isin(
        case_connections["person_uid"]
    )]
    [[
        "person_uid",
        "CASE_NUMBER",
        "DESCRIPTION"
    ]]
    .groupby([ "person_uid" ])
    ["DESCRIPTION"]
    .value_counts()
    .unstack()
    .fillna(0)
    .astype(int)
)

connected_person_counts.head()
# -

log((
    connected_person_counts
    .loc[lambda df: df["SUSPECT"] > 0]
    .loc[lambda df: df["VICTIM"] > 0]
    .pipe(len)
), "People listed as both victims and suspects")

# ### Victims

log((
    connected_person_counts
    .loc[lambda df: df["VICTIM"] > 1]
    .pipe(len)
), "People listed as victims in multiple cases")

# For people listed as victims in multiple cases, for how many was their final recorded case a homicide vs. non-fatal shooting?:

(
    person_victims_and_suspects
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    .loc[lambda df: df["person_uid"].isin(
        connected_person_counts
        .loc[lambda df: df["VICTIM"] > 1]
        .index
    )]
    .sort_values("CASE_DATE")
    .groupby("person_uid")
    ["VICT_TYPE"]
    .last()
    .value_counts()
    .pipe(log, "Final victim type for persons shot in multiple incidents")
)

# As a sanity-check, the *initial* case types for these victims should all be non-fatal — and they are:

# Same code as above, but replacing `.last()` with `.first()`:
(
    person_victims_and_suspects
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    .loc[lambda df: df["person_uid"].isin(
        connected_person_counts
        .loc[lambda df: df["VICTIM"] > 1]
        .index
    )]
    .sort_values("CASE_DATE")
    .groupby("person_uid")
    ["VICT_TYPE"]
    .first()
    .value_counts()
    .pipe(log, "Initial victim type for persons shot in multiple incidents")
)

# ### Suspects

log((
    connected_person_counts
    .loc[lambda df: df["SUSPECT"] > 1]
    .pipe(len)
), "People listed as suspects in multiple cases")

# Here, we identify cases where at least one of the suspects was also named a suspect in another case:

# +
cases_with_repeat_suspects = (
    person_victims_and_suspects
    .loc[lambda df: (
        df["person_uid"].isin(
            connected_person_counts
            .loc[lambda df: df["SUSPECT"] > 1]
            .index
        ) &
        (df["DESCRIPTION"] == "SUSPECT")
    )]    
    ["CASE_NUMBER"]
    .drop_duplicates()
)

cases_with_repeat_suspects.head()
# -

len(cases_with_repeat_suspects)

(
    person_victims_and_suspects
    .loc[lambda df: df["CASE_NUMBER"].isin(cases_with_repeat_suspects)]
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    ["VICT_TYPE"]
    .value_counts()
    .pipe(log, "Victim counts for cases with multiple-case suspects")
)

# ## Demographics of people named in multiple cases

# There's one person in the data with conflicting information about their sex:

(
    persons
    .loc[lambda df: df["DESCRIPTION"] != "PERSON OF INTEREST"]
    .replace("UNKNOWN", pd.np.nan)        
    .groupby([
        "person_uid"
    ])
    [[ "RACE", "SEX" ]]
    .nunique()
    .sort_values([ "RACE", "SEX" ])
    .loc[lambda df: (df > 1).sum(axis = 1) > 0]
)    

# To handle this exception, and to use the person's age at the time of their *last recorded* case, we use the most recent case on file for each person:

# +
cx_latest_demographics = (
    persons
    # Select only persons who serve as connections between cases
    .loc[lambda df: df["person_uid"].isin(case_connections["person_uid"])]
    .sort_values("CASE_DATE", ascending = False)
    .drop_duplicates(subset = [ "person_uid" ])
    .replace("UNKNOWN", pd.np.nan)
    .groupby([
        "person_uid"
    ])
    [[ "RACE", "SEX", "age_at_case", "CASE_DATE" ]]
    .last()
    .reset_index()
)

cx_latest_demographics.head()
# -

len(cx_latest_demographics)

cx_latest_demographics["age_at_case"].value_counts().sort_index()

cx_latest_demographics.isnull().sum()

assert cx_latest_demographics.isnull().sum().sum() == 0


def classify_age(age):
    if age < 18:
        return "group0: <18"
    if age < 28:
        return "group1: 18–27"
    if age < 38:
        return "group2: 28-37"
    return "group3: 38+"


# +
cx_demographic_groups = (
    cx_latest_demographics
    .assign(
        age_group = lambda df: df["age_at_case"].apply(classify_age)
    )
    .groupby([
        "age_group",
        "SEX",
        "RACE"
    ])
    .size()
    .unstack()
    .fillna(0)
    .astype(int)
    .reset_index()
    .assign(
        age_group = lambda df: df["age_group"].str.replace(r".*:", "")
    )
    .set_index("age_group")
    .pipe(log, "Demographic groupings of mulitply-involved persons")
)

cx_demographic_groups
# -

# # Spatial analysis
#
# In this section, we identify connected cases that, based on the addresses provided, appear to cross the Baltimore Police Department's police-district boundaries.

# +
districts = (
    gpd
    .read_file("../inputs/Police_Districts/Police_Districts.shp")
    .to_crs(epsg = 4326)
)

districts
# -

districts.plot();

# ## Join districts to other spatial data

# +
geocoded_addresses = (
    pd.read_csv(
        "../inputs/baltimore-case-addresses.csv"
    )
    .assign(
        lng = lambda x: x["lng"].fillna(0.0),
        lat = lambda x: x["lat"].fillna(0.0),
    )
    .assign(
        geometry = lambda x: x.fillna(0.0).apply(
            lambda row: Point(row["lng"], row["lat"]),
            axis = 1
        )
    )
    .pipe(gpd.GeoDataFrame)
)

geocoded_addresses.crs = districts.crs

geocoded_addresses.head()

# +
incident_locations = (    
    geocoded_addresses
    .pipe(
        gpd.sjoin,
        districts
        [[
            "Dist_Name",
            "geometry"
        ]].rename(columns = {
            "Dist_Name": "police_district",
        }),
        how = "left",
    )
    .drop(columns = [
        "index_right"
    ])
    [[
        "agency_incident_id",
        "location_address",
        "latlng_method",
        "geocoded_location_interpreted",
        "geocoded_location_type",
        "geocoding_imprecise",
        "geometry",
        "police_district"
    ]]
    .drop_duplicates(subset = [
        "agency_incident_id",
        "geocoded_location_interpreted"
    ])
)

incident_locations.head()
# -

(
    person_victims_and_suspects
    .loc[lambda df: df["DESCRIPTION"] == "VICTIM"]
    .pipe(len)
)

assert incident_locations["agency_incident_id"].nunique() == len(incident_locations)

(
    incident_locations
    ["police_district"]
    .fillna("[blank]")
    .value_counts()
)

(
    geocoded_addresses
    ["geocoding_imprecise"]
    .fillna("[blank]")
    .value_counts()
)

# Note: Per above, there is one incident not geocoded to within a Baltimore police district. The address appears to be just outside district lines. (It is not a case connected to others in this analysis.) There are also two addresses that Google's API reported could only be imprecisely geocoded; all plausible locations for each of those addresses, however, are contained within a single police district.

# ## Analyze geography of connected cases

# +
# Identify the districts for each of the cases,
# for any two connected cases
connected_cases_districts = (
    connected_cases
    .merge(
        groups[[ "CASE_NUMBER", "group" ]],
        how = "left",
        left_on = [ "CASE_NUMBER_x" ],
        right_on = [ "CASE_NUMBER" ],
        validate = "m:1",        
    )
    .drop(columns = [ "CASE_NUMBER" ])
    .merge(
        incident_locations[[ "agency_incident_id", "geometry", "police_district" ]],
        how = "left",
        left_on = [ "CASE_NUMBER_x" ],
        right_on = [ "agency_incident_id" ],
        validate = "m:1",
    )
    .drop(columns = [ "agency_incident_id" ])
    .merge(
        incident_locations[[ "agency_incident_id", "geometry", "police_district" ]],
        how = "left",
        left_on = [ "CASE_NUMBER_y" ],
        right_on = [ "agency_incident_id" ],
        validate = "m:1",        
    )
    .drop(columns = [ "agency_incident_id" ])
)

connected_cases_districts.head()
# -

log(
    connected_cases_districts
    .pipe(lambda df: df["police_district_x"] != df["police_district_y"])
    .mean()
    .round(3),
    "Proportion of *case connections* that are multi-district"
)

(    connected_cases_districts
    .assign(
        same_district = lambda df: df["police_district_x"] == df["police_district_y"]
    )
    .melt(
        id_vars = [ "same_district" ],
        value_vars = [ "CASE_NUMBER_x", "CASE_NUMBER_y" ]
    )
     .groupby(["value"])
    ["same_district"]
    .mean()
    .pipe(lambda x: (x < 1).mean())
    .round(3),
)

log(
    connected_cases_districts
    .assign(
        same_district = lambda df: df["police_district_x"] == df["police_district_y"]
    )
    .melt(
        id_vars = [ "same_district" ],
        value_vars = [ "CASE_NUMBER_x", "CASE_NUMBER_y" ]
    )
    .groupby(["value"])
    ["same_district"]
    .mean()
    .pipe(lambda x: (x < 1).mean())
    .round(3),
    "Proportion of *connected cases* that have at least one out-of-district connection"
)

# The most common district combinations:

(
    connected_cases_districts
    .groupby([
        "police_district_x",
        "police_district_y",
    ])
    .size()
    .unstack()
    .fillna(0)
    .astype(int)
)

(
    connected_cases_districts
    .loc[lambda df: df["police_district_x"] != df["police_district_y"]]
    .assign(
        combo = lambda df: df[[ "police_district_x", "police_district_y" ]].apply(sorted, axis = 1).apply(tuple)
    )
    ["combo"]
    .value_counts()
    .loc[lambda x: x >= 10]
)


# # Mapping

def tee(df, fn, *args, **kwargs):
    fn(df, *args, **kwargs)
    return df


# +
multi_district_lines = gpd.GeoDataFrame(
    connected_cases_districts
    .loc[lambda df: df["police_district_x"] != df["police_district_y"]]
    .assign(
        geometry = lambda df: df[[ "geometry_x", "geometry_y" ]].apply(LineString, axis = 1),
        combo = lambda df: df[[ "police_district_x", "police_district_y" ]].apply(sorted, axis = 1).apply(tuple)        
    )
)

multi_district_lines.head()


# -

def plot_connected_shootings(df, ax):
    districts.plot(
        ax = ax,
        color = "white",
        edgecolor = "black",
    )
    
    (
        df
        .pipe(tee, lambda df: (
            gpd.GeoDataFrame({
                "geometry": df["geometry_x"]
            })
            .plot(
                ax = ax,
                marker = "o",
                markersize = 8,
                color = "black",
                alpha = 0.3
            )
        ))
        .pipe(tee, lambda df: (
            gpd.GeoDataFrame({
                "geometry": df["geometry_y"]
            })
            .plot(
                ax = ax,
                marker = "o",
                markersize = 8,
                color = "black",
                alpha = 0.3
            )
        ))
        .pipe(tee, lambda df: (
            df[["geometry"]]
            .plot(
                ax = ax,
                alpha = 0.3
            )
        ))    
    )
    
    ax.set_axis_off()


# +
fig, ax = plt.subplots(
    figsize = (10, 10)
)

plot_connected_shootings(multi_district_lines, ax)

ax.set_title("Multi-District Connected Shootings")
pass
# -

# ---
#
# ---
#
# ---

StarcoderdataJupyterScriptsDedupFiltered