From 78f72d25e6524e01a9dfa5270c781a047bdaa30a Mon Sep 17 00:00:00 2001 From: Andy Eschbacher Date: Thu, 30 Mar 2017 15:11:12 -0400 Subject: [PATCH] refactoring --- src/pg/sql/25_optimization.sql | 10 +-- .../crankshaft/optimization/optim.py | 72 ++++++++++++++----- 2 files changed, 60 insertions(+), 22 deletions(-) diff --git a/src/pg/sql/25_optimization.sql b/src/pg/sql/25_optimization.sql index 5636a19..be1d065 100644 --- a/src/pg/sql/25_optimization.sql +++ b/src/pg/sql/25_optimization.sql @@ -6,14 +6,16 @@ CDB_OptimAssignments(drain text, marginal_cost text, waste_per_person numeric DEFAULT 0.01, recycle_rate numeric DEFAULT 0.0, - dist_rate numeric DEFAULT 0.15) + dist_rate numeric DEFAULT 0.15, + dist_threshold numeric DEFAULT null) RETURNS table(drain_id bigint, source_id int, cost numeric) AS $$ from crankshaft.optimization import Optim -params = {'waste_per_person': float(waste_per_person), - 'recycle_rate': float(recycle_rate), - 'dist_rate': float(dist_rate)} +params = {'waste_per_person': waste_per_person, + 'recycle_rate': recycle_rate, + 'dist_rate': dist_rate, + 'dist_threshold': dist_threshold} optim = Optim(drain, source, drain_capacity, source_production, marginal_cost, **params) diff --git a/src/py/crankshaft/crankshaft/optimization/optim.py b/src/py/crankshaft/crankshaft/optimization/optim.py index d446b8e..5bacf6b 100644 --- a/src/py/crankshaft/crankshaft/optimization/optim.py +++ b/src/py/crankshaft/crankshaft/optimization/optim.py @@ -24,7 +24,9 @@ class Optim(object): 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)} + 'recycle_rate': kwargs.get('recycle_rate', 0.0), + 'dist_threshold': kwargs.get('dist_threshold', None)} + self._check_model_params() # model data self.model_data = { @@ -35,7 +37,11 @@ class Optim(object): self.data_provider.get_column(source_table, production_column)), 'marginal_cost': self.data_provider.get_column(drain_table, - marginal_column)} + marginal_column), + 'distance': self.data_provider.get_pairwise_distances(source_table, + drain_table), + 'cost': self.calc_cost() + } # database ids self.ids = { @@ -45,20 +51,46 @@ class Optim(object): 'source': self.data_provider.get_column(source_table, 'cartodb_id', dtype=int)} - # derivative data self.n_sources = len(self.ids['source']) self.n_drains = len(self.ids['drain']) - self.cost = self.calc_cost(source_table, - drain_table) def _check_constraints(self): """Check if inputs are within constraints""" - if self.model_data['source_amount'].sum() > self.model_data['drain_capacity'].sum(): + if (self.model_data['source_amount'].sum() > + self.model_data['drain_capacity'].sum()): plpy.error("Solution not possible. Drain capacity is smaller " "than total source production.") + return None + + 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'] is None or + self.model_params['dist_threshold'] < 0): + 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") + + return None def output(self): - """...""" + """Output the calculated 'optimal' assignments if solution is not infeasible. + + :returns: List of source id/drain id pairs and the associated cost of + transport from source to drain + :rtype: List of tuples + """ # n_drains x n_sources matrix (row, column) assignments = self.optim() @@ -80,8 +112,8 @@ class Optim(object): # assigned_costs = [(drain_id_crosswalk[drain_index[i]], source_id_crosswalk[source_index[i]], - self.cost[drain_index[i], - source_index[i]]) + self.model_data['cost'][drain_index[i], + source_index[i]]) for i in range(len(source_index))] return assigned_costs @@ -103,26 +135,24 @@ class Optim(object): """ return waste * (marginal + self.model_params['dist_cost'] * distance) - def calc_cost(self, source_table, drain_table): + def calc_cost(self): """ 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 """ - distances = self.data_provider.get_pairwise_distances(source_table, - drain_table) 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(distances)]) - return costs.reshape(distances.shape) + for pair, distance in np.ndenumerate(self.model_data['distance'])]) + return costs.reshape(self.model_data['distance'].shape) def optim(self): """solve linear optimization problem Equations of the form: - minimize c'*x + minimize c'*x by assigning x values subject to G*x <= h A*x = b x[k] is binary @@ -132,7 +162,7 @@ class Optim(object): # --- # costs # elements chosen to minimize sum - cost = cvxopt.matrix(self.cost.ravel('F')) + cost = cvxopt.matrix(self.model_data['cost'].ravel('F')) # --- # equality constraint variables @@ -143,6 +173,12 @@ class Optim(object): range(self.n_drains * self.n_sources)) b = cvxopt.matrix(np.ones((self.n_sources, 1)), tc='d') + # 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 + # --- # inequality constraint variables # each plant never goes over capacity @@ -160,8 +196,8 @@ class Optim(object): A=A, b=b, B=binary_entries) if sol != 'optimal': raise Exception("No solution possible: {}".format(sol)) - assign_shape = (self.cost.shape[1], self.cost.shape[0]) + assign_shape = (self.model_data['cost'].shape[1], self.model_data['cost'].shape[0]) - # Note: assignments needs to be shaped like self.cost.T + # Note: assignments needs to be shaped like self.model_data['cost'].T return np.array(assignments, dtype=int).flatten().reshape(assign_shape)