cartodb4/crankshaft
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

[mvp] fractional optimization

Browse Source
This commit is contained in:
Andy Eschbacher
2017-05-08 23:47:17 -04:00
parent ddfdd87044
commit b0614aefc6
3 changed files with 137 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/pg/sql/25_optimization.sql
View File

@@ -1,11 +1,10 @@
CREATE OR REPLACE FUNCTION
CDB_OptimAssignments(drain text,
source text,
CREATE OR REPLACE FUNCTION
CDB_OptimAssignments(source text,
drain text,
drain_capacity text,
source_production text,
marginal_cost text,
waste_per_person numeric DEFAULT 0.01,
recycle_rate numeric DEFAULT 0.0,
dist_matrix_query text,
dist_rate numeric DEFAULT 0.15,
dist_threshold numeric DEFAULT null)
RETURNS table(drain_id bigint, source_id int, cost numeric) AS $$
@@ -15,13 +14,12 @@ from crankshaft.optimization import Optim
def cast_val(val):
return float(val) if val is not None else None
params = {'waste_per_person': cast_val(waste_per_person),
'recycle_rate': cast_val(recycle_rate),
'dist_rate': cast_val(dist_rate),
params = {'dist_rate': cast_val(dist_rate),
'dist_threshold': cast_val(dist_threshold)}
optim = Optim(drain, source, drain_capacity, source_production, marginal_cost,
**params)
optim = Optim(source, drain, dist_matrix_query, drain_capacity,
source_production, marginal_cost, **params)
x = optim.output()
return x

40
src/py/crankshaft/crankshaft/analysis_data_provider.py
View File

@@ -4,6 +4,8 @@ import pysal_utils as pu
import numpy as np
class AnalysisDataProvider(object):
"""Analysis providers for crankshaft functions. These rely on database
access through `plpy`"""
def get_getis(self, w_type, params):
"""fetch data for getis ord's g"""
try:
@@ -91,6 +93,44 @@ class AnalysisDataProvider(object):
resp = plpy.execute(query)
return np.array(resp[0]['col'], dtype=dtype)
def get_distance_matrix(self, query, origin_ids, destination_ids):
"""Transforms a SQL table origin-destination table into a distance
matrix.
:param query: Query that exposes the data needed for building the
distance matrix. Query should have the following columns:
- origin_id (int)
- destination_id (int)
- length_km (numeric)
:type query: str
:param origin_ids: List of origin IDs
:type origin_ids: list of ints
:param destination_ids: List of origin IDs
:type destination_ids: list of ints
:returns: 2D array of distances from all origins to all destinations
:rtype: numpy.array
"""
try:
resp = plpy.execute('''
SELECT "origin_id", "destination_id", "length_km"
FROM ({query}) as _wrap
'''.format(query=query))
except plpy.SPIError as err:
plpy.error("Failed to build distance matrix: {}".format(err))
pairs = {(row['origin_id'], row['destination_id']): row['length_km']
for row in resp}
distance_matrix = np.array([
pairs[(origin, destination)]
for destination in destination_ids
for origin in origin_ids
])
return np.array(distance_matrix,
dtype=float).reshape((len(destination_ids),
len(origin_ids)))
def get_pairwise_distances(self, drain_query, source_query,
id_col='cartodb_id'):
"""Retuns the pairwise distances between row i and j for all i in

160
src/py/crankshaft/crankshaft/optimization/optim.py
View File

@@ -9,47 +9,44 @@ class Optim(object):
"""Linear optimization class for logistics cost minimization
Optimization for logistics
based on models:
- amount_per_unit * (1 - recycle_rate) * population
- source_amount * (marginal_cost + transport_cost * distance)
That is, `cost ~ population * distance`
"""
def __init__(self, drain_query, source_table, capacity_column, # pylint: disable=too-many-arguments
production_column, marginal_column, **kwargs):
def __init__(self, source_query, drain_query, dist_matrix_query,
capacity_column, production_column, marginal_column,
**kwargs):
# set data provider - defaults to SQL database access
self.data_provider = kwargs.get('data_provider',
AnalysisDataProvider())
# model parameters
self.model_params = {
'amount_per_unit': kwargs.get('amount_per_unit', 0.01),
'dist_cost': kwargs.get('dist_cost', 0.15),
'recycle_rate': kwargs.get('recycle_rate', 0.0),
'dist_threshold': kwargs.get('dist_threshold', None)}
'dist_threshold': kwargs.get('dist_threshold', None),
'solver': kwargs.get('solver', 'glpk')}
self._check_model_params()
# model data
self.model_data = {
'drain_capacity': self.data_provider.get_column(drain_query,
capacity_column),
'source_amount': (self.model_params['amount_per_unit'] *
(1. - self.model_params['recycle_rate']) *
self.data_provider.get_column(source_table,
production_column)),
'marginal_cost': self.data_provider.get_column(drain_query,
marginal_column),
'distance': self.data_provider.get_pairwise_distances(source_table,
drain_query)
}
self.model_data['cost'] = self.calc_cost()
# database ids
self.ids = {
'drain': self.data_provider.get_column(drain_query,
'cartodb_id',
dtype=int),
'source': self.data_provider.get_column(source_table,
'source': self.data_provider.get_column(source_query,
'cartodb_id',
dtype=int)}
# model data
self.model_data = {
'drain_capacity': self.data_provider.get_column(drain_query,
capacity_column),
'source_amount': self.data_provider.get_column(source_query,
production_column),
'marginal_cost': self.data_provider.get_column(drain_query,
marginal_column),
'distance':
self.data_provider.get_distance_matrix(dist_matrix_query,
self.ids['source'],
self.ids['drain'])}
self.model_data['cost'] = self.calc_cost()
self.n_sources = len(self.ids['source'])
self.n_drains = len(self.ids['drain'])
@@ -73,22 +70,17 @@ class Optim(object):
def _check_model_params(self):
"""Ensure model parameters are well formed"""
if (self.model_params['recycle_rate'] is None or
self.model_params['recycle_rate'] < 0 or
self.model_params['recycle_rate'] > 1):
raise ValueError("`recycle_rate` must be between 0 and 1.")
if (self.model_params['amount_per_unit'] is None or
self.model_params['amount_per_unit'] < 0):
raise ValueError("`amount_per_unit` must be greater than zero.")
if (self.model_params['dist_threshold'] <= 0 and
self.model_params['dist_threshold'] is not None):
raise ValueError("`dist_threshold` must be greater than zero")
if (self.model_params['dist_cost'] is None or
self.model_params['dist_cost'] < 0):
raise ValueError("`dist_cost must be greater than zero")
raise ValueError("`dist_cost` must be greater than zero")
if self.model_params['solver'] not in (None, 'glpk'):
raise ValueError("`solver` must be one of 'glpk' (default) "
"or None.")
return None
@@ -99,8 +91,7 @@ class Optim(object):
transport from source to drain
:rtype: List of tuples
"""
# n_drains x n_sources matrix (row, column)
# retrieve fractional assignments
assignments = self.optim()
# crosswalks for matrix index -> cartodb_id
@@ -118,11 +109,11 @@ class Optim(object):
nonzeros = np.nonzero(assignments)
source_index, drain_index = nonzeros[0], nonzeros[1]
#
assigned_costs = [(drain_id_crosswalk[drain_index[i]],
source_id_crosswalk[source_index[i]],
self.model_data['cost'][drain_index[i],
source_index[i]])
for i in range(len(source_index))]
assigned_costs = [(drain_id_crosswalk[drain_index[idx]],
source_id_crosswalk[source_val],
self.model_data['cost'][drain_index[idx],
source_val])
for idx, source_val in enumerate(source_index)]
return assigned_costs
def cost_func(self, distance, waste, marginal):
@@ -148,12 +139,13 @@ class Optim(object):
Populate an d x s matrix according to the cost equation
:returns: d x s matrix of costs from area i to plant j
:rtype: NumPy matrix
:rtype: numpy.array
"""
costs = np.array([self.cost_func(distance,
self.model_data['source_amount'][pair[1]],
self.model_data['marginal_cost'][pair[0]])
for pair, distance in np.ndenumerate(self.model_data['distance'])])
costs = np.array(
[self.cost_func(distance,
self.model_data['source_amount'][pair[1]],
self.model_data['marginal_cost'][pair[0]])
for pair, distance in np.ndenumerate(self.model_data['distance'])])
return costs.reshape(self.model_data['distance'].shape)
def optim(self):
@@ -163,51 +155,77 @@ class Optim(object):
minimize c'*x by assigning x values
subject to G*x <= h
A*x = b
x[k] is binary
:returns: Assignments array (of 1s and 0s) of shape c.T
:rtype: NumPy array
0 <= x[k] <= 1
:returns: Fractional assignments array (of 1s and 0s) of shape c.T.
Value at position (i, j) corresponds to the fraction of source
`i`'s supply to drain `j`.
:rtype: numpy.array
"""
n_pairings = self.n_sources * self.n_drains
# ---
# costs
# elements chosen to minimize sum
cost = cvxopt.matrix(self.model_data['cost'].ravel('F'))
cost = np.nan_to_num(self.model_data['cost'])
cost = cvxopt.matrix(cost.ravel('F'))
# ---
# equality constraint variables
# each area is serviced once
A = cvxopt.spmatrix(1., # pylint: disable=invalid-name
A = cvxopt.spmatrix(1.,
[i // self.n_drains
for i in range(self.n_drains * self.n_sources)],
range(self.n_drains * self.n_sources))
b = cvxopt.matrix(np.ones((self.n_sources, 1)), tc='d') # pylint: disable=invalid-name
for i in range(n_pairings)],
range(n_pairings), tc='d')
b = cvxopt.matrix([1.] * self.n_sources, tc='d')
# make nan's in cost impossible
if np.isnan(self.model_data['distance']).any():
i_vals, j_vals = np.where(np.isnan(self.model_data['distance']))
for idx, i_val in enumerate(i_vals):
i = int(i_val)
j = int(i_val * self.n_drains + j_vals[idx])
A[i, j] = 0
# knock out values above distance threshold
if self.model_params['dist_threshold']:
j_locs, i_locs = np.where(self.model_data['distance'] > 100)
for idx, ival in enumerate(i_locs):
A[int(ival), int(ival * 10 + j_locs[idx])] = 0
j_vals, i_vals = np.where(self.model_data['distance'] >
self.model_params['dist_threshold'])
for idx, ival in enumerate(i_vals):
A[int(ival), int(ival * self.n_drains + j_vals[idx])] = 0
# ---
# inequality constraint variables
# each plant never goes over capacity
drain_capacity = cvxopt.matrix(self.model_data['drain_capacity'],
tc='d')
source_amounts = cvxopt.spmatrix(
np.repeat(self.model_data['source_amount'], self.n_drains),
[i % self.n_drains for i in range(self.n_drains * self.n_sources)],
range(self.n_drains * self.n_sources))
drain_capacity = cvxopt.matrix([
cvxopt.matrix(self.model_data['drain_capacity'], tc='d'),
cvxopt.matrix([1.] * n_pairings, tc='d'),
cvxopt.matrix([0.] * n_pairings, tc='d')
])
binary_entries = set(range(self.n_drains * self.n_sources))
# inequality maxima
ineq_maxs = cvxopt.sparse([
cvxopt.spmatrix(
np.repeat(self.model_data['source_amount'], self.n_drains),
[i % self.n_drains for i in range(n_pairings)],
range(n_pairings), tc='d'),
cvxopt.spmatrix(1.,
range(n_pairings),
range(n_pairings)),
cvxopt.spmatrix(-1.,
range(n_pairings),
range(n_pairings))
], tc='d')
# solve
(sol, assignments) = ilp(c=cost, G=source_amounts, h=drain_capacity,
A=A, b=b, B=binary_entries)
if sol != 'optimal':
sol = cvxopt.solvers.lp(c=cost, G=ineq_maxs, h=drain_capacity,
A=A, b=b, solver=self.model_params['solver'])
if sol['status'] != 'optimal':
raise Exception("No solution possible: {}".format(sol))
assign_shape = (self.model_data['cost'].shape[1],
self.model_data['cost'].shape[0])
# Note: assignments needs to be shaped like self.model_data['cost'].T
return np.array(assignments,
dtype=int).flatten().reshape(assign_shape)
# NOTE: assignments needs to be shaped like self.model_data['cost'].T
return np.array(sol['x'],
dtype=float)\
.flatten()\
.reshape((self.model_data['cost'].shape[1],
self.model_data['cost'].shape[0]))