diff --git a/README.md b/README.md index eed0cacc..486c10e0 100644 --- a/README.md +++ b/README.md @@ -97,6 +97,7 @@ All examples can be run interactively by launching this repository as a [![Binde - [scikit-learn](https://scikit-learn.org/stable/) - [geopandas](https://geopandas.org/) - [pulp](https://coin-or.github.io/pulp/) +- [spaghetti](https://github.com/pysal/spaghetti) ## Installation spopt is available on the [Python Package Index](https://pypi.org/). Therefore, you can either install directly with pip from the command line: diff --git a/docs/api.rst b/docs/api.rst index f4e3dd4e..2aff5f73 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -24,4 +24,15 @@ Model based approaches for aggregating a large set of geographic units (with sma region.Spenc region.WardSpatial +Locate Methods +-------------- + +.. autosummary:: + :toctree: generated/ + + locate.coverage.LSCP + locate.coverage.MCLP + locate.PCenter + locate.PMedian + diff --git a/requirements.txt b/requirements.txt index 1d340818..8107d3cf 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,4 @@ pandas>=1 pulp scikit-learn>=0.22 scipy>=0.11 +spaghetti diff --git a/spopt/locate/base.py b/spopt/locate/base.py index 9c9bd5ab..9cfe0679 100644 --- a/spopt/locate/base.py +++ b/spopt/locate/base.py @@ -41,7 +41,7 @@ class BaseOutputMixin: def client_facility_array(self) -> None: """ - Create an array 2d $m$ x $n$, where m is number of clients and n is number of facilities. + Create an array 2d MxN, where m is number of clients and n is number of facilities. """ if hasattr(self, "fac2cli"): self.cli2fac = [[] for i in range(self.aij.shape[0])] @@ -63,7 +63,7 @@ def uncovered_clients(self) -> None: class CoveragePercentageMixin: """ - Mixin for calculate the percentage of area covered + Mixin to calculate the percentage of area covered """ def get_percentage(self): @@ -74,7 +74,22 @@ def get_percentage(self): class MeanDistanceMixin: + """ + Mixin to calculate the mean distance between demand and facility sites chosen + """ + def get_mean_distance(self, weight: np.array): + """ + Calculate the mean distance + Parameters + ---------- + weight: np.array + weight of all demand points + + Returns + ------- + None + """ self.mean_dist = self.problem.objective.value() / weight.sum() diff --git a/spopt/locate/coverage.py b/spopt/locate/coverage.py index f07e3c84..c93a0638 100644 --- a/spopt/locate/coverage.py +++ b/spopt/locate/coverage.py @@ -61,6 +61,51 @@ def from_cost_matrix( Returns ------- LSCP object + + Examples + -------- + >>> from spopt.locate.coverage import LSCP + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Calculate the cost matrix + >>> cost_matrix = ntw.allneighbordistances( + ... sourcepattern=ntw.pointpatterns["clients"], + ... destpattern=ntw.pointpatterns["facilities"]) + + Create LSCP instance from cost matrix + + >>> lscp_from_cost_matrix = LSCP.from_cost_matrix(cost_matrix, max_coverage=8) + >>> lscp_from_cost_matrix = lscp_from_cost_matrix.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> lscp_from_cost_matrix.facility_client_array() + >>> lscp_from_cost_matrix.fac2cli + """ r_fac = range(cost_matrix.shape[1]) @@ -116,6 +161,48 @@ def from_geodataframe( Returns ------- LSCP object + + Examples + -------- + >>> from spopt.locate.coverage import LSCP + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Create LSCP instance from cost matrix + + >>> lscp_from_geodataframe = LSCP.from_geodataframe(clients_snapped, facilities_snapped, + ... "geometry", "geometry", + ... max_coverage=8, distance_metric="euclidean") + >>> lscp_from_geodataframe = lscp_from_geodataframe.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> lscp_from_geodataframe.facility_client_array() + >>> lscp_from_geodataframe.fac2cli + """ dem = gdf_demand[demand_col] @@ -152,7 +239,7 @@ def from_geodataframe( def facility_client_array(self) -> None: """ - Create an array 2d $m$ x $n$, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the $facility_0$ cover the entire array. + Create an array 2d MxN, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the facility-i cover the entire array. Returns ------- @@ -250,6 +337,56 @@ def from_cost_matrix( Returns ------- MCLP object + + Examples + -------- + + >>> from spopt.locate.coverage import MCLP + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Calculate the cost matrix + + >>> cost_matrix = ntw.allneighbordistances( + ... sourcepattern=ntw.pointpatterns["clients"], + ... destpattern=ntw.pointpatterns["facilities"]) + + Simulate demand weights from 1 to 12 + + >>> ai = numpy.random.randint(1, 12, 100) + + Create MCLP instance from cost matrix + + >>> mclp_from_cost_matrix = MCLP.from_cost_matrix(cost_matrix, ai, max_coverage=7, p_facilities=4) + >>> mclp_from_cost_matrix = mclp_from_cost_matrix.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> mclp_from_cost_matrix.facility_client_array() + >>> mclp_from_cost_matrix.fac2cli """ r_fac = range(cost_matrix.shape[1]) r_cli = range(cost_matrix.shape[0]) @@ -313,6 +450,60 @@ def from_geodataframe( Returns ------- MCLP object + + Examples + -------- + >>> from spopt.locate.coverage import MCLP + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Simulate demand weights from 1 to 12 + + >>> ai = numpy.random.randint(1, 12, 100) + >>> clients_snapped['weights'] = ai + + Create MCLP instance from geodataframe + + >>> mclp_from_geodataframe = MCLP.from_geodataframe( + ... clients_snapped, + ... facilities_snapped, + ... "geometry", + ... "geometry", + ... "weights", + ... max_coverage=7, + ... p_facilities=4, + ... distance_metric="euclidean" + ... ) + + >>> mclp_from_geodataframe = mclp_from_geodataframe.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> mclp_from_geodataframe.facility_client_array() + >>> mclp_from_geodataframe.fac2cli """ service_load = gdf_demand[weights_cols].to_numpy() dem = gdf_demand[demand_col] @@ -351,7 +542,7 @@ def from_geodataframe( def facility_client_array(self) -> None: """ - Create an array 2d $m$ x $n$, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the $facility_0$ cover the entire array. + Create an array 2d MxN, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the facility-i cover the entire array. Returns ------- diff --git a/spopt/locate/p_center.py b/spopt/locate/p_center.py index f6587019..11ff3135 100644 --- a/spopt/locate/p_center.py +++ b/spopt/locate/p_center.py @@ -64,6 +64,53 @@ def from_cost_matrix( Returns ------- PCenter object + + Examples + -------- + + >>> from spopt.locate import PCenter + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Calculate the cost matrix + + >>> cost_matrix = ntw.allneighbordistances( + ... sourcepattern=ntw.pointpatterns["clients"], + ... destpattern=ntw.pointpatterns["facilities"]) + + Create PCenter instance from cost matrix + + >>> pcenter_from_cost_matrix = PCenter.from_cost_matrix(cost_matrix, p_facilities=4) + >>> pcenter_from_cost_matrix = pcenter_from_cost_matrix.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> pcenter_from_cost_matrix.facility_client_array() + >>> pcenter_from_cost_matrix.fac2cli + """ r_cli = range(cost_matrix.shape[0]) r_fac = range(cost_matrix.shape[1]) @@ -130,6 +177,52 @@ def from_geodataframe( Returns ------- PCenter object + + Examples + -------- + >>> from spopt.locate import PCenter + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Create PCenter instance from cost matrix + + >>> pcenter_from_geodataframe = PCenter.from_geodataframe( + ... clients_snapped, + ... facilities_snapped, + ... "geometry", + ... "geometry", + ... p_facilities=P_FACILITIES, + ... distance_metric="euclidean" + ... ) + >>> pcenter_from_geodataframe = pcenter_from_geodataframe.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> pcenter_from_geodataframe.facility_client_array() + >>> pcenter_from_geodataframe.fac2cli """ dem = gdf_demand[demand_col] @@ -166,7 +259,7 @@ def from_geodataframe( def facility_client_array(self) -> None: """ - Create an array 2d $m$ x $n$, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the $facility_0$ cover the entire array. + Create an array 2d MxN, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the facility-i cover the entire array. Returns ------- diff --git a/spopt/locate/p_median.py b/spopt/locate/p_median.py index 888df70d..5eb240e8 100644 --- a/spopt/locate/p_median.py +++ b/spopt/locate/p_median.py @@ -88,6 +88,56 @@ def from_cost_matrix( Returns ------- PMedian object + + Examples + -------- + + >>> from spopt.locate import PMedian + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Calculate the cost matrix + + >>> cost_matrix = ntw.allneighbordistances( + ... sourcepattern=ntw.pointpatterns["clients"], + ... destpattern=ntw.pointpatterns["facilities"]) + + Simulate demand weights from 1 to 12 + + >>> ai = numpy.random.randint(1, 12, 100) + + Create PMedian instance from cost matrix + + >>> pmedian_from_cost_matrix = PMedian.from_cost_matrix(cost_matrix, ai, p_facilities=4) + >>> pmedian_from_cost_matrix = pmedian_from_cost_matrix.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> pmedian_from_cost_matrix.facility_client_array() + >>> pmedian_from_cost_matrix.fac2cli """ r_cli = range(cost_matrix.shape[0]) r_fac = range(cost_matrix.shape[1]) @@ -156,6 +206,58 @@ def from_geodataframe( Returns ------- PMedian object + + Examples + -------- + + >>> from spopt.locate import PMedian + >>> from spopt.locate.util import simulated_geo_points + >>> import pulp + >>> import spaghetti + + Create regular lattice + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + Simulate points belong to lattice + + >>> demand_points = simulated_geo_points(street_buffered, needed=100, seed=5) + >>> facility_points = simulated_geo_points(street_buffered, needed=5, seed=6) + + Snap points to the network + + >>> ntw.snapobservations(demand_points, "clients", attribute=True) + >>> clients_snapped = spaghetti.element_as_gdf(ntw, pp_name="clients", snapped=True) + >>> ntw.snapobservations(facility_points, "facilities", attribute=True) + >>> facilities_snapped = spaghetti.element_as_gdf(ntw, pp_name="facilities", snapped=True) + + Simulate demand weights from 1 to 12 + + >>> ai = numpy.random.randint(1, 12, 100) + >>> clients_snapped['weights'] = ai + + Create PMedian instance from cost matrix + + >>> pmedian_from_geodataframe = PMedian.from_geodataframe( + ... clients_snapped, + ... facilities_snapped, + ... "geometry", + ... "geometry", + ... "weights", + ... p_facilities=P_FACILITIES, + ... distance_metric="euclidean") + >>> pmedian_from_geodataframe = pmedian_from_geodataframe.solve(pulp.PULP_CBC_CMD(msg=False)) + + Get facility lookup demand coverage array + + >>> pmedian_from_geodataframe.facility_client_array() + >>> pmedian_from_geodataframe.fac2cli """ service_load = gdf_demand[weights_cols].to_numpy() @@ -193,7 +295,7 @@ def from_geodataframe( def facility_client_array(self) -> None: """ - Create an array 2d $m$ x $n$, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the $facility_0$ cover the entire array. + Create an array 2d MxN, where m is number of facilities and n is number of clients. Each row represent a facility and has an array containing clients index meaning that the facility-i cover the entire array. Returns ------- diff --git a/spopt/locate/util.py b/spopt/locate/util.py index 0e201a75..6a4f17f5 100644 --- a/spopt/locate/util.py +++ b/spopt/locate/util.py @@ -20,6 +20,24 @@ def simulated_geo_points(in_data, needed, seed) -> geopandas.GeoDataFrame: ------- geopandas.GeoDataFrame + Examples + -------- + >>> import spaghetti + >>> from spopt.locate.util import simulated_geo_points + + >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) + >>> ntw = spaghetti.Network(in_data=lattice) + >>> street = spaghetti.element_as_gdf(ntw, arcs=True) + + >>> street_buffered = geopandas.GeoDataFrame( + ... geopandas.GeoSeries(street["geometry"].buffer(0.3).unary_union), + ... crs=street.crs, + ... columns=["geometry"]) + + >>> points_simulated = simulated_geo_points(street_buffered, needed=10, seed=1) + >>> type(points_simulated) + + """ geoms = in_data.geometry