Source code for spopt.locate.coverage

# ruff: noqa: B009

import warnings

import numpy as np
import pulp
from geopandas import GeoDataFrame
from scipy.spatial.distance import cdist

from .base import (
    BackupPercentageMixinMixin,
    BaseOutputMixin,
    CoveragePercentageMixin,
    FacilityModelBuilder,
    LocateSolver,
)


[docs] class LSCP(LocateSolver, BaseOutputMixin): r""" Implement the Location Set Covering Problem (*LSCP*) optimization model and solve it. A capacitated version, the Capacitated Location Set Covering Problem – System Optimal (*CLSCP-SO*), can also be solved by passing in client client demand and facility capacities. The standard *LSCP*, as adapted from :cite:`church_murray_2018`, can be formulated as: .. math:: \begin{array}{lllll} \displaystyle \textbf{Minimize} & \displaystyle \sum_{j \in J}{Y_j} && & (1) \\ \displaystyle \textbf{Subject To} & \displaystyle \sum_{j \in N_i}{Y_j} \geq 1 && \forall i \in I & (2) \\ & Y_j \in \{0,1\} && \forall j \in J & (3) \\ & && & \\ \displaystyle \textbf{Where} && i & = & \textrm{index of demand points/areas/objects in set } I \\ && j & = & \textrm{index of potential facility sites in set } J \\ && S & = & \textrm{maximum acceptable service distance or time standard} \\ && d_{ij} & = & \textrm{shortest distance or travel time between locations } i \textrm{ and } j \\ && N_i & = & \{j | d_{ij} < S\} \\ && Y_j & = & \begin{cases} 1, \textrm{if a facility is sited at location } j \\ 0, \textrm{otherwise} \\ \end{cases} \end{array} The *CLSCP-SO*, as adapted from :cite:`church_murray_2018` (see also :cite:`current1988capacitated`), can be formulated as: .. math:: \begin{array}{lllll} \displaystyle \textbf{Minimize} & \displaystyle \sum_{j \in J}{Y_j} && & (1) \\ \displaystyle \textbf{Subject to} & \displaystyle \sum_{j \in N_i}{z_{ij}} = 1 && \forall i \in I & (2) \\ & \displaystyle \sum_{i \in I} a_i z_{ij} \leq C_jx_j && \forall j \in J & (3) \\ & Y_j \in {0,1} && \forall j \in J & (4) \\ & z_{ij} \geq 0 && \forall i \in I \quad \forall j \in N_i & (5) \\ & && & \\ \displaystyle \textbf{Where:} && i & = & \textrm{index of demand points/areas/objects in set } I \\ && j & = & \textrm{index of potential facility sites in set } J \\ && S & = & \textrm{maximal acceptable service distance or time standard} \\ && d_{ij} & = & \textrm{shortest distance or travel time between nodes } i \textrm{ and } j \\ && N_i & = & \{j | d_{ij} < S\} \\ && a_i & = & \textrm{amount of demand at } i \\ && C_j & = & \textrm{capacity of potential facility } j \\ && z_{ij} & = & \textrm{fraction of demand } i \textrm{ that is assigned to facility } j \\ && Y_j & = & \begin{cases} 1, \text{if a facility is located at node } j \\ 0, \text{otherwise} \\ \end{cases} \end{array} Parameters ---------- name : str The problem name. problem : pulp.LpProblem A ``pulp`` instance of an optimization model that contains constraints, variables, and an objective function. Attributes ---------- name : str The problem name. problem : pulp.LpProblem A ``pulp`` instance of an optimization model that contains constraints, variables, and an objective function. fac2cli : numpy.array A 2D array storing facility to client relationships where each row represents a facility and contains an array of client indices with which it is associated. An empty client array indicates the facility is associated with no clients. cli2fac : numpy.array The inverse of ``fac2cli`` where client to facility relationships are shown. aij : numpy.array A cost matrix in the form of a 2D array between origins and destinations. """ # noqa: E501
[docs] def __init__(self, name: str, problem: pulp.LpProblem): super().__init__(name, problem)
def __add_obj(self) -> None: """ Add the objective function to the model. Minimize y1 + y2 + ... + yj Returns ------- None """ fac_vars = getattr(self, "fac_vars") self.problem += pulp.lpSum(fac_vars), "objective function"
[docs] @classmethod def from_cost_matrix( cls, cost_matrix: np.array, service_radius: float, predefined_facilities_arr: np.array = None, demand_quantity_arr: np.array = None, facility_capacity_arr: np.array = None, name: str = "lscp", ): """ Create an ``LSCP`` object based on a cost matrix. A capacitated version of the *LSCP* (the *CLSCP-SO*) can be solved by passing in the ``demand_quantity_arr`` and ``facility_capacity_arr`` keyword arguments. Parameters ---------- cost_matrix : numpy.array A cost matrix in the form of a 2D array between origins and destinations. service_radius : float Maximum acceptable service distance. predefined_facilities_arr : numpy.array (default None) A binary 1D array of service facilities that must appear in the solution. For example, consider 3 facilites ``['A', 'B', 'C']``. If facility ``'B'`` must be in the model solution, then the passed in array should be ``[0, 1, 0]``. demand_quantity_arr : numpy.array (default None) Amount of demand at each client location. facility_capacity_arr : numpy.array (default None) Capacity for service at each facility location. name : str, (default 'lscp') The name of the problem. Returns ------- spopt.locate.coverage.LSCP Examples -------- >>> from spopt.locate import LSCP >>> from spopt.locate.util import simulated_geo_points >>> import geopandas >>> import numpy >>> import pulp >>> import spaghetti Create a regular lattice. >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) >>> ntw = spaghetti.Network(in_data=lattice) >>> streets = spaghetti.element_as_gdf(ntw, arcs=True) >>> streets_buffered = geopandas.GeoDataFrame( ... geopandas.GeoSeries(streets["geometry"].buffer(0.2).unary_union), ... crs=streets.crs, ... columns=["geometry"] ... ) Simulate points about the lattice. >>> demand_points = simulated_geo_points(streets_buffered, needed=100, seed=5) >>> facility_points = simulated_geo_points(streets_buffered, needed=5, seed=6) Snap the points to the network of lattice edges. >>> 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 from origins to destinations. >>> cost_matrix = ntw.allneighbordistances( ... sourcepattern=ntw.pointpatterns["clients"], ... destpattern=ntw.pointpatterns["facilities"]) Create and solve an ``LSCP`` instance from the cost matrix. >>> lscp_from_cost_matrix = LSCP.from_cost_matrix( ... cost_matrix, service_radius=8 ... ) >>> lscp_from_cost_matrix = lscp_from_cost_matrix.solve( ... pulp.PULP_CBC_CMD(msg=False) ... ) Get the facility lookup demand coverage array. >>> for fac, cli in enumerate(lscp_from_cost_matrix.fac2cli): ... print(f"facility {fac} serving {len(cli)} clients") facility 0 serving 0 clients facility 1 serving 63 clients facility 2 serving 85 clients facility 3 serving 0 clients facility 4 serving 58 clients """ n_cli = cost_matrix.shape[0] r_cli = range(n_cli) r_fac = range(cost_matrix.shape[1]) model = pulp.LpProblem(name, pulp.LpMinimize) lscp = LSCP(name, model) if demand_quantity_arr is not None and facility_capacity_arr is None: raise ValueError( "Demand quantities supplied with no facility capacities. " "Model cannot satisfy clients with different " "demands without facility capacities." ) lscp.aij = np.zeros(cost_matrix.shape) lscp.aij[cost_matrix <= service_radius] = 1 if (demand_quantity_arr is None) and (facility_capacity_arr is not None): demand_quantity_arr = np.ones(n_cli) FacilityModelBuilder.add_facility_integer_variable(lscp, r_fac, "y[{i}]") if predefined_facilities_arr is not None: FacilityModelBuilder.add_predefined_facility_constraint( lscp, predefined_facilities_arr ) if demand_quantity_arr is not None: sum_demand = demand_quantity_arr.sum() sum_capacity = facility_capacity_arr.sum() if sum_demand > sum_capacity: raise ValueError( f"Infeasible model. Demand greater than capacity " f"({sum_demand} > {sum_capacity})." ) FacilityModelBuilder.add_client_assign_variable( lscp, r_cli, r_fac, "z[{i}_{j}]", up_bound=None, lp_category=pulp.LpContinuous, ) FacilityModelBuilder.add_facility_capacity_constraint( lscp, demand_quantity_arr, facility_capacity_arr, r_cli, r_fac ) FacilityModelBuilder.add_client_demand_satisfaction_constraint( lscp, r_cli, r_fac ) else: FacilityModelBuilder.add_set_covering_constraint(lscp, r_cli, r_fac) lscp.__add_obj() return lscp
[docs] @classmethod def from_geodataframe( cls, gdf_demand: GeoDataFrame, gdf_fac: GeoDataFrame, demand_col: str, facility_col: str, service_radius: float, predefined_facility_col: str = None, demand_quantity_col: str = None, facility_capacity_col: str = None, distance_metric: str = "euclidean", name: str = "lscp", ): """ Create an ``LSCP`` object from ``geopandas.GeoDataFrame`` objects. Calculate the cost matrix between demand and facility locations before building the problem within the ``from_cost_matrix()`` method. A capacitated version of the *LSCP* (the *CLSCP-SO*) can be solved by passing in the ``demand_quantity_col`` and ``facility_capacity_col`` keyword arguments. Parameters ---------- gdf_demand : geopandas.GeoDataFrame Demand locations. gdf_fac : geopandas.GeoDataFrame Facility locations. demand_col : str Demand sites geometry column name. facility_col : str Facility candidate sites geometry column name. service_radius : float Maximum acceptable service distance. predefined_facility_col : str (default None) Column name representing facilities are already defined. This a binary assignment per facility. For example, consider 3 facilites ``['A', 'B', 'C']``. If facility ``'B'`` must be in the model solution, then the column should be ``[0, 1, 0]``. demand_quantity_col : str Column name representing amount of demand at each client location. facility_capacity_arr : str Column name representing capacity for service at each facility location. distance_metric : str (default 'euclidean') A metric used for the distance calculations supported by `scipy.spatial.distance.cdist <https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html>`_. name : str (default 'lscp') The name of the problem. Returns ------- spopt.locate.coverage.LSCP Examples -------- >>> from spopt.locate import LSCP >>> from spopt.locate.util import simulated_geo_points >>> import geopandas >>> import numpy >>> import pulp >>> import spaghetti Create a regular lattice. >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) >>> ntw = spaghetti.Network(in_data=lattice) >>> streets = spaghetti.element_as_gdf(ntw, arcs=True) >>> streets_buffered = geopandas.GeoDataFrame( ... geopandas.GeoSeries(streets["geometry"].buffer(0.2).unary_union), ... crs=streets.crs, ... columns=["geometry"] ... ) Simulate points about the lattice. >>> demand_points = simulated_geo_points(streets_buffered, needed=100, seed=5) >>> facility_points = simulated_geo_points(streets_buffered, needed=5, seed=6) Snap the points to the network of lattice edges and extract as ``GeoDataFrame`` objects. >>> 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 and solve an ``LSCP`` instance from the ``GeoDataFrame`` objects. >>> lscp_from_geodataframe = LSCP.from_geodataframe( ... clients_snapped, ... facilities_snapped, ... "geometry", ... "geometry", ... service_radius=8, ... distance_metric="euclidean" ... ) >>> lscp_from_geodataframe = lscp_from_geodataframe.solve( ... pulp.PULP_CBC_CMD(msg=False) ... ) Get the facility lookup demand coverage array. >>> for fac, cli in enumerate(lscp_from_geodataframe.fac2cli): ... print(f"facility {fac} serving {len(cli)} clients") facility 0 serving 0 clients facility 1 serving 0 clients facility 2 serving 100 clients facility 3 serving 0 clients facility 4 serving 0 clients """ # noqa: E501 demand_quantity_arr = None if demand_quantity_col is not None: demand_quantity_arr = gdf_demand[demand_quantity_col].to_numpy() facility_capacity_arr = None if facility_capacity_col is not None: facility_capacity_arr = gdf_fac[facility_capacity_col].to_numpy() if demand_quantity_arr is not None and facility_capacity_arr is None: raise ValueError( "Demand quantities supplied with no facility capacities. " "Model cannot satisfy clients with different " "demands without facility capacities." ) predefined_facilities_arr = None if predefined_facility_col is not None: predefined_facilities_arr = gdf_fac[predefined_facility_col].to_numpy() dem = gdf_demand[demand_col] fac = gdf_fac[facility_col] dem_type_geom = dem.geom_type.unique() fac_type_geom = fac.geom_type.unique() _msg = ( " geodataframe contains mixed type geometries or is not a point. Be " "sure deriving centroid from geometries doesn't affect the results." ) if len(dem_type_geom) > 1 or "Point" not in dem_type_geom: warnings.warn(f"Demand{_msg}", UserWarning, stacklevel=2) dem = dem.centroid if len(fac_type_geom) > 1 or "Point" not in fac_type_geom: warnings.warn(f"Facility{_msg}", UserWarning, stacklevel=2) fac = fac.centroid dem_data = np.array([dem.x.to_numpy(), dem.y.to_numpy()]).T fac_data = np.array([fac.x.to_numpy(), fac.y.to_numpy()]).T if gdf_demand.crs != gdf_fac.crs: raise ValueError( "Geodataframes crs are different: " f"gdf_demand-{gdf_demand.crs}, gdf_fac-{gdf_fac.crs}" ) distances = cdist(dem_data, fac_data, distance_metric) return cls.from_cost_matrix( distances, service_radius, predefined_facilities_arr, demand_quantity_arr, facility_capacity_arr, name=("capacitated-" + name if facility_capacity_col is not None else name), )
[docs] def facility_client_array(self) -> None: """ Create a 2D array storing **facility to client relationships** where each row represents a facility and contains an array of client indices with which it is associated. An empty client array indicates the facility is associated with no clients. Returns ------- None """ fac_vars = getattr(self, "fac_vars") len_fac_vars = len(fac_vars) self.fac2cli = [] for j in range(len_fac_vars): array_cli = [] if fac_vars[j].value() > 0: for i in range(self.aij.shape[0]): if self.aij[i, j] > 0: array_cli.append(i) self.fac2cli.append(array_cli)
[docs] def solve(self, solver: pulp.LpSolver, results: bool = True): """ Solve the ``LSCP`` model. Parameters ---------- solver : pulp.apis.LpSolver A solver supported by ``pulp``. results : bool (default True) If ``True`` it will create metainfo (which facilities cover which demand) and vice-versa, and the uncovered demand. Returns ------- spopt.locate.coverage.LSCP """ self.problem.solve(solver) self.check_status() if results: self.facility_client_array() self.client_facility_array() return self
[docs] class LSCPB(LocateSolver, BaseOutputMixin, BackupPercentageMixinMixin): r""" Implement the Location Set Covering Problem - Backup (*LSCP-B*) optimization model and solve it. The *LSCP-B*, as adapted from :cite:`church_murray_2018`, can be formulated as: .. math:: \begin{array}{lllll} \displaystyle \textbf{Maximize} & \displaystyle \sum_{i \in I}{U_i} && & (1) \\ \displaystyle \textbf{Subject To} & \displaystyle \sum_{j \in J}{a_{ij}}{Y_j} \geq 1 + U_i && \forall i \in I & (2) \\ & \displaystyle \sum_{j \in J}{Y_j} = p && & (3) \\ & U_i \leq 1 && \forall i \in I & (4) \\ & Y_j \in \{0, 1\} && \forall j \in J & (5) \\ & && & \\ \displaystyle \textbf{Where} && i & = & \textrm{index of demand points/areas/objects in set } I \\ && j & = & \textrm{index of potential facility sites in set } J \\ && p & = & \textrm{objective value identified by using the } LSCP \\ && U_i & = & \begin{cases} 1, \textrm{if demand location is covered twice} \\ 0, \textrm{if demand location is covered once} \\ \end{cases} \\ && a_{ij} & = & \begin{cases} 1, \textrm{if facility location } j \textrm{ covers demand location } i \\ 0, \textrm{otherwise} \\ \end{cases} \\ && Y_j & = & \begin{cases} 1, \textrm{if a facility is sited at location } j \\ 0, \textrm{otherwise} \\ \end{cases} \end{array} Parameters ---------- name : str The problem name. problem : pulp.LpProblem A ``pulp`` instance of an optimization model that contains constraints, variables, and an objective function. solver : pulp.LpSolver A solver supported by ``pulp``. Attributes ---------- name : str The problem name. problem : pulp.LpProblem A ``pulp`` instance of an optimization model that contains constraints, variables, and an objective function. solver : pulp.LpSolver A solver supported by ``pulp``. lscp_obj_value : float The objective value returned from a solved ``LSCP`` instance. fac2cli : numpy.array A 2D array storing facility to client relationships where each row represents a facility and contains an array of client indices with which it is associated. An empty client array indicates the facility is associated with no clients. cli2fac : numpy.array The inverse of ``fac2cli`` where client to facility relationships are shown. aij : numpy.array A cost matrix in the form of a 2D array between origins and destinations. """ # noqa: E501
[docs] def __init__( self, name: str, problem: pulp.LpProblem, solver: pulp.LpSolver, ): self.solver = solver super().__init__(name, problem)
def __add_obj(self) -> None: """ Add the objective function to the model. Maximize U1 + U2 + ... + Uj Returns ------- None """ cov_vars = getattr(self, "cli_vars") self.problem += pulp.lpSum(cov_vars), "objective function"
[docs] @classmethod def from_cost_matrix( cls, cost_matrix: np.array, service_radius: float, solver: pulp.LpSolver, predefined_facilities_arr: np.array = None, name: str = "lscp-b", ): """ Create an ``LSCPB`` object based on a cost matrix. Parameters ---------- cost_matrix : numpy.array A cost matrix in the form of a 2D array between origins and destinations. service_radius : float Maximum acceptable service distance. solver : pulp.LpSolver A solver supported by ``pulp``. predefined_facilities_arr : numpy.array (default None) A binary 1D array of service facilities that must appear in the solution. For example, consider 3 facilites ``['A', 'B', 'C']``. If facility ``'B'`` must be in the model solution, then the passed in array should be ``[0, 1, 0]``. name : str (default 'lscp-b') The problem name. Returns ------- spopt.locate.coverage.LSCPB Examples -------- >>> from spopt.locate import LSCPB >>> from spopt.locate.util import simulated_geo_points >>> import geopandas >>> import numpy >>> import pulp >>> import spaghetti Create a regular lattice. >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) >>> ntw = spaghetti.Network(in_data=lattice) >>> streets = spaghetti.element_as_gdf(ntw, arcs=True) >>> streets_buffered = geopandas.GeoDataFrame( ... geopandas.GeoSeries(streets["geometry"].buffer(0.2).unary_union), ... crs=streets.crs, ... columns=["geometry"] ... ) Simulate points about the lattice. >>> demand_points = simulated_geo_points(streets_buffered, needed=100, seed=5) >>> facility_points = simulated_geo_points(streets_buffered, needed=5, seed=6) Snap the points to the network of lattice edges. >>> 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 from origins to destinations. >>> cost_matrix = ntw.allneighbordistances( ... sourcepattern=ntw.pointpatterns["clients"], ... destpattern=ntw.pointpatterns["facilities"] ... ) Create and solve an ``LSCPB`` instance from the cost matrix. >>> lscpb_from_cost_matrix = LSCPB.from_cost_matrix( ... cost_matrix, service_radius=8, solver=pulp.PULP_CBC_CMD(msg=False) ... ) >>> lscpb_from_cost_matrix = lscpb_from_cost_matrix.solve() Get the facility lookup demand coverage array. >>> for fac, cli in enumerate(lscpb_from_cost_matrix.fac2cli): ... print(f"facility {fac} serving {len(cli)} clients") facility 0 serving 0 clients facility 1 serving 63 clients facility 2 serving 85 clients facility 3 serving 92 clients facility 4 serving 0 clients Get the percentage of clients covered by more than one facility. >>> round(lscpb_from_cost_matrix.backup_perc, 3) 88.0 88% of clients are covered by more than 1 facility """ if predefined_facilities_arr is not None: lscp = LSCP.from_cost_matrix( cost_matrix, service_radius, predefined_facilities_arr ) else: lscp = LSCP.from_cost_matrix(cost_matrix, service_radius) lscp.solve(solver) r_cli = range(cost_matrix.shape[0]) r_fac = range(cost_matrix.shape[1]) model = pulp.LpProblem(name, pulp.LpMaximize) lscpb = LSCPB(name, model, solver) lscpb.lscp_obj_value = lscp.problem.objective.value() FacilityModelBuilder.add_facility_integer_variable(lscpb, r_fac, "y[{i}]") FacilityModelBuilder.add_client_integer_variable(lscpb, r_cli, "x[{i}]") lscpb.aij = np.zeros(cost_matrix.shape) lscpb.aij[cost_matrix <= service_radius] = 1 if predefined_facilities_arr is not None: FacilityModelBuilder.add_predefined_facility_constraint( lscpb, predefined_facilities_arr ) lscpb.__add_obj() FacilityModelBuilder.add_facility_constraint(lscpb, lscpb.lscp_obj_value) FacilityModelBuilder.add_backup_covering_constraint(lscpb, r_fac, r_cli) return lscpb
[docs] @classmethod def from_geodataframe( cls, gdf_demand: GeoDataFrame, gdf_fac: GeoDataFrame, demand_col: str, facility_col: str, service_radius: float, solver: pulp.LpSolver, predefined_facility_col: str = None, distance_metric: str = "euclidean", name: str = "lscp-b", ): """ Create an ``LSCPB`` object from ``geopandas.GeoDataFrame`` objects. Calculate the cost matrix between demand and facility locations before building the problem within the ``from_cost_matrix()`` method. Parameters ---------- gdf_demand : geopandas.GeoDataFrame Demand locations. gdf_fac : geopandas.GeoDataFrame Facility locations. demand_col : str Demand sites geometry column name. facility_col : str Facility candidate sites geometry column name. service_radius : float Maximum acceptable service distance. solver : pulp.LpSolver A solver supported by ``pulp``. predefined_facility_col : str (default None) Column name representing facilities are already defined. This a binary assignment per facility. For example, consider 3 facilites ``['A', 'B', 'C']``. If facility ``'B'`` must be in the model solution, then the column should be ``[0, 1, 0]``. distance_metric : str (default 'euclidean') A metric used for the distance calculations supported by `scipy.spatial.distance.cdist <https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html>`_. name : str (default 'lscp-b') The name of the problem. Returns ------- spopt.locate.coverage.LSCPB Examples -------- >>> from spopt.locate import LSCPB >>> from spopt.locate.util import simulated_geo_points >>> import geopandas >>> import numpy >>> import pulp >>> import spaghetti Create a regular lattice. >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) >>> ntw = spaghetti.Network(in_data=lattice) >>> streets = spaghetti.element_as_gdf(ntw, arcs=True) >>> streets_buffered = geopandas.GeoDataFrame( ... geopandas.GeoSeries(streets["geometry"].buffer(0.2).unary_union), ... crs=streets.crs, ... columns=["geometry"] ... ) Simulate points about the lattice. >>> demand_points = simulated_geo_points(streets_buffered, needed=100, seed=5) >>> facility_points = simulated_geo_points(streets_buffered, needed=5, seed=6) Snap the points to the network of lattice edges and extract as ``GeoDataFrame`` objects. >>> 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 and solve an ``LSCPB`` instance from the ``GeoDataFrame`` objects. >>> lscpb_from_geodataframe = LSCPB.from_geodataframe( ... clients_snapped, ... facilities_snapped, ... "geometry", ... "geometry", ... service_radius=8, ... solver=pulp.PULP_CBC_CMD(msg=False), ... distance_metric="euclidean" ... ) >>> lscpb_from_geodataframe = lscpb_from_geodataframe.solve() Get the facility lookup demand coverage array. >>> for fac, cli in enumerate(lscpb_from_geodataframe.fac2cli): ... print(f"facility {fac} serving {len(cli)} clients") facility 0 serving 0 clients facility 1 serving 0 clients facility 2 serving 100 clients facility 3 serving 0 clients facility 4 serving 0 clients Get the percentage of clients covered by more than one facility. >>> round(lscpb_from_geodataframe.backup_perc, 3) 0.0 All clients are covered by 1 facility because only one facility is needed to solve the LSCP. """ # noqa: E501 predefined_facilities_arr = None if predefined_facility_col is not None: predefined_facilities_arr = gdf_fac[predefined_facility_col].to_numpy() dem = gdf_demand[demand_col] fac = gdf_fac[facility_col] dem_type_geom = dem.geom_type.unique() fac_type_geom = fac.geom_type.unique() _msg = ( " geodataframe contains mixed type geometries or is not a point. Be " "sure deriving centroid from geometries doesn't affect the results." ) if len(dem_type_geom) > 1 or "Point" not in dem_type_geom: warnings.warn(f"Demand{_msg}", UserWarning, stacklevel=2) dem = dem.centroid if len(fac_type_geom) > 1 or "Point" not in fac_type_geom: warnings.warn(f"Facility{_msg}", UserWarning, stacklevel=2) fac = fac.centroid dem_data = np.array([dem.x.to_numpy(), dem.y.to_numpy()]).T fac_data = np.array([fac.x.to_numpy(), fac.y.to_numpy()]).T if gdf_demand.crs != gdf_fac.crs: raise ValueError( "Geodataframes crs are different: " f"gdf_demand-{gdf_demand.crs}, gdf_fac-{gdf_fac.crs}" ) distances = cdist(dem_data, fac_data, distance_metric) return cls.from_cost_matrix( distances, service_radius, solver, predefined_facilities_arr, name )
[docs] def facility_client_array(self) -> None: """ Create a 2D array storing **facility to client relationships** where each row represents a facility and contains an array of client indices with which it is associated. An empty client array indicates the facility is associated with no clients. Returns ------- None """ fac_vars = getattr(self, "fac_vars") len_fac_vars = len(fac_vars) self.fac2cli = [] for j in range(len_fac_vars): array_cli = [] if fac_vars[j].value() > 0: for i in range(self.aij.shape[0]): if self.aij[i, j] > 0: array_cli.append(i) self.fac2cli.append(array_cli)
[docs] def solve(self, results: bool = True): """ Solve the ``LSCPB`` model. Parameters ---------- results : bool (default True) If ``True`` it will create metainfo (which facilities cover which demand) and vice-versa, and the uncovered demand. Returns ------- spopt.locate.coverage.LSCPB """ self.problem.solve(self.solver) self.check_status() if results: self.facility_client_array() self.client_facility_array() self.get_percentage() return self
[docs] class MCLP(LocateSolver, BaseOutputMixin, CoveragePercentageMixin): r""" Implement the Maximal Coverage Location Problem (*MCLP*) optimization model and solve it. The *MCLP*, as adapted from :cite:`church_murray_2018`, can be formulated as: .. math:: \begin{array}{lllll} \displaystyle \textbf{Maximize} & \displaystyle \sum_{i \in I}{a_iX_i} && & (1) \\ \displaystyle \textbf{Subject To} & \displaystyle \sum_{j \in N_i}{Y_j \geq X_i} && \forall i \in I & (2) \\ & \displaystyle \sum_{j \in J}{Y_j} = p && & (3) \\ & X_i \in \{0, 1\} && \forall i \in I & (4) \\ & Y_j \in \{0, 1\} && \forall j \in J & (5) \\ & && & \\ \displaystyle \textbf{Where} && i & = & \textrm{index of demand points/areas/objects in set } I \\ && j & = & \textrm{index of potential facility sites in set } J \\ && p & = & \textrm{the number of facilities to be sited} \\ && S & = & \textrm{maximum acceptable service distance or time standard} \\ && d_{ij} & = & \textrm{shortest distance or travel time between locations } i \textrm{ and } j \\ && N_i & = & \{j | d_{ij} < S\} \\ && X_i & = & \begin{cases} 1, \textrm{if client location } i \textrm{ is covered within service standard } S \\ 0, \textrm{otherwise} \\ \end{cases} \\ && Y_j & = & \begin{cases} 1, \textrm{if a facility is sited at location } j \\ 0, \textrm{otherwise} \\ \end{cases} \end{array} Parameters ---------- name : str The problem name. problem : pulp.LpProblem A ``pulp`` instance of an optimization model that contains constraints, variables, and an objective function. Attributes ---------- name : str The problem name. problem : pulp.LpProblem A ``pulp`` instance of an optimization model that contains constraints, variables, and an objective function. fac2cli : numpy.array A 2D array storing facility to client relationships where each row represents a facility and contains an array of client indices with which it is associated. An empty client array indicates the facility is associated with no clients. cli2fac : numpy.array The inverse of ``fac2cli`` where client to facility relationships are shown. aij : numpy.array A cost matrix in the form of a 2D array between origins and destinations. n_cli_uncov : int The number of uncovered client locations. """ # noqa: E501
[docs] def __init__(self, name: str, problem: pulp.LpProblem): super().__init__(name, problem)
def __add_obj(self, weights: np.array, range_clients: range) -> None: """ Add the objective function to the model. Maximize w1 * y1 + w2 * y2 + ... + wi * yi Returns ------- None """ dem_vars = getattr(self, "cli_vars") self.problem += ( pulp.lpSum([weights.flatten()[i] * dem_vars[i] for i in range_clients]), "objective function", )
[docs] @classmethod def from_cost_matrix( cls, cost_matrix: np.array, weights: np.array, service_radius: float, p_facilities: int, predefined_facilities_arr: np.array = None, name: str = "mclp", ): """ Create a ``MCLP`` object based on cost matrix. Parameters ---------- cost_matrix : numpy.array A cost matrix in the form of a 2D array between origins and destinations. weights : numpy.array A 1D array of service load or population demand. service_radius : float Maximum acceptable service distance. p_facilities : int The number of facilities to be located. predefined_facilities_arr : numpy.array (default None) A binary 1D array of service facilities that must appear in the solution. For example, consider 3 facilites ``['A', 'B', 'C']``. If facility ``'B'`` must be in the model solution, then the passed in array should be ``[0, 1, 0]``. name : str (default 'mclp') The problem name. Returns ------- spopt.locate.coverage.MCLP Examples -------- >>> from spopt.locate import MCLP >>> from spopt.locate.util import simulated_geo_points >>> import geopandas >>> import numpy >>> import pulp >>> import spaghetti Create a regular lattice. >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) >>> ntw = spaghetti.Network(in_data=lattice) >>> streets = spaghetti.element_as_gdf(ntw, arcs=True) >>> streets_buffered = geopandas.GeoDataFrame( ... geopandas.GeoSeries(streets["geometry"].buffer(0.2).unary_union), ... crs=streets.crs, ... columns=["geometry"] ... ) Simulate points about the lattice. >>> demand_points = simulated_geo_points(streets_buffered, needed=100, seed=5) >>> facility_points = simulated_geo_points(streets_buffered, needed=5, seed=6) Snap the points to the network of lattice edges. >>> 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 from origins to destinations. >>> 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 and solve an ``MCLP`` instance from the cost matrix. >>> mclp_from_cost_matrix = MCLP.from_cost_matrix( ... cost_matrix, ai, service_radius=7, p_facilities=4 ... ) >>> mclp_from_cost_matrix = mclp_from_cost_matrix.solve( ... pulp.PULP_CBC_CMD(msg=False) ... ) Get the facility lookup demand coverage array. >>> for fac, cli in enumerate(mclp_from_cost_matrix.fac2cli): ... print(f"facility {fac} serving {len(cli)} clients") facility 0 serving 44 clients facility 1 serving 54 clients facility 2 serving 75 clients facility 3 serving 77 clients facility 4 serving 0 clients Get the number of clients uncovered and percentage covered. >>> mclp_from_cost_matrix.n_cli_uncov 1 >>> mclp_from_cost_matrix.perc_cov 99.0 """ n_cli = cost_matrix.shape[0] r_cli = range(n_cli) r_fac = range(cost_matrix.shape[1]) model = pulp.LpProblem(name, pulp.LpMaximize) mclp = MCLP(name, model) FacilityModelBuilder.add_facility_integer_variable(mclp, r_fac, "x[{i}]") FacilityModelBuilder.add_client_integer_variable(mclp, r_cli, "y[{i}]") mclp.aij = np.zeros(cost_matrix.shape) mclp.aij[cost_matrix <= service_radius] = 1 weights = np.reshape(weights, (n_cli, 1)) mclp.__add_obj(weights, r_cli) if predefined_facilities_arr is not None: FacilityModelBuilder.add_predefined_facility_constraint( mclp, predefined_facilities_arr ) FacilityModelBuilder.add_maximal_coverage_constraint(mclp, r_fac, r_cli) FacilityModelBuilder.add_facility_constraint(mclp, p_facilities) return mclp
[docs] @classmethod def from_geodataframe( cls, gdf_demand: GeoDataFrame, gdf_fac: GeoDataFrame, demand_col: str, facility_col: str, weights_cols: str, service_radius: float, p_facilities: int, predefined_facility_col: str = None, distance_metric: str = "euclidean", name: str = "mclp", ): """ Create an ``MCLP`` object from ``geopandas.GeoDataFrame`` objects. Calculate the cost matrix between demand and facility locations before building the problem within the ``from_cost_matrix()`` method. Parameters ---------- gdf_demand : geopandas.GeoDataFrame Demand locations. gdf_fac : geopandas.GeoDataFrame Facility locations. demand_col : str Demand sites geometry column name. facility_col : str Facility candidate sites geometry column name. weights_cols : str The weight column name representing service load or demand. service_radius : float Maximum acceptable service distance. p_facilities: int The number of facilities to be located. predefined_facility_col : str (default None) Column name representing facilities are already defined. This a binary assignment per facility. For example, consider 3 facilites ``['A', 'B', 'C']``. If facility ``'B'`` must be in the model solution, then the column should be ``[0, 1, 0]``. distance_metric : str (default 'euclidean') A metric used for the distance calculations supported by `scipy.spatial.distance.cdist <https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html>`_. name : str (default 'mclp') The name of the problem. Returns ------- spopt.locate.coverage.MCLP Examples -------- >>> from spopt.locate import MCLP >>> from spopt.locate.util import simulated_geo_points >>> import geopandas >>> import pulp >>> import numpy >>> import spaghetti Create a regular lattice. >>> lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True) >>> ntw = spaghetti.Network(in_data=lattice) >>> streets = spaghetti.element_as_gdf(ntw, arcs=True) >>> streets_buffered = geopandas.GeoDataFrame( ... geopandas.GeoSeries(streets["geometry"].buffer(0.2).unary_union), ... crs=streets.crs, ... columns=["geometry"] ... ) Simulate points about the lattice. >>> demand_points = simulated_geo_points(streets_buffered, needed=100, seed=5) >>> facility_points = simulated_geo_points(streets_buffered, needed=5, seed=6) Snap the points to the network of lattice edges and extract as ``GeoDataFrame`` objects. >>> 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 and solve an ``MCLP`` instance from the ``GeoDataFrame`` objects. >>> mclp_from_geodataframe = MCLP.from_geodataframe( ... clients_snapped, ... facilities_snapped, ... "geometry", ... "geometry", ... "weights", ... service_radius=7, ... p_facilities=4, ... distance_metric="euclidean" ... ) >>> mclp_from_geodataframe = mclp_from_geodataframe.solve( ... pulp.PULP_CBC_CMD(msg=False) ... ) Get the facility lookup demand coverage array. >>> for fac, cli in enumerate(mclp_from_geodataframe.fac2cli): ... print(f"facility {fac} serving {len(cli)} clients") facility 0 serving 63 clients facility 1 serving 80 clients facility 2 serving 96 clients facility 3 serving 0 clients facility 4 serving 58 clients Get the number of clients uncovered and percentage covered. >>> mclp_from_geodataframe.n_cli_uncov 0 >>> mclp_from_geodataframe.perc_cov 100.0 """ # noqa E501 predefined_facilities_arr = None if predefined_facility_col is not None: predefined_facilities_arr = gdf_fac[predefined_facility_col].to_numpy() service_load = gdf_demand[weights_cols].to_numpy() dem = gdf_demand[demand_col] fac = gdf_fac[facility_col] dem_type_geom = dem.geom_type.unique() fac_type_geom = fac.geom_type.unique() _msg = ( " geodataframe contains mixed type geometries or is not a point. Be " "sure deriving centroid from geometries doesn't affect the results." ) if len(dem_type_geom) > 1 or "Point" not in dem_type_geom: warnings.warn(f"Demand{_msg}", UserWarning, stacklevel=2) dem = dem.centroid if len(fac_type_geom) > 1 or "Point" not in fac_type_geom: warnings.warn(f"Facility{_msg}", UserWarning, stacklevel=2) fac = fac.centroid dem_data = np.array([dem.x.to_numpy(), dem.y.to_numpy()]).T fac_data = np.array([fac.x.to_numpy(), fac.y.to_numpy()]).T if gdf_demand.crs != gdf_fac.crs: raise ValueError( "Geodataframes crs are different: " f"gdf_demand-{gdf_demand.crs}, gdf_fac-{gdf_fac.crs}" ) distances = cdist(dem_data, fac_data, distance_metric) return cls.from_cost_matrix( distances, service_load, service_radius, p_facilities, predefined_facilities_arr, name, )
[docs] def facility_client_array(self) -> None: """ Create a 2D array storing **facility to client relationships** where each row represents a facility and contains an array of client indices with which it is associated. An empty client array indicates the facility is associated with no clients. Returns ------- None """ fac_vars = getattr(self, "fac_vars") cli_vars = getattr(self, "cli_vars") len_fac_vars = len(fac_vars) self.fac2cli = [] for j in range(len_fac_vars): array_cli = [] if fac_vars[j].value() > 0: for i in range(self.aij.shape[0]): if cli_vars[i].value() > 0 and self.aij[i, j] > 0: array_cli.append(i) self.fac2cli.append(array_cli)
[docs] def solve(self, solver: pulp.LpSolver, results: bool = True): """ Solve the ``MCLP`` model Parameters ---------- solver : pulp.LpSolver A solver supported by ``pulp``. results : bool (default True) If ``True`` it will create metainfo (which facilities cover which demand) and vice-versa, and the uncovered demand. Returns ------- spopt.locate.coverage.MCLP """ self.problem.solve(solver) self.check_status() if results: self.facility_client_array() self.client_facility_array() self.uncovered_clients() self.get_percentage() return self