diff --git a/spopt/route/__init__.py b/spopt/route/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/spopt/route/engine.py b/spopt/route/engine.py new file mode 100644 index 00000000..636c8551 --- /dev/null +++ b/spopt/route/engine.py @@ -0,0 +1,224 @@ +try: + import osrm + has_bindings = True +except (ImportError,ModuleNotFoundError) as e: + has_bindings = False +import os +import numpy +import requests +import warnings +import geopandas +import shapely +from sklearn import metrics + +# TODO: needs to be configurable by site +_OSRM_DATABASE_FILE = "" + +def build_route_table(demand_sites, candidate_depots, cost='distance', http=not has_bindings, database_path=_OSRM_DATABASE_FILE, port=5000): + """ + Build a route table using OSRM, either over http or over py-osrm bindings + """ + if isinstance(demand_sites, (geopandas.GeoSeries, geopandas.GeoDataFrame)): + demand_sites = demand_sites.geometry.get_coordinates().values + if isinstance(candidate_depots, (geopandas.GeoSeries, geopandas.GeoDataFrame)): + candidate_depots = candidate_depots.geometry.get_coordinates().values + if cost not in ("distance", "duration", "both"): + raise ValueError(f"cost option '{cost}' not one of the supported options, ('distance', 'duration', 'both')") + if http: + try: + distances, durations = _build_route_table_http(demand_sites, candidate_depots, cost=cost, port=port) + except (requests.ConnectionError, requests.JSONDecodeError): + warnings.warn( + "Failed to connect to routing engine... using haversine distance" + " and (d/500)**.75 for durations" + ) + distances = metrics.pairwise_distances( + numpy.fliplr(numpy.deg2rad(demand_sites)), + numpy.fliplr(numpy.deg2rad(candidate_depots)), + metric="haversine" + ) * 6371000 + durations = numpy.ceil((distances / 10) ** .75) + else: + distances, durations = _build_route_table_pyosrm( + demand_sites, candidate_depots, database_path=database_path + ) + for D in (distances, durations): + if D is None: + continue + n_row, n_col = D.shape + assert n_row == len(candidate_depots) + assert n_col == len(demand_sites) + no_route_available = numpy.isnan(D) + D[no_route_available] = D[~no_route_available].sum() + if cost == 'distance': + return distances + elif cost == 'duration': + return durations + elif cost == 'both': + return distances, durations + +def build_specific_route(waypoints, port=5000, http=not has_bindings, return_durations=True, database_path=_OSRM_DATABASE_FILE): + """ + Build a route over the road network from each waypoint to each other waypoint. If the routing engine is not found, this builds straight-line + routes, and measures their duration as a nonlinear function of the + haversine distance between input points. + """ + if isinstance(waypoints, (geopandas.GeoSeries, geopandas.GeoDataFrame)): + waypoints = waypoints.geometry.get_coordinates().values + if http: + try: + out = _build_specific_route_http(waypoints, port=port, return_durations=return_durations) + except (requests.ConnectionError, requests.JSONDecodeError): + warnings.warn( + "Failed to connect to routing engine... constructed routes" + " will be straight lines and may not follow the road network." + ) + route = shapely.LineString(waypoints) + prep_points = numpy.fliplr(numpy.deg2rad(waypoints)) + durations = [ + (metrics.pairwise.haversine_distances([prep_points[i]], [prep_points[i+1]]) + * 637000 / 10)**.75 + for i in range(len(prep_points)-1) + ] + out = (route, durations) if return_durations else route + else: + route = _build_specific_route_pyosrm(waypoints, database_path=database_path, return_durations=return_durations) + if return_durations: + route, durations = out + return route, durations + else: + route = out + return route + +def _build_specific_route_http(waypoints, return_durations=True, port=5000): + + # TODO: needs to be configurable by site + baseurl = f"http://127.0.0.1:{int(port)}/route/v1/driving/" + + point_string = ";".join( + map( + lambda x: "{},{}".format(*x), + waypoints, + ) + ) + + request_url = ( + baseurl + + point_string + + "?" + + "steps=true" + + "&" + + f"geometries=geojson" + + "&" + + "annotations=true" + ) + routes = requests.get(request_url).json()['routes'] + assert len(routes) == 1 + route = routes[0] + #sub_coordinates = numpy.empty(shape=(0,2)) + route_shape = shapely.geometry.shape(route['geometry']) + leg_durations = numpy.array([leg['duration'] for leg in route['legs']]) + """ + for leg_i, leg in enumerate(route['legs']): + durations[i] = leg['duration'] + for steps in leg['steps']: + assert steps['geometry']['type'] == "LineString" + sub_coordinates = numpy.row_stack((sub_coordinates, + numpy.asarray(steps['geometry']['coordinates'])[:-1] + )) + """ + #route_shape = shapely.LineString(sub_coordinates) + numpy.testing.assert_array_equal( + shapely.get_num_geometries(route_shape), + numpy.ones((len(waypoints),)) + ) + if return_durations: + return route_shape, leg_durations + else: + return route_shape + +def _build_specific_route_pyosrm(waypoints, database_path=_OSRM_DATABASE_FILE, return_durations=False): + raise NotImplementedError() + +def _build_route_table_http(demand_sites, candidate_depots, cost='distance', port=5000): + """ + Build a route table using the http interface to the OSRM engine + """ + request_url = _create_route_request(demand_sites, candidate_depots, cost=cost, port=port) + request = requests.get(request_url) + content = request.json() + if cost == 'distance': + D = numpy.asarray(content["distances"]).astype(float) + output = (D,None) + elif cost == 'duration': + D = numpy.asarray(content["durations"]).astype(float) + output = (None,D) + elif cost == 'both': + distances = numpy.asarray(content["distances"]).astype(float) + durations = numpy.asarray(content["durations"]).astype(float) + output = (distances, durations) + else: + raise ValueError(f"cost option '{cost}' not one of the supported options, ('distance', 'duration', 'both')") + return output + + +def _create_route_request(demand_sites, candidate_depots, cost='distance', port=5000): + point_string = ";".join( + map( + lambda x: "{},{}".format(*x), + numpy.row_stack((candidate_depots, demand_sites)), + ) + ) + n_demands = len(demand_sites) + n_supplys = len(candidate_depots) + source_string = "sources=" + ";".join(numpy.arange(n_supplys).astype(str)) + destination_string = "destinations=" + ";".join( + numpy.arange(n_supplys, n_demands + n_supplys).astype(str) + ) + # TODO: needs to be configurable by site + baseurl = f"http://127.0.0.1:{int(port)}/table/v1/driving/" + if cost=='distance': + annotation = "&annotations=distance" + elif cost=='duration': + annotation = "&annotations=duration" + elif cost=='both': + annotation = "&annotations=duration,distance" + else: + annotation = "" + + request_url = ( + baseurl + + point_string + + "?" + + source_string + + "&" + + destination_string + + annotation + + "&exclude=ferry" + ) + return request_url + + +def _build_route_table_pyosrm(demand_sites, candidate_depots, database_path=_OSRM_DATABASE_FILE): + """ + build a route table using py-osrm + https://github.com/gis-ops/py-osrm + """ + engine = osrm.OSRM( + storage_config=database_path, + use_shared_memory=False + ) + n_demands = len(demand_sites) + n_supplys = len(candidate_depots) + query_params = osrm.TableParameters( # noqa: F821 + coordinates=[ + (float(lon), float(lat)) + for (lon, lat) + in numpy.row_stack((demand_sites, candidate_depots)) + ], + sources=list(numpy.arange(n_demands)), + destinations=list(numpy.arange(n_demands, n_demands + n_supplys)), + annotations=["distance"], + ) + res = engine.Table(query_params) + return numpy.asarray(res["distances"]).astype(float).T diff --git a/spopt/route/heuristic.py b/spopt/route/heuristic.py new file mode 100644 index 00000000..6be3e0b5 --- /dev/null +++ b/spopt/route/heuristic.py @@ -0,0 +1,464 @@ +import geopandas +import numpy +import pandas +import shapely + +import pyvrp +from . import engine +from . import utils + +class SpecificationError(Exception): + pass + +_MAX_INT = numpy.iinfo(numpy.int64).max + +class LastMile: + def __init__( + self, + depot_location=None, + depot_open=pandas.to_datetime("2030-01-02 07:45:00"), + depot_close=None, + depot_name=None, + cost_unit=1e-4, + ): + """ + Initialize a LastMile problem. + + Arguments + --------- + depot_location : tuple + longitude and latitude of depot for whom routes must be drawn. + depot_open : pandas.Datetime + The time from which trucks can start leaving the depot. + depot_close : pandas.Datetime + The time by which trucks must return to the depot. + depot_name : str + the name of the depot facility being sited + cost_unit : float + the fraction of cost values to round to. This is set to hundreths of + a cent by default (1e-4), assuming inputted costs are in typical + decimal notation for euros and cents. + """ + self.model = pyvrp.Model() + self.depot_location = depot_location + self.depot_open = depot_open + self.depot_close = depot_close + self.depot_name = depot_name + self.cost_unit = cost_unit + + def add_truck_type( + self, + name=None, + n_truck=None, + capacity=7, + time_windows=None, + fixed_cost=None, + cost_per_meter=None, + cost_per_minute=30 / 60, + max_duration=pandas.Timedelta(hours=8, minutes=00), + max_distance=None, + ): + """ + Add a single vehicle type to the LastMile problem. Must have + added clients first before setting up trucks, since the truck capacity + will be re-scaled to provide the count of minimum delivery values the + truck can carry. + + Parameters + ---------- + name : str + name of the truck type being added + n_truck : int + how many of this truck type are available to use + capacity : float + how big the truck is. Must be measured in the same units as + the client demand input + time_windows : pandas.DataFrame + data frame containing open and close times for this truck type's + routes. If not provided, these are set to the depot open and close + times. But, you may pass a custom value here if you want to, say, + force big noisy trucks to come back to the depot early + fixed_cost : float + the fixed cost per day of using a truck of this type + cost_per_meter : float + the variable cost per meter distance of using this truck + cost_per_minute : + the variable cost per minute time of using this truck + max_duration : pandas.Timedelta + the total allowed length of the route + max_distance : float + the total allowed distance (in meters) of the route. + + Returns + ------- + this LastMile object with these new truck types + """ + if not hasattr(self, "demand_unit_"): + raise SpecificationError( + "You must set the clients before adding truck types." + ) + if time_windows is None: + time_windows = pandas.Series([self.depot_open, self.depot_close], index=['open_1', 'close_1']).to_frame(self.depot_name).T + depot_ints = utils.integralize_time_windows(time_windows) + if name is None: + if hasattr(self, "trucks_"): + name = str(len(self.trucks_) + 1) + else: + name = "0" + v_ = self.model.add_vehicle_type( + name=str(name), + num_available=int(n_truck), + capacity=int(capacity / self.demand_unit_), + fixed_cost=int(fixed_cost / self.cost_unit), + tw_early=depot_ints[0,0], + tw_late=depot_ints[0,1], + max_duration=int(max_duration.total_seconds() / 60), + max_distance=max_distance if max_distance is not None else _MAX_INT, + unit_distance_cost=int(cost_per_meter / self.cost_unit), + unit_duration_cost=int(cost_per_minute / self.cost_unit), + ) + if hasattr(self, "trucks_"): + self.trucks_.append(v_) + else: + self.trucks_ = [v_] + return self + + def add_trucks_from_frame( + self, + truck_frame, + n_trucks=None + ): + """ + Add a single vehicle type to the LastMile problem. Must have + added clients first before setting up trucks, since the truck capacity + will be re-scaled to provide the count of minimum delivery values the + truck can carry. + + New truck types will be added to the self.trucks_ list. + + Parameters + ---------- + name : str + name of the truck type being added + n_trucks : int + how many trucks in total are allowed across all the inputted truck types. + If this option is provided, then the dataframe's n_truck column will govern + the *fraction* of trucks for each type, rounded down so that trucks do not + exceed the n_trucks limit. The fraction is calculated as + `truck_frame.n_truck / truck_frame.n_truck.sum()` + capacity : float + how big the truck is. Must be measured in the same units as + the client demand input + time_windows : pandas.DataFrame + data frame containing open and close times for this truck type's + routes. If not provided, these are set to the depot open and close + times. But, you may pass a custom value here if you want to, say, + force big noisy trucks to come back to the depot early + fixed_cost : float + the fixed cost per day of using a truck of this type + cost_per_meter : float + the variable cost per meter distance of using this truck + cost_per_minute : + the variable cost per minute time of using this truck + max_duration : pandas.Timedelta + the total allowed length of the route + max_distance : float + the total allowed distance (in meters) of the route. + + Returns + ------- + A LastMile() object modified in place to add a new truck type + """ + if n_trucks is not None: + raise NotImplementedError("will not yet calculate truck fractions") + keep_cols = truck_frame.columns.isin( + ['name', + 'n_truck', + 'capacity', + 'time_windows', + 'fixed_cost', + 'cost_per_meter', + 'cost_per_minute', + 'max_duration', + 'max_distance'] + ) + for name, row in truck_frame.loc[:,keep_cols].iterrows(): + truck_spec = row.to_dict() + if 'name' not in truck_spec.keys(): + truck_spec['name'] = name + self.add_truck_type(**truck_spec) + return self + + def add_clients( + self, + locations, + delivery=None, + pickup=None, + time_windows=None, + names=None, + service_times=None, + ): + """ + Add delivery targets to the model + + Parameters + ---------- + locations : geopandas.GeoSeries/geopandas.GeoDataFrame + the locations to add to the problem + delivery : numpy.ndarray + the demand values used to calculate load to be delivered to clients + pickup : numpy.ndarray + the demand values used to calculate load to be picked up from clients along a route + time_windows : pandas.DataFrame + open and close windows for each client. If every client has an open and a close, then + this dataframe should have two columns. If some clients have two open periods, + then this dataframe should have four columns, with open and close times interleaved. + names : numpy.ndarray + names to use for each client. this should serve as a unique ID for the client. + service_times : numpy.ndarray + the amount of time it takes to service the client, either as a string that is passed + directly to pandas.to_timedelta() or an existing iterable of timedelta objects. + + Returns + ------- + this LastMile() object with clients set to the clients_ attribute, + a depot objecti assigned as the depot_ attribute, and the unit of + demand assigned to the demand_unit_ attribute. + """ + coords = shapely.get_coordinates(locations.geometry) + all_coords = numpy.row_stack((self.depot_location, coords)) + n_clients = coords.shape[0] + + locints = utils.integralize_lonlat(all_coords) + depot_xy, client_xy = locints[0, :], locints[1:, :] + + self.depot_ = self.model.add_depot( + x=depot_xy[0], y=depot_xy[1], name=self.depot_name + ) + if time_windows is None: + time_windows = pandas.DataFrame.from_dict( + {"open_1":[self.depot_open]*len(locations.geometry), + "close_1": [self.depot_close]*len(locations.geometry) + }, + ) + time_windows.index = locations.geometry.index + if names is None: + names = locations.geometry.index + if service_times is None: + service_times = pandas.Series(numpy.zeros((n_clients,)), index=names).astype(int) + time_ints = utils.integralize_time_windows(time_windows) + if (delivery is None) | (pickup is None): + if (delivery is None): + delivery = pandas.Series(numpy.zeros((n_clients,)), index=names).astype(int) + if (pickup is None): + pickup = pandas.Series(numpy.zeros((n_clients,)), index=names).astype(int) + self.demand_unit_ = utils.calculate_minimum_unit(numpy.hstack((delivery, pickup))) + delivery_ints,_ = utils.integralize_demand(delivery, unit=self.demand_unit_) + pickup_ints,_ = utils.integralize_demand(pickup, unit=self.demand_unit_) + + service_timeblocks = numpy.ceil( + pandas.to_timedelta(service_times).dt.total_seconds().values / 60 + ) + clients = [] + + for i, name in enumerate(names): + x, y = client_xy[i] + x, y = x.item(), y.item() + delivery = delivery_ints.iloc[i].item() + pickup = pickup_ints.iloc[i].item() + service_duration = service_timeblocks[i].item() + tws = time_ints[i] + if len(tws) == 2: + lunchbreak=False + else: + lunchbreak = tws[1] > tws[2] + if not lunchbreak: # no lunchbreak + client_ = self.model.add_client( + x=x, + y=y, + delivery=delivery, + pickup=pickup, + tw_early=tws[0].item(), + tw_late=tws[1].item(), + name=str(name), + service_duration=service_duration, + ) + clients.append(client_) + else: + cg_ = self.model.add_client_group() + c1 = self.model.add_client( + x=x, + y=y, + delivery=delivery, + pickup=pickup, + tw_early=tws[0].item(), + tw_late=tws[1].item(), + name=name, + required=False, + service_duration=service_duration, + group=cg_, + ) + c2 = self.model.add_client( + x=x, + y=y, + delivery=delivery, + pickup=pickup, + tw_early=tws[2].item(), + tw_late=tws[3].item(), + name=name, + required=False, + service_duration=service_duration, + group=cg_, + ) + clients.extend([c1, c2]) + timedict = dict(open_1=time_windows.iloc[:, 0], close_1=time_windows.iloc[:, 1]) + if time_windows.shape[1] == 4: + timedict.update(dict(open_2=time_windows.iloc[:, 2], close_2=time_windows.iloc[:, 3])) + self.clients_ = geopandas.GeoDataFrame( + pandas.DataFrame.from_dict( + dict(delivery=delivery, pickup=pickup, service_time=service_times) + ).assign(**timedict), + index=names, + geometry=locations.geometry, + crs=locations.crs, + ) + return self + + def solve(self, stop=pyvrp.stop.NoImprovement(1e6), *args, **kwargs): + """ + Solve a LastMile() instance according to the existing specification. + + Parameters + ---------- + stop : pyvrp.stop.StoppingCriterion + A stopping rule that governs when the simulation will be ended. + Set to terminate solving after one million iterations with no improvement. + + + Returns + ------- + This LastMile() object, having added the results object to self.result_, as well + as the routes and stops found to routes_ and stops_, respectively + + Notes + ----- + other arguments and keyword arguments are passed directly to the pyvrp.Model.solve() method + """ + if (not hasattr(self, "clients_")) | (not hasattr(self, "trucks_")): + raise SpecificationError( + "must assign both clients and trucks to" " solve a problem instance." + ) + all_lonlats = numpy.row_stack( + (self.depot_location, shapely.get_coordinates(self.clients_.geometry)) + ) + self._setup_graph(all_lonlats=all_lonlats) + self.result_ = self.model.solve(stop=stop, *args, **kwargs) + self.routes_, self.stops_ = utils.routes_and_stops( + self.result_.best, self.model, self.clients_, self.depot_location, cost_unit=self.cost_unit + ) + return self + + solve.__doc__ = pyvrp.Model.solve.__doc__ + + def explore(self): + """ + Make a webmap of the solution, colored by the route name. + """ + if not hasattr(self, "routes_"): + raise SpecificationError("must have solved the model to show the result") + m = self.routes_.sort_values("route_name").explore( + "route_name", categorical=True, tiles="CartoDB positron" + ) + stops_for_map = self.stops_.copy() + stops_for_map["eta"] = self.stops_.eta.astype(str) + stops_for_map.explore( + "route_name", + m=m, + legend=False, + style_kwds=dict(color="black", radius=3, weight=1.5), + ) + geopandas.GeoDataFrame( + geometry=geopandas.points_from_xy( + x=[self.depot_location[0]], y=[self.depot_location[1]], crs="epsg:4326" + ) + ).explore(m=m, color="black", marker_type="marker") + return m + + def _setup_graph(self, all_lonlats): + """ + This sets up the graph pertaining to an inputted set of longitude and latitude coordinates. + + Note that this assumes that there is a single vehicle profile. + + TODO: For multiple vehicle profiles, we would need to identify + the restricted and the base profiles, then update the model + with an edge for each profile. + """ + raw_distances, raw_durations = engine.build_route_table( + all_lonlats, all_lonlats, cost="both" + ) + # how many minutes does it take to get from place to place? + durations_by_block = numpy.ceil(raw_durations / 60) + ##### WARNING!!!!!!! THIS IS A BUG IN OSRM #5855 + durations = numpy.clip(durations_by_block, 0, durations_by_block.max()) + distances = numpy.clip(raw_distances, 0, raw_distances.max()).round(0) + + duration_df = pandas.DataFrame( + durations, + index=[self.depot_name] + self.clients_.index.tolist(), + columns=[self.depot_name] + self.clients_.index.tolist(), + ) + distance_df = pandas.DataFrame( + distances, + index=[self.depot_name] + self.clients_.index.tolist(), + columns=[self.depot_name] + self.clients_.index.tolist(), + ) + for source_ix, source in enumerate(self.model.locations): + for sink_ix, sink in enumerate(self.model.locations): + self.model.add_edge( + source, + sink, + distance=distance_df.loc[source.name, sink.name].item(), + duration=duration_df.loc[ + source.name, sink.name + ].item(), # TODO: nogo zones + ) + + def write_result( + self, filestem=None, write_geometries=False + ): + """ + Write the result of a solution out to three files: + 1. the routes_ table is written to a file describing the route efficiency + 2. the stops_ table is written to a file describing each route + 3. the folium map from the .explore() method is written to a file + + This method requires that the model is solved first. + + Parameters + ---------- + filestem : str + start of the name for all output files. If not provided, then the + depot name is used as the default. + write_geometries : bool + whether or not to write the geometries out to a file. If True, the + output formats are geopackages. If False, as is default, then the + output formats are csvs. The folium map is always written to an html file. + + Returns + ------- + the operation writes to disk and returns None. + """ + if not hasattr(self, "routes_"): + raise SpecificationError("Model must be solved before results can be written.") + if write_geometries: + def writer(df, filename): + df.to_file(filename+".gpkg") + else: + def writer(df, filename): + df.drop("geometry", axis=1, errors='ignore').to_csv(filename + ".csv") + if filestem is None: + filestem = self.depot_name.replace(" ", "_") + writer(self.routes_, filestem+"_routes") + writer(self.stops_, filestem+"_stops") + self.explore().save(filestem+"_map.html") \ No newline at end of file diff --git a/spopt/route/utils.py b/spopt/route/utils.py new file mode 100644 index 00000000..79ed3e3f --- /dev/null +++ b/spopt/route/utils.py @@ -0,0 +1,514 @@ +import numpy +import pandas +import routing +import copy +import geopandas +import shapely +import randomname +MIDNIGHT = pandas.to_datetime( + "2030-01-02 00:00:00", format='%Y-%m-%d %H:%M:%S' +) + +def integralize_time_windows(windows, freq=pandas.Timedelta(minutes=1)): + """ + Convert time windows to a count of time units since midnight. + + Arguments + --------- + windows : pandas.DataFrame + a dataframe containing the time windows + freq : pandas.TimeDelta + the base time unit to use for counting. Defaults to one minute. + """ + timeblocks = pandas.Series(pandas.date_range( + start="2030-01-02 00:00:00", + end="2030-01-02 23:59:59", + freq=freq + )) + output = [] + for col in windows.columns: + dt = pandas.to_datetime(windows[col]) + itime = dt.apply( + lambda t: _integralize_time(t, timeblocks) + ) + output.append(itime) + return numpy.column_stack(output) + +def _integralize_time(t, timeblocks): + """ + Find the time block in timeblocks to which time t corresponds. + If time t is null, return -1, which is assumed to be handled elsewhere. + If time t is is not Null, then the time block into which t falls is + returned. + If time t does not fall into any block, then it's assumed to fall into + a block past the final block. + """ + if pandas.isnull(t): + return -1 + in_blocks = timeblocks.gt(t, fill_value=pandas.NaT) + if in_blocks.any(): + return in_blocks.argmax(skipna=False) + return len(timeblocks) + +def integralize_lonlat(lonlat, offset=True): + """ + Convert the lonlat values to integer meters relative to a UTM projection. + + Arguments + --------- + lonlat : numpy.ndarray or geopandas.GeoSeries/geopandas.GeoDataFrame + input geometries or coordinate values to convert to integer positions + offset : bool + whether or not to remove the false northing/easting on the UTM + reference frame, making the bottom left of the reference frame correspond to location 0,0. This is done by default. + + Returns + ------- + numpy.ndarray of integer coordinate locations to the nearest UTM-meter. + """ + if isinstance(lonlat, numpy.ndarray): + lon, lat = lonlat.T + lonlat = geopandas.GeoSeries(geopandas.points_from_xy( + x=lon, + y=lat, + crs='EPSG:4326' + ) + ) + elif isinstance(lonlat, (geopandas.GeoSeries,geopandas.GeoDataFrame)): + if lonlat.crs is None: + lonlat = lonlat.set_crs("epsg:4326") + utm = lonlat.estimate_utm_crs() + geoms = lonlat.to_crs(utm).geometry + raw_coords = numpy.column_stack((geoms.x, geoms.y)) + if not offset: + raw_coords -= geoms.total_bounds[:2] + return raw_coords.round(0).astype(int) + +def calculate_minimum_unit(vector): + """ + Calculate the smallest positive value in the vector. + Used throughout to estimate the integer "unit" for some of the + integralize_* functions. + """ + min_unit = vector[vector>0].min() + return min_unit + +def integralize_demand(demands, unit=None): + """ + Convert the demand values to integer units in terms of the minimum + nonzero demand. This converts the demand units to a count of + minimum nonzero demand values required to satisfy the client. Thus, + integralize(demands, unit) * unit will always be at least as big as + demands themselves, and may be larger by at most "unit" + + Arguments + --------- + demands : numpy.ndarray + demands used in the problem that must be integralized to a count. + unit : float + demand value to use as the unit count for the demands. + """ + if unit is None: + unit = calculate_minimum_unit(demands) + return numpy.ceil(demands/unit).astype(int), unit + +def routes_and_stops( + solution, + model, + target_geoms, + depot_location, + cost_unit=1e-4 + ): + """ + Calculate route geometries and stop etas/waypoint numbers from an input + solution. + + Arguments + --------- + solution : pyvrp.Solution + routing solution from running the pyvrp solver that describes the + best available routes to satisfy the constraints and specifications + recorded in the `model` + model : pyvrp.Model + the model reflecting the problem specification that has been solved + and recored in `solution` + target_geoms : geopandas.GeoSeries/geopandas.GeoDataFrame + the real-world longitude and latitude that correspond to the clients + recorded in the model. This should *not* include the depot location + which is provided as a separate argument, unless the depot is also + located at a client. + depot_location : tuple + the longitude and latitude values that correspond to the location + of the delivery depot + + Returns + ------- + two dataframes containing the routes and stops. the routes dataframe + will have one row per route, while the stops dataframe will be the same length + as the target_geoms input + """ + assert solution.is_feasible(), "solution is infeasible, routes cannot be used." + assert solution.is_complete(), "solution does not visit all required clients, routes cannot be used." + n_routes = solution.num_routes() + route_names = list(randomname.sample("adj/", "names/surnames/french", + n=n_routes + )) + + problem_data = model.data() + + # problem assumes all trucks have the same departure time + # problem assumes that this is in minutes since 00:00:00 + + route_lut = dict(zip(route_names, solution.routes())) + stops = [ + (route_name, r.visits()) + for route_name, r in route_lut.items() + ] + + stops = pandas.DataFrame( + stops + ).rename( + columns={0:"route_name", 1:"stop_idx"} + ).set_index("route_name") + + # calculate visit time, + # distances and durations are assumed constant over + # vehicle type + duration_matrix, = problem_data.duration_matrices() + distance_matrix, = problem_data.distance_matrices() + # TODO: would it be helpful to have the running capacity? + def timedelta_from_visits( + route, + duration_matrix=duration_matrix, + locations=model.locations + ): + """ + This is a private function to estimate the time changes + that evolve over a route using the model specific information, + rather than using the osrm-provided durations on demand. + This is to account for any waiting that occurs at the stops. + """ + full_visits = [0, *route.visits(), 0] + arrival_minutes = [route.start_time()] + for stop_number, stop_idx in enumerate(full_visits[:-1]): + next_stop_idx = full_visits[stop_number + 1] + travel_duration = duration_matrix[stop_idx, next_stop_idx] + # if service duration is not recorded, we assume + # there is no service time (like, for a waypoint) + service_duration = getattr( + locations[stop_idx], "service_duration", 0 + ) + # once you're at stop_idx, you spend service_duration + # there, and then spend travel_duration to get to the + # next spot. So, the deltas should be + # [0, service_duration[1] + travel_duration[0,1], ...] + # since the depot has service duration 0 + arrival_time = arrival_minutes[stop_number] + service_duration + travel_duration + # if you arrive at the target before it's open, then you have to wait + arrival_time = numpy.maximum( + getattr( + locations[stop_idx], + "tw_early", + -numpy.inf + ), arrival_time + ) + arrival_minutes.append(arrival_time) + tds = pandas.to_timedelta(arrival_minutes, unit='minutes') + return tds + + stops['eta'] = pandas.Series( + {name:timedelta_from_visits(r)[1:-1] + MIDNIGHT + for name,r in route_lut.items()} + ) + stops['stop_number'] = stops.stop_idx.apply(lambda x: numpy.arange(len(x))+1) + + big_stops = stops.explode( + ["stop_idx", "stop_number", "eta"] + ) + big_stops['target_uid'] = [ + model.locations[s].name for s in big_stops.stop_idx + ] + big_stops['stop_number'] = big_stops.groupby("route_name").cumcount().astype(int) + 1 + + stop_output = target_geoms.copy(deep=True) + stop_output = big_stops.reset_index().merge( + target_geoms, left_on='target_uid', right_index=True, + how='right' + ) + stop_output['route_name'] = stop_output.route_name.fillna("unassigned") + stop_output['stop_number'] = stop_output.stop_number.fillna(-1) + stop_output = stop_output.sort_values(["route_name","stop_number"]) + stop_output = geopandas.GeoDataFrame( + stop_output, + geometry='geometry', + crs=target_geoms.crs + ) + + route_data = [] + + for name, group in stop_output.groupby("route_name"): + route_obj = route_lut[name] + group = group.sort_values("stop_number") + coordinates = shapely.get_coordinates(group.geometry) + shape, durations = routing.build_specific_route( + numpy.vstack( + ( + depot_location, + coordinates, + depot_location + ) + ) + ) + route_truck_type = route_obj.vehicle_type() + truck_obj = model.vehicle_types[route_truck_type] + deptime, rettime = pandas.to_timedelta([ + route_obj.start_time(), + route_obj.end_time() + ], unit="minutes" + ) + MIDNIGHT + + route_data.append(( + name, + truck_obj.name, + route_obj.duration(), + route_obj.distance(), + route_obj.distance_cost() * cost_unit, + route_obj.duration_cost() * cost_unit, + truck_obj.fixed_cost * cost_unit, + ( route_obj.distance_cost() + + route_obj.duration_cost() + + truck_obj.fixed_cost + ) * cost_unit, + deptime, + rettime, + round(route_obj.duration() / truck_obj.max_duration * 100, 2), + round(route_obj.delivery() / truck_obj.capacity * 100, 2), + round(route_obj.distance() / truck_obj.max_distance * 100, 2), + shape + )) + + route_output = geopandas.GeoDataFrame( + pandas.DataFrame( + route_data, + columns = [ + 'route_name', + 'truck_type', + 'duration_min', + 'distance_m', + 'fuel_cost_€', + 'labor_cost_€', + 'truck_cost_€', + 'total_cost_€', + 'departure', + 'arrival', + 'utilization_time', + 'utilization_load', + 'utilization_rangelimit', + 'geometry' + ] + ), + geometry='geometry', + crs=target_geoms.crs + ) + + return route_output, stop_output + +def route_webmap( + problem, model, locs, depot_location, return_dataframes=False + ): + """ + Create a webmap from input data. + + Arguments + --------- + solution : pyvrp.Solution + routing solution from running the pyvrp solver that describes the + best available routes to satisfy the constraints and specifications + recorded in the `model` + model : pyvrp.Model + the model reflecting the problem specification that has been solved + and recored in `solution` + target_geoms : geopandas.GeoSeries/geopandas.GeoDataFrame + the real-world longitude and latitude that correspond to the clients + recorded in the model. This should *not* include the depot location + which is provided as a separate argument, unless the depot is also + located at a client. + depot_location : tuple + the longitude and latitude values that correspond to the location + of the delivery depot + return_dataframes : bool + whether or not to return the routes and stops dataframes as well + as the route map. If True, then the output is returned as a tuple + containing (map, routes, stops). Otherwise, the output is just the map. + + Returns + ------- + a folium.Map, and (if return_dataframes==True), two dataframes containing + the routes and stops. the routes dataframe will have one row per route, + while the stops dataframe will be the same length as the target_geoms input + """ + routes, stops = routes_and_stops( + problem, model, locs, depot_location + ) + m = routes.sort_values("route_name").explore( + "route_name", + categorical=True, + tiles="CartoDB positron" + ) + stops_for_map = stops.copy() + stops_for_map['eta'] = stops.eta.astype(str) + stops_for_map[ + ["target_uid", 'stop_number', 'route_name', 'eta', + '1rst_opening_hours', '1rst_closing_hours', + '2nd_opening_hours', '2nd_closing_hours', + 'sum_delivered_volume', 'geometry'] + ].explore("route_name", + m=m, legend=False, + style_kwds=dict(color='black', radius=3, weight=1.5) + ) + geopandas.GeoDataFrame( + geometry=geopandas.points_from_xy( + x=[depot_location[0]], y=[depot_location[1]], + crs="epsg:4326" + ) + ).explore(m=m, color="black", marker_type="marker") + if return_dataframes: + return m, routes, stops + return m + +def _to_dict(input, indices=None): + """ + construct a dictionary from an iterable for use in networkx updates + """ + if isinstance(input, dict): + return input + if indices is None: + indices = range(len(input)) + try: + return input.to_dict() + except AttributeError: + return dict(zip(indices, input)) + +def restrict_by_zone( + targets, + restrictions, + distances, + depot_to_points, + points_to_depot, + penalty_value=None + ): + """ + the restriction dataframe must be indexed by + the truck type that cannot go into that zone. + So, if truck type 2 cannot go into restriction + zone 7, the seventh row should have "2" as its index, + *not* "large petrol". + """ + distances = copy.deepcopy(distances) + depot_to_points = copy.deepcopy(depot_to_points) + points_to_depot = copy.deepcopy(points_to_depot) + try: + restrictions.index.astype(int) + last_index = restrictions.index.max() + distances[:,:,last_index] + except (ValueError, IndexError): + raise ValueError( + "the restrictions dataframe must be indexed by the integer" + " the restricted truck. So, if truck type 2 cannot go into" + " restriction zone 7, then the seventh row of `restrictions`" + " must have index 2, and distances[2] should be the cost/distance" + " matrix for that truck type." + ) + if penalty_value is None: + # works for any number of distance matrices + penalty_value = numpy.sum(distances) + + rzone_i, target_i = targets.sindex.query( + restrictions.geometry, + predicate='intersects' + ) + + for i, truck_i in enumerate(restrictions.index): + mask = (rzone_i == i) + if mask.sum() > 0: + # this is the entire row space; can't enter + distances[target_i[mask], :, truck_i] = penalty_value + # this is the entire column space; can't leave + distances[:, target_i[mask], truck_i] = penalty_value + # shouldn't start there + depot_to_points[target_i[mask], truck_i] = penalty_value + # or return from there + points_to_depot[target_i[mask], truck_i] = penalty_value + else: + continue + + return distances, depot_to_points, points_to_depot + + + + +def build_clients( + model, + xy_ints, + demands, + time_windows, + names, + service_durations, + ): + """ + Legacy function to construct a collection of clients + from an input set of demands, time windows, names, and service times. + """ + time_ints = integralize_time_windows(time_windows) + demand_ints, demand_units = integralize_demand( + demands + ) + service_timeblock_durations = numpy.ceil( + pandas.to_timedelta( + service_durations + ).dt.total_seconds().values / 60 + ) + clients = [] + for i,name in enumerate(names): + x, y = xy_ints[i] + x, y = x.item(), y.item() + delivery = demand_ints[i].item() + service_duration=service_timeblock_durations[i].item() + tws = time_ints[i] + if (tws[1]>tws[2]): # no lunchbreak + client_ = model.add_client( + x=x, + y=y, + delivery=delivery, + tw_early=tws[0].item(), + tw_late=tws[1].item(), + name=name, + service_duration = service_duration + ) + clients.append(client_) + else: + cg_ = model.add_client_group() + c1 = model.add_client( + x=x, + y=y, + delivery=delivery, + tw_early=tws[0].item(), + tw_late=tws[1].item(), + name=name, + required=False, + service_duration = service_duration, + group=cg_ + ) + c2 = model.add_client( + x=x, + y=y, + delivery=delivery, + tw_early=tws[2].item(), + tw_late=tws[3].item(), + name=name, + required=False, + service_duration = service_duration, + group=cg_ + ) + clients.extend([c1, c2]) + return clients \ No newline at end of file