greedy.py

import numpy as np
import scipy as sp
from numpy import random
from numpy import matlib
import matplotlib.pyplot as plt
from gurobipy import GRB
import gurobipy as grb


class Problem:
    def __init__(self, C=20, F=15):
        self.C = C
        self.F = F
        
        self.clients = np.random.rand(2, C)  # client positions
        self.facilities = np.random.rand(2, F)  # facility positions
        
        # maximum number of clients per facility
        self.capacities = np.ones((F,), dtype=np.int32) * 4
        
        # assignment cost is defined as the squared distance between a client and a facility
        # print(np.matlib.repmat(self.clients[0, :], F, 1) )
        dx = \
            np.matlib.repmat(self.clients[0, :], F, 1) - \
            np.matlib.repmat(self.facilities[0, :], C, 1).transpose()
        # print("dx\n",dx)
        dy = \
            np.matlib.repmat(self.clients[1, :], F, 1) - \
            np.matlib.repmat(self.facilities[1, :], C, 1).transpose()
        # print("dy\n", dy)
        self.assignment_costs = 3 * (dx * dx + dy * dy)
        print(self.assignment_costs)
        print("=====")
        # self.assignment_costs.sort()
        
        self.opening_costs = np.ones((F,))
    
    def assign_random_capacities(self):
        """
		Assign more or less random capacities to facilities.

		This is one of the possible ways to change the problem configuration.
		In other words, use this function when testing your solution!
		"""
        while True:
            self.capacities = \
                np.random.randint(2 * self.C // self.F, size=self.F) + 1
            if sum(self.capacities) > self.C * 1.3:
                break
    
    def assign_random_opening_costs(self):
        """
		Assign more or less random opening costs to facilities.

		Same as above -- use this for your report.
		"""
        # could be float, but let it be simple
        self.opening_costs = \
            np.random.randint((self.C + self.F - 1) // self.F, size=self.F) + 1
    
    def plot(self, y, assignments, fig=plt):
        """
		Plot the given solution (y, assignments)

		Arguments:
		y, assignments -- see Problem.objective().
		fig            -- an instance of matplotlib.axes._axes.Axes to draw on.
						  Also, can be matplotlib.pyplot, in this case use the default Axes.
						  This is useful to compare your results (see "Results comparison" section).
		"""
        
        # y = y.astype(np.int32)
        # assignments = assignments.astype(np.int32)
        
        y = np.array(y, dtype=np.int32)
        assignments = np.array(assignments, np.int32)
        
        for cli, fac in enumerate(assignments):
            fig.plot([self.clients[0, cli], self.facilities[0, fac]], \
                     [self.clients[1, cli], self.facilities[1, fac]], c=(.7, .7, .7))
        
        fig.scatter(self.clients[0, :], self.clients[1, :], s=15.0, c=assignments, \
                    vmin=0, vmax=self.F - 1)
        
        fig.scatter(self.facilities[0, :], self.facilities[1, :], s=54.0, \
                    c=range(self.F), linewidth=[1 * el for el in y])
    
    def objective(self, y, assignments):
        """
		Return objective function value given a solution.
		If the solution is infeasible, return infinity.

		Arguments:
		y           -- a binary 1D array of size F. y[i] is 1 iff i-th facility is open.
		assignments -- an integer 1D array of size C. assignments[i] is index of facility
					   that i-th client is assigned to.
		"""
        assert len(y) == self.F
        assert len(assignments) == self.C
        
        # y = y.astype(np.int32)
        # assignments = assignments.astype(np.int32)
        
        y = np.array(y, dtype=np.int32)
        assignments = np.array(assignments, np.int32)
        
        retval = sum(is_opened * opening_cost \
                     for is_opened, opening_cost in zip(y, self.opening_costs))
        
        assignment_counts = np.zeros_like(y)
        
        for cli, fac in enumerate(assignments):
            if not y[fac]:
                return np.inf
            else:
                retval += self.assignment_costs[fac, cli]
                assignment_counts[fac] += 1
        
        if any(assignment_counts > self.capacities):
            return np.inf
        
        return retval
    
    def solve_gurobi(self, verbose=False):
        """
		Solve the problem using mixed integer program solver.
		Return `y, assignments` (see Problem.objective() docstring for format).

		Arguments:
		verbose -- controls Gurobi output.
		"""
        m = grb.Model("facility")
        
        y = []
        for i_f in range(self.F):
            y.append(m.addVar(vtype=GRB.BINARY))
        
        x = []
        for i_f in range(self.F):
            x.append([])
            for i_c in range(self.C):
                x[i_f].append(m.addVar(vtype=GRB.BINARY))
        
        # the objective is to minimize the total fixed and variable costs
        m.modelSense = GRB.MINIMIZE
        
        # update model to integrate new variables
        m.update()
        
        # set optimization objective - minimize sum of fixed costs
        obj_summands = []
        for i_f in range(self.F):
            obj_summands.append(self.opening_costs[i_f] * y[i_f])
        
        for i_f in range(self.F):
            for i_c in range(self.C):
                obj_summands.append(self.assignment_costs[i_f][i_c] * x[i_f][i_c])
        
        m.setObjective(grb.quicksum(obj_summands))
        
        # set constraints
        for i_c in range(self.C):
            client_constr_summands = [x[i_f][i_c] for i_f in range(self.F)]
            m.addConstr(sum(client_constr_summands), GRB.EQUAL, 1.0)
        
        for i_f in range(self.F):
            facility_constr_summands = [x[i_f][i_c] for i_c in range(self.C)]
            m.addConstr(sum(facility_constr_summands), \
                        GRB.LESS_EQUAL, self.capacities[i_f] * y[i_f])
        
        for i_f in range(self.F):
            for i_c in range(self.C):
                m.addConstr(x[i_f][i_c], GRB.LESS_EQUAL, y[i_f])
        
        # optimize
        m.setParam(GRB.Param.OutputFlag, verbose)
        m.optimize()
        
        facilities_opened = [y[i_f].X for i_f in range(self.F)]
        clients_assignment = \
            [i_f for i_c in range(self.C) for i_f in range(self.F) if x[i_f][i_c].X != 0]
        
        return facilities_opened, clients_assignment
    
    def solve_greedy(self):
        y = np.empty(self.F, dtype=np.int32)
        assignments = np.empty(self.C, dtype=np.int32)
        
        # your code here
        # ================
        for index, item in enumerate(self.assignment_costs.T):
            assignments[index] = np.where(item == np.min(item))[0][0]
            #print(index, '--', np.where(item == np.min(item))[0][0])
        
        print(assignments)
        
        # Open facilities based on the assigments
        for i in set(assignments):
            y[int(i)] = 1
        
        return y, assignments


# to make results reproducible, be sure to put this line before creating a Problem()
np.random.seed(666)

problem = Problem()
problem.assign_random_capacities()
problem.assign_random_opening_costs()

# completely random solution
# open_idx = list(range(problem.F))
# random.shuffle(open_idx)
# open_idx = open_idx[:problem.F * 4 // 5] # open 80% of facilities
#
# y = np.zeros(problem.F)
# y[open_idx] = 1
#
# from itertools import cycle
# assignments = np.empty(problem.C)
# for cli, fac in zip(range(problem.C), cycle(open_idx)):
#     assignments[cli] = fac
# print("completely random solution")
# print("----")
# print(assignments)
# print("END")
# ax = plt.figure().gca()
# problem.plot(y, assignments, ax)
# ax.set_title('Totally random solution, objective = %.4f' % problem.objective(y, assignments))
#
print("Gurobi")
ax = plt.figure().gca()
y, assignments = problem.solve_gurobi()
print(y)
print('------')
print(assignments)
problem.plot(y, assignments, ax)
ax.set_title('Most likely optimal solution, objective = %.4f' % problem.objective(y, assignments))


print("Greedy")
ax = plt.figure().gca()
y2, assignments2 = problem.solve_greedy()
print(y2)
print('------')
print(assignments2)
problem.plot(y2, assignments2, ax)
ax.set_title('Greedy solution, objective = %.4f' % problem.objective(y2, assignments2))