<jupyter_start><jupyter_text>Real World Facility Location Problem*Authors:* [<NAME>](https://github.com/gegen07), [<NAME>](https://github.com/jGaboardi), [<NAME>](https://github.com/ljwolf), [<NAME>](https://github.com/qszhao)This tutorial aims to show a facility location application. To achieve this goal the tuorial will solve a real world facility location problem using a dataset describing an area of census tract 205, San Francisco. The problem can be written as: store sites should supply the demand in this census tract considering the distance between the two sites, a demand point and a supply candidate site. This tutorial model this problem using four models already developed in literature: LSCP (Location Set Covering Problem), MCLP (Maximal Coverage Location Problem), P-Median and P-Center. If you don't know any of these models, it's recommended to see the notebooks that explain more deeply about each model: [LSCP](./lscp.ipynb), [MCLP](./mclp.ipynb), [P-Median](./p-median.ipynb), [P-Center](./p-center.ipynb)Besides that, the tutorial show how to use different solvers that [PULP](https://github.com/coin-or/pulp) support.<jupyter_code>import numpy
import geopandas
import pandas
import pulp
from shapely.geometry import Point
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>The problem is composed by 4 datafiles- network_distance is the calculus of distance between facility candidate sites calculated by ArcGIS Network Analyst Extension- demand_points represents the demand points with some important features for the facility location problem like population- facility_points represents the stores that are candidate facility sites- tract is the polygon of census tract 205.All datasets are online on this [repository](https://github.com/huanfachen/Open_source_location_cover_models/tree/master/Data/San_Francisco_store)<jupyter_code>DIRPATH = "../spopt/tests/data/"

network_distance = pandas.read_csv(DIRPATH+"SF_network_distance_candidateStore_16_censusTract_205_new.csv")
demand_points = pandas.read_csv(DIRPATH+ "SF_demand_205_centroid_uniform_weight.csv", index_col=0)
facility_points = pandas.read_csv(DIRPATH + "SF_store_site_16_longlat.csv", index_col=0)
tract = geopandas.read_file(DIRPATH + "ServiceAreas_4.shp")<jupyter_output><empty_output><jupyter_text>Display network_distance dataframe<jupyter_code>display(network_distance)<jupyter_output><empty_output><jupyter_text>Display demand points dataframe<jupyter_code>display(demand_points)<jupyter_output><empty_output><jupyter_text>Display facility_points dataframe<jupyter_code>display(facility_points)<jupyter_output><empty_output><jupyter_text>Plot tract<jupyter_code>tract.plot()<jupyter_output><empty_output><jupyter_text>To start modeling the problem assuming the arguments expected by `spopt.locate`, we should pass a numpy 2D array as a cost_matrix. So, first we pivot the network_distance dataframe. _Note that the columns and rows are in alphabetically order._<jupyter_code>ntw_dist_piv = network_distance.pivot_table(values="distance", index="DestinationName", columns="name")
ntw_dist_piv<jupyter_output><empty_output><jupyter_text>Here the pivot table is transformed to numpy 2D array such as `spopt.locate` expected. The matrix has a shape of 205x16.<jupyter_code>cost_matrix = ntw_dist_piv.to_numpy()
cost_matrix<jupyter_output><empty_output><jupyter_text>Now, as the rows and columns of our cost_matrix are sorted, so we have to sort our facility points and demand points geodataframes too.<jupyter_code>facility_points_gdf = geopandas.GeoDataFrame(
    facility_points,
    geometry=geopandas.points_from_xy(
        facility_points.long, facility_points.lat
    ),
).sort_values(by=['NAME']).reset_index()

demand_points_gdf = geopandas.GeoDataFrame(
    demand_points,
    geometry=geopandas.points_from_xy(demand_points.long, demand_points.lat),
).sort_values(by=['NAME']).reset_index()<jupyter_output><empty_output><jupyter_text>Here the rest of parameters are set. The ai is the demand weight and in this case we model as population in year 2000 of this tract.<jupyter_code>ai = demand_points_gdf['POP2000'].to_numpy()

service_dist = 5000
p_facility = 4<jupyter_output><empty_output><jupyter_text>Below, the method is used to plot the results of the four models that we will prepare to solve the problem.<jupyter_code>from matplotlib.patches import Patch
import matplotlib.lines as mlines

dv_colors = [
    "darkcyan",
    "mediumseagreen",
    "cyan",
    "darkslategray",
    "lightskyblue",
    "limegreen",
    "darkgoldenrod",
    "peachpuff",
    "coral",
    "mediumvioletred",
    "blueviolet",
    "fuchsia",
    "thistle",
    "lavender",
    "saddlebrown",
] 

def plot_results(model, facility_points_gdf, demand_points_gdf, facility_count, title):
    arr_points = []
    fac_sites = []
    
    for i in range(facility_count):
        if model.fac2cli[i]:

            geom = demand_points_gdf.iloc[model.fac2cli[i]]['geometry']
            arr_points.append(geom)
            fac_sites.append(i)

    fig, ax = plt.subplots(figsize=(12, 12))
    legend_elements = []

    tract.plot(ax=ax, alpha=1, color='tan', zorder=1)
    legend_elements.append(mlines.Line2D(
        [],
        [],
        color='tan',
        label='tract',
    ))

    facility_points_gdf.plot(ax=ax, color='brown', marker="*", markersize=80, zorder=2)
    legend_elements.append(mlines.Line2D(
        [],
        [],
        color='brown',
        marker="*",
        linewidth=0,
        label=f'facility sites ($n$={facility_count})'
    ))

    demand_points_gdf.plot(ax=ax, edgecolor='k', color='white', zorder=2)
    legend_elements.append(mlines.Line2D(
        [],
        [],
        color='white',
        marker="o",
        markeredgecolor='k',
        linewidth=0,
        label=f'clients not covered'
    ))
 
    for i in range(len(arr_points)):
        gdf = geopandas.GeoDataFrame(arr_points[i])

        label = f"points covered by {facility_points_gdf.iloc[[fac_sites[i]]]['NAME'].squeeze()}"
        legend_elements.append(Patch(facecolor=dv_colors[i], edgecolor="k", label=label))

        gdf.plot(ax=ax, zorder=3, alpha=0.7, edgecolor="k", color=dv_colors[i], label=label)
        facility_points_gdf.iloc[[fac_sites[i]]].plot(ax=ax,
                                marker="*",
                                markersize=200 * 3.0,
                                alpha=0.8,
                                zorder=4,
                                edgecolor="k",
                                facecolor=dv_colors[i])
        
        legend_elements.append(mlines.Line2D(
            [],
            [],
            color=dv_colors[i],
            marker="*",
            ms=20 / 2,
            markeredgecolor="k",
            linewidth=0,
            alpha=0.8,
            label=f"{facility_points_gdf.iloc[[fac_sites[i]]]['NAME'].squeeze()}",
        ))

    plt.title(title, fontweight="bold")
    plt.legend(handles = legend_elements, loc='upper left', bbox_to_anchor=(1.05, 1))<jupyter_output><empty_output><jupyter_text>LSCP<jupyter_code>from spopt.locate.coverage import LSCP

lscp = LSCP.from_cost_matrix(cost_matrix, service_dist)

lscp = lscp.solve(pulp.GLPK(msg=False))

lscp.facility_client_array()
plot_results(lscp, facility_points_gdf, demand_points_gdf, facility_points_gdf.shape[0], "LSCP")
lscp.problem.objective.value()<jupyter_output><empty_output><jupyter_text>MCLP<jupyter_code>from spopt.locate.coverage import MCLP

mclp = MCLP.from_cost_matrix(
    cost_matrix,
    ai,
    max_coverage=service_dist,
    p_facilities=4,
)

mclp = mclp.solve(pulp.GLPK(msg=False))

mclp.facility_client_array()
plot_results(mclp, facility_points_gdf, demand_points_gdf, facility_points_gdf.shape[0], "MCLP")
mclp.problem.objective.value()<jupyter_output><empty_output><jupyter_text>MCLP Percentage Covered<jupyter_code>mclp.uncovered_clients()
mclp.get_percentage()
mclp.percentage<jupyter_output><empty_output><jupyter_text>P-Center<jupyter_code>from spopt.locate import PCenter
pcenter = PCenter.from_cost_matrix(
    cost_matrix, p_facilities=p_facility
)

pcenter = pcenter.solve(pulp.GLPK(msg=False))
pcenter.facility_client_array()
plot_results(pcenter, facility_points_gdf, demand_points_gdf, facility_points_gdf.shape[0], "P-Center")
pcenter.problem.objective.value()<jupyter_output><empty_output><jupyter_text>P-Median<jupyter_code>from spopt.locate import PMedian
pmedian = PMedian.from_cost_matrix(
    cost_matrix, ai, p_facilities=p_facility
)

pmedian = pmedian.solve(pulp.GLPK(msg=False))
pmedian.facility_client_array()
plot_results(pmedian, facility_points_gdf, demand_points_gdf, facility_points_gdf.shape[0], "P-Median")
pmedian.problem.objective.value()
pmedian.get_mean_distance(weight=ai)
pmedian.mean_dist<jupyter_output><empty_output><jupyter_text>Using different solver for one modelFor this task we get the MCLP model as example.<jupyter_code>from spopt.locate.coverage import MCLP

mclp = MCLP.from_cost_matrix(
            cost_matrix,
            ai,
            max_coverage=service_dist,
            p_facilities=4,
        )<jupyter_output><empty_output><jupyter_text>First, it's good to see which solver you have installed in your machine.<jupyter_code>solvers = pulp.listSolvers(onlyAvailable=True)
solvers<jupyter_output><empty_output><jupyter_text>Above we can see that it returns a list with different solvers that are available. So, it's up to user choose the best solver that fits the model. Let's get the percentage as a metric to evaluate which solver is the best or improves the model.<jupyter_code>import time

for i in range(len(solvers)):
    solver = pulp.get_solver(solvers[i], msg=False)

    start = time.time()
    mclp = mclp.solve(solver)
    end = time.time()

    mclp.facility_client_array()
    mclp.uncovered_clients()
    mclp.get_percentage()

    percentage = mclp.percentage
    time_spent = end-start

    print(f"Solver: {solvers[i]} - Percentage: {percentage} - Time Spent: {time_spent}")<jupyter_output>Solver: GLPK_CMD - Percentage: 0.8975609756097561 - Time Spent: 0.06669068336486816
Solver: CPLEX_CMD - Percentage: 0.8975609756097561 - Time Spent: 0.08563423156738281
Solver: GUROBI_CMD - Percentage: 0.8975609756097561 - Time Spent: 0.07414436340332031
Solver: PULP_CBC_CMD - Percentage: 0.8975609756097561 - Time Spent: 0.06678557395935059
Solver: COIN_CMD - Percentage: 0.8975609756097561 - Time Spent: 0.05470633506774902
{'fraction_non_alphanumeric': 0.10621414913957936, 'fraction_numerical': 0.025047801147227534, 'mean_word_length': 4.573255194459244, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 80, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '425563', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3766, 'n_tokens_neox': 3515, 'n_words': 722, 'path': 'docs/notebooks/facloc-real-world.ipynb', 'repo_name': 'gis-ops/spopt', 'size': '10490'}