diff --git a/src/pg/sql/25_optimization.sql b/src/pg/sql/25_optimization.sql index 79a2110..79e3b4d 100644 --- a/src/pg/sql/25_optimization.sql +++ b/src/pg/sql/25_optimization.sql @@ -4,12 +4,12 @@ CDB_OptimAssignments(drain text, drain_capacity text, source_production text, marginal_cost text) -RETURNS setof int AS $$ +RETURNS table(drain_id bigint, source_id int, cost numeric) AS $$ from crankshaft.optimization import Optim optim = Optim(drain, source, drain_capacity, source_production, marginal_cost) -x = optim.optim() +x = optim.output() print(x) return x diff --git a/src/py/crankshaft/crankshaft/analysis_data_provider.py b/src/py/crankshaft/crankshaft/analysis_data_provider.py index 22ac1c6..e2da80e 100644 --- a/src/py/crankshaft/crankshaft/analysis_data_provider.py +++ b/src/py/crankshaft/crankshaft/analysis_data_provider.py @@ -66,15 +66,25 @@ class AnalysisDataProvider(object): except plpy.SPIError, err: plpy.error('Analysis failed: %s' % err) - def get_column(self, table, column): + def get_column(self, table, column, dtype=float): """ + Retrieve the column from the specified table + + :param table: table to retrieve column from + :type table: text + :param column: column to retrieve + :type column: text + :param dtype: data type in column (e.g, float, int, str) + :type dtype: type + :returns: column from table as a NumPy array + :rtype: NumPy array """ query = ''' SELECT array_agg("{column}" ORDER BY "cartodb_id" ASC) as col FROM "{table}" '''.format(table=table, column=column) resp = plpy.execute(query) - return np.array(resp[0]['col'], dtype=float) + return np.array(resp[0]['col'], dtype=dtype) def get_pairwise_distances(self, drain, source): """retuns the pairwise distances between row i and j for all i in table1 and j in table1""" diff --git a/src/py/crankshaft/crankshaft/optimization/optim.py b/src/py/crankshaft/crankshaft/optimization/optim.py index ae9f94b..01124a9 100644 --- a/src/py/crankshaft/crankshaft/optimization/optim.py +++ b/src/py/crankshaft/crankshaft/optimization/optim.py @@ -1,9 +1,8 @@ """optimization""" -import sys +import plpy +import numpy as np import cvxopt from cvxopt.glpk import ilp -import numpy as np -import plpy from crankshaft.analysis_data_provider import AnalysisDataProvider class Optim(object): @@ -33,6 +32,13 @@ class Optim(object): production_column)) self.marginal_cost = self.data_provider.get_column(drain_table, marginal_column) + # database ids + self.drain_ids = self.data_provider.get_column(drain_table, + 'cartodb_id', + dtype=int) + self.source_ids = self.data_provider.get_column(source_table, + 'cartodb_id', + dtype=int) # derivative data self.distances = self.data_provider.get_pairwise_distances(source_table, drain_table) @@ -40,6 +46,45 @@ class Optim(object): self.n_plants = len(self.distances) self.cost = self.calc_cost() + def output(self): + """...""" + + # n_drains x n_sources matrix (row, column) + assignments = self.optim() + # + plpy.notice("self.cost.shape: {}".format(str(self.cost.shape))) + plpy.notice("self.cost: {}".format(self.cost)) + # + plpy.notice(assignments) + plpy.notice("assignments.shape: {}".format(str(assignments.shape))) + + # crosswalks for matrix index -> cartodb_id + drain_id_crosswalk = {} + for idx, cid in enumerate(self.drain_ids): + # matrix index -> cartodb_id + drain_id_crosswalk[idx] = cid + # plpy.notice(drain_id_crosswalk) + + source_id_crosswalk = {} + for idx, cid in enumerate(self.source_ids): + # matrix index -> cartodb_id + source_id_crosswalk[idx] = cid + # plpy.notice(source_id_crosswalk) + + # find non-zero entries + nonzeros = np.nonzero(assignments) + plpy.notice("nonzeros: {}".format(str(nonzeros))) + source_index, drain_index = nonzeros[0], nonzeros[1] + # + plpy.notice(source_index) + plpy.notice(drain_index) + assigned_costs = [(drain_id_crosswalk[drain_index[i]], + source_id_crosswalk[source_index[i]], + self.cost[drain_index[i], + source_index[i]]) + for i in range(len(source_index))] + return assigned_costs + def test(self): """ just plpy.notice the stored information @@ -75,7 +120,7 @@ class Optim(object): self.waste_in_area[pair[1]], self.marginal_cost[pair[0]]) for pair, distance in np.ndenumerate(self.distances)]) - return costs + return costs.reshape(self.distances.shape) def optim(self): """solve linear optimization problem @@ -88,27 +133,30 @@ class Optim(object): """ # costs # elements chosen to minimize sum + # NOTE: used to be ravel('F') c = cvxopt.matrix(self.cost.ravel('F')) # equality constraint variables # each area is serviced once - A = cvxopt.spmatrix(1., [i // self.n_plants - for i in range(self.n_plants * self.n_areas)], - range(self.n_plants * self.n_areas)) + A = cvxopt.spmatrix(1., + [i // self.n_plants + for i in range(self.n_plants * self.n_areas)], + range(self.n_plants * self.n_areas)) b = cvxopt.matrix(np.ones((self.n_areas, 1)), tc='d') # inequality constraint variables # each plant never goes over capacity - h = cvxopt.matrix(self.plant_capacity) + h = cvxopt.matrix(self.plant_capacity, tc='d') G = cvxopt.spmatrix(np.repeat(self.waste_in_area, self.n_plants), [i % self.n_plants for i in range(self.n_plants * self.n_areas)], range(self.n_plants * self.n_areas)) - + binary_entries = set(range(len(c))) # solve - sol, x = ilp(c=c, G=G, h=h, A=A, b=b) + (sol, x) = ilp(c=c, G=G, h=h, A=A, b=b, B=binary_entries) # assignment = np.array(x).reshape((self.n_areas, self.n_plants)) if sol != 'optimal': - raise Exception("Solution not soluble: {}".format(sol)) - - return np.array(x) + raise Exception("No solution possible: {}".format(sol)) + x_shape = (self.cost.shape[1], self.cost.shape[0]) + # Note: x needs to be shaped like self.cost.T + return np.array(x, dtype=int).flatten().reshape(x_shape)