import spopt, libpysal, geopandas, numpy
mexico = geopandas.read_file(libpysal.examples.get_path("mexicojoin.shp"))
mexico["count"] = 1
attrs = [f"PCGDP{year}" for year in range(1950, 2010, 10)]
w = libpysal.weights.Queen.from_dataframe(mexico)
mexico["count"], threshold_name, threshold, top_n = 1, "count", 4, 2
numpy.random.seed(123456)
model = spopt.region.MaxPHeuristic(mexico, w, attrs, threshold_name, threshold, top_n)
model.solve()
mexico["maxp_new"] = model.labels_
mexico.plot(column="maxp_new", categorical=True, figsize=(12,8), ec="w");

from spopt.locate import MCLP
from spopt.locate.util import simulated_geo_points
import numpy, geopandas, pulp, spaghetti

solver = pulp.PULP_CBC_CMD(msg=False, warmStart=True)
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.5).unary_union),
    crs=street.crs,
    columns=["geometry"],
)
client_points = simulated_geo_points(street_buffered, needed=100, seed=5)
ntw.snapobservations(client_points, "clients", attribute=True)
clients_snapped = spaghetti.element_as_gdf(
    ntw, pp_name="clients", snapped=True
)
facility_points = simulated_geo_points(street_buffered, needed=10, seed=6)
ntw.snapobservations(facility_points, "facilities", attribute=True)
facilities_snapped = spaghetti.element_as_gdf(
    ntw, pp_name="facilities", snapped=True
)
cost_matrix = ntw.allneighbordistances(
    sourcepattern=ntw.pointpatterns["clients"],
    destpattern=ntw.pointpatterns["facilities"],
)
numpy.random.seed(0)
ai = numpy.random.randint(1, 12, 100)
mclp_from_cost_matrix = MCLP.from_cost_matrix(cost_matrix, ai, 4, p_facilities=4)
mclp_from_cost_matrix = mclp_from_cost_matrix.solve(solver)