import numpy as np
import pandas as pd
def iterate_raam(demand, supply, travel,
max_cycles = 151, initial_step = 0.2, min_step = 0.005, half_life = 50,
limit_initial = 20, verbose = False):
norig, ndest = travel.shape
assignment = np.zeros((norig, ndest))
assignment[range(norig), travel.argmin(axis = 1)] = demand
for i in range(max_cycles):
demand_at_supply = assignment.sum(axis = 0)
congestion_cost = demand_at_supply / supply
total_cost = (congestion_cost + travel)
max_locations = np.ma.masked_array(total_cost, assignment == 0).argmax(axis = 1)
min_locations = total_cost.argmin(axis = 1)
slmin = supply[min_locations]
slmax = supply[max_locations]
trlmin = travel[range(norig), min_locations]
trlmax = travel[range(norig), max_locations]
drlmin = assignment[range(norig), min_locations]
drlmax = assignment[range(norig), max_locations]
dr = drlmin + drlmax
drotherlmin = demand_at_supply[min_locations] - drlmin
drotherlmax = demand_at_supply[max_locations] - drlmax
drlmin_new = ((slmin * slmax) / (slmin + slmax)) * \
((trlmax - trlmin) + (dr + drotherlmax) / slmax - drotherlmin / slmin)
delta = drlmin_new - drlmin
delta = np.minimum(delta, drlmax)
delta = np.where(max_locations == min_locations, 0, delta)
if type(initial_step) is float:
step_size = initial_step * 0.5 ** (i / half_life)
if step_size < min_step: step_size = min_step
delta = np.minimum(delta, step_size * demand).astype(int)
else:
step_size = int(np.round(initial_step * 0.5 ** (i / half_life)))
if step_size < min_step: step_size = min_step
delta = np.minimum(delta, step_size).astype(int)
## We don't want "attractive locations" getting mobbed.
## This will only happen in the first 10-20 cycles.
## So only do these (somewhat costly checks) then.
if i < limit_initial:
delta_mat = np.zeros(travel.shape)
delta_mat[range(norig), min_locations] += delta
naive_assignment = delta_mat.sum(axis = 0) / (supply)# * rho)
scale_factor = np.maximum(naive_assignment, 1)
delta_mat = (delta_mat / scale_factor).round().astype(int)
delta = delta_mat.sum(axis = 1)
assignment[range(norig), min_locations] += delta
assignment[range(norig), max_locations] -= delta
assert((assignment.sum(axis = 1) == demand).all())
if not (i % 25):
raam_cost = (total_cost * assignment).sum(axis = 1) / assignment.sum(axis = 1)
if verbose:
print("{:d} {:.2f} {:d} {:.3f}".format(i, raam_cost.mean(), delta.sum(), step_size), end = " || ")
raam_cost = (total_cost * assignment).sum(axis = 1) / assignment.sum(axis = 1)
return raam_cost
[docs]def raam(demand_df, supply_df, cost_df,
demand_index = True, demand_name = "demand",
supply_index = True, supply_name = "supply",
cost_origin = "origin", cost_dest = "dest", cost_name = "cost",
tau = 60, rho = None,
max_cycles = 150, initial_step = 0.2, min_step = 0.005, half_life = 50,
verbose = False):
"""Calculate the rational agent access model's total cost --
a weighted travel and congestion cost.
The balance of the two costs is expressed by the
:math:`\\tau` parameter, which corresponds to the travel time
required to accept of congestion by 100% of the mean demand to supply ratio
in the study area.
Parameters
----------
demand_df : `pandas.DataFrame <https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html>`_
The origins dataframe, containing a location index and a total demand.
demand_origin : str
is the name of the column of `demand` that holds the origin ID.
demand_value : str
is the name of the column of `demand` that holds the aggregate demand at a location.
supply_origin : str
is the name of the column of `demand` that holds the origin ID.
supply_df : `pandas.DataFrame <https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html>`_
The origins dataframe, containing a location index and level of supply
cost_df : `pandas.DataFrame <https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html>`_
This dataframe contains a link between neighboring demand locations, and a cost between them.
cost_origin : str
The column name of the locations of users or consumers.
cost_dest : str
The column name of the supply or resource locations.
cost_name : str
The column name of the travel cost between origins and destinations
weight_fn : function
This fucntion will weight the value of resources/facilities,
as a function of the raw cost.
max_cycles : int
Max number of cycles.
max_shift : int
This is the maximum number to shift in each cycle.
max_cost : float
This is the maximum cost to consider in the weighted sum;
note that it applies along with the weight function.
Returns
-------
access : pandas.Series
A -- potentially-weighted -- Rational Agent Access Model cost.
"""
if demand_index is not True: demand_df = demand_df.set_index(demand_index)
if supply_index is not True: supply_df = supply_df.set_index(supply_index)
demand_df = demand_df[demand_df[demand_name] > 0].copy()
supply_df = supply_df[supply_df[supply_name] > 0].copy()
demand_locations = list(set(cost_df[cost_origin]) & set(demand_df.index))
supply_locations = list(set(cost_df[cost_dest]) & set(supply_df.index))
cost_pivot = cost_df.pivot(index=cost_origin, columns=cost_dest, values=cost_name)
try:
travel_np = cost_pivot.loc[demand_locations, supply_locations].to_numpy().copy()
except:
travel_np = cost_pivot.loc[demand_locations, supply_locations].values.copy()
travel_np = travel_np / tau
travel_np = np.ma.masked_array(travel_np, np.isnan(travel_np))
# If it is not specified, rho is the average demand to supply ratio.
if rho is None: rho = demand_df[demand_name].sum() / supply_df[supply_name].sum()
try:
supply_np = supply_df.loc[supply_locations, supply_name].to_numpy().copy()
except:
supply_np = supply_df.loc[supply_locations, supply_name].values.copy()
supply_np = supply_np * rho
# Change this -- should be
try:
demand_np = demand_df.loc[demand_locations, demand_name].to_numpy().copy()
except:
demand_np = demand_df.loc[demand_locations, demand_name].values.copy()
raam_cost = iterate_raam(demand_np, supply_np, travel_np, verbose = verbose,
max_cycles = max_cycles, initial_step = initial_step,
min_step = min_step, half_life = half_life)
rs = pd.Series(name = "RAAM", index = demand_locations, data = raam_cost)
return rs