#C++ version
import numpy as np
import matplotlib.pyplot as plt
from sortedcontainers import *
from timeit import default_timer as timer
import seaborn as sns; sns.set() # for plot styling
from matplotlib.gridspec import GridSpec
def facility_location(clients, facilities, distances_list, c, opening_cost, draw=False):
"""
Implementation of algorithm outline found on page 183 of
"The Design of Approximation Algorithms" by Williamson and Shmoys,
and page 282 of "Approximation Algorithms for Metric Facility
Location and k-Median Problems Using the Primal-Dual Schema and
Lagrangian Relaxation" by Jain and Vazirani
Takes as input a lists representing the facilites and clients,
a (sorted) list of distances, and an integer for the facility opening cost.
Executes the primal-dual uncapacitated facility location algorithm
and finds a good approximate assignment for clients to facilities.
"""
print("Starting primal-dual algorithm")
start = timer()
# Variable c_ij -- the distance between facility i and client j
# Assume this is already in sorted order
client_distances = distances_list
num_client_distances = len(client_distances)
# Variable w -- how much a client contributes to a particular facility
# client_w[i][j] stores the time when client j starts to contribute to facility i
client_w = [[-1 for j in range(len(clients))] for i in range(len(facilities))]
# Record contribution when clients are declared connected
client_v = [-1 for j in range(len(clients))]
# We assume all facilitiy costs are equal
facility_opening_cost = [opening_cost] * len(facilities)
# Variable S -- gets a copy of the list of clients
# Each client stores a list of facilities that it's assigned to
# Coordinate 2 says if the client is currently open
clients_copy = [[x, []] for x in clients]
num_open_clients = len(clients_copy)
# Variable T -- starts empty, but will accept facilities
open_facilities = [-1] * len(facilities)
# Initial time each facility will be paid for is one more than lambda
# Second coordinate keeps track of the facility index
maximum_distance_plus_one = opening_cost + 1
facility_pay_schedule = [[maximum_distance_plus_one, x] for x in range(len(facilities))]
next_paid_facility = SortedList(facility_pay_schedule)
# Keep track of the index of the edge we service next
current_edge_index = 0
# Keep track of how much time has passed
current_time = 0
# Want to track how many times we update facility cost
counter = [0]
# Here we flag "closed" facilities with the value -1
while num_open_clients > 0:
if current_edge_index == num_client_distances:
print("All edges have been evaluated. Continuing algorithm")
# What if the algorithm needs more time?
while not next_paid_facility == []:
next_facility = next_paid_facility[0]
if next_facility[0] >= 0:
current_time = next_facility[0]
num_open_clients = update_facilities(next_facility[1], open_facilities[next_facility[1]][1],
open_facilities, facility_pay_schedule, facility_opening_cost,
next_paid_facility, clients_copy, num_open_clients, client_w, client_v, current_time)
if num_open_clients <= 0:
break
break
distance = client_distances[current_edge_index]
next_edge = distance[0]
# Decide which event happens next
if not next_paid_facility == []:
next_facility = next_paid_facility[0]
else:
# Don't go here!
next_facility = [next_edge + 1, 0]
if next_edge <= next_facility[0]:
# An edge goes tight
# if not clients_copy[distance[2]] == -1:
num_open_clients = service_tight_edge(open_facilities, facilities, clients_copy, distance, client_w, facility_opening_cost, facility_pay_schedule, next_paid_facility, num_open_clients, counter, client_v)
current_time = distance[0]
current_edge_index += 1
else:
# A facility is now paid for
if next_facility[0] >= 0:
current_time = next_facility[0]
num_open_clients = update_facilities(next_facility[1], open_facilities[next_facility[1]][1], open_facilities, facility_pay_schedule, facility_opening_cost, next_paid_facility, clients_copy, num_open_clients, client_w, client_v, current_time)
end = timer()
print("Finished primal-dual algorithm. Time:", end - start)
print("current_time:",current_time)
# Some facilities might be paid for now, but haven't officially opened.
# We fix that here.
for i in range(len(facilities)):
if not facility_pay_schedule[i] == -1:
# Update client contribution to facility i
# Only update pay from clients that are not yet connected
num_clients_facility_i = len(open_facilities[i][1]) #ここでエラーするので、try exceptにする？ TODO
update_client_pay = num_clients_facility_i * (current_time - open_facilities[i][2])
open_facilities[i][3] += update_client_pay
current_client_pay = open_facilities[i][3]
if current_client_pay >= facility_opening_cost[i]:
# Clients connected to facility i are removed from contributing to other facilities
num_open_clients = update_facilities(i, open_facilities[i][1], open_facilities, facility_pay_schedule, facility_opening_cost, next_paid_facility, clients_copy, num_open_clients, client_w, client_v, current_time)
pruning_primal_obj, pruning_opened_facilities, extra_cost = cluster_assignment_via_pruning(clients, facilities, client_w, c, facility_pay_schedule, client_v, opening_cost)
print("num open facilities:", len(pruning_opened_facilities))
# plot_final_points(clients, client_assignments)
for x in pruning_opened_facilities:
x[0] = facilities[x[0]]
if draw:
plot_final_clusters(pruning_opened_facilities)
end = timer()
print("Finsished! Time:", end - start)
def service_tight_edge(open_facilities, facilities, clients_copy, distance, client_w, facility_opening_cost, facility_pay_schedule, next_paid_facility, num_open_clients, counter, client_v):
"""
Subroutine for the facility location algorithm for the event that
an edge i,j goes tight
"""
i = distance[1]
j = distance[2]
# print("Checking connection", i, j, "at time", distance[0])
# Set the time for when client j starts to contribute
if not clients_copy[j] == -1:
client_w[i][j] = distance[0]
# check to see if client j neighbors facility i
# Assign and remove it if it does, provided the facility is paid for
if facility_pay_schedule[i] == -1:
if not clients_copy[j] == -1:
open_facilities[i][1].add(j)
num_open_clients = update_facilities(i, [j], open_facilities, facility_pay_schedule, facility_opening_cost, next_paid_facility, clients_copy, num_open_clients, client_w, client_v, distance[0])
# Otherwise, client j contributes to the cost of facility i
# Update contributions of other clients too
else:
if not clients_copy[j] == -1:
clients_copy[j][1] += [i]
# If facility i currently has no contributing clients
if open_facilities[i] == -1:
# Second coordinate -- track clients contributing to facility i
# Third coordinate tracks times of last cost update
# Fourth coordinate tracks running contribution of clients so far
open_facilities[i] = [facilities[i], SortedList(), distance[0], 0]
open_facilities[i][1].add(j)
# Update when facility i will be paid for.
next_paid_facility.remove([facility_pay_schedule[i][0], i])
facility_pay_schedule[i][0] = facility_opening_cost[i]
next_paid_facility.add([facility_pay_schedule[i][0], i])
# Facility i already has some clients
else:
# Update client contribution to facility i
num_clients_facility_i = len(open_facilities[i][1])
update_client_pay = num_clients_facility_i * (distance[0] - open_facilities[i][2])
open_facilities[i][3] += update_client_pay
# open_facilities[i][3] += new_client_pay
current_client_pay = open_facilities[i][3]
# print("Facility",i,"now will be paid in",facility_pay_schedule[i][0] - current_client_pay)
# Add client j to facility i
if not clients_copy[j] == -1:
open_facilities[i][1].add(j)
open_facilities[i][2] = distance[0]
if current_client_pay >= facility_opening_cost[i]:
# Clients connected to facility i are removed from contributing to other facilities
num_open_clients = update_facilities(i, open_facilities[i][1], open_facilities, facility_pay_schedule, facility_opening_cost, next_paid_facility, clients_copy, num_open_clients, client_w, client_v, distance[0])
else:
# Update when facility i will be paid for
counter[0] += 1
next_paid_facility.remove([facility_pay_schedule[i][0], i])
if not len(open_facilities[i][1]) == 0:
remaining_time = (facility_opening_cost[i] - current_client_pay) / len(open_facilities[i][1])
else:
# No currently contributing clients, so it should not get paid for anytime soon
remaining_time = facility_opening_cost[i] - distance[0] + 1
# print("new pay schedule is", facility_pay_schedule[i][0])
facility_pay_schedule[i][0] = distance[0] + remaining_time
next_paid_facility.add([facility_pay_schedule[i][0], i])
return num_open_clients
def update_facilities(i, update_clients, open_facilities, facility_pay_schedule, facility_opening_cost,
next_paid_facility, clients_copy, num_open_clients, client_w, client_v, current_time):
"""
Facility i is now paid for, so we declare its clients to be connected to it.
We only update clients specified in update_clients
"""
# iterate over clients contributing to facility i
clients_to_add = []
for j in update_clients:
if not clients_copy[j] == -1:
# declare client j to be connected
# print("Client",j,"connects to facility", i,"at time", current_time)
# Remove client j from contributing to anymore facilities
# Update when other facilities are going to be paid for
for f in clients_copy[j][1]:
if not open_facilities[f] == -1 and not f == i:
if not facility_pay_schedule[f] == -1:
if not client_w[f][j] == -1:
next_paid_facility.remove([facility_pay_schedule[f][0], f])
# print("Updating", f, "old val", facility_pay_schedule[f][0])
# Record and update contribution of client j to facility f
num_clients_facility_f = len(open_facilities[f][1])
update_client_pay = num_clients_facility_f * (current_time - open_facilities[f][2])
open_facilities[f][3] += update_client_pay
current_client_pay = open_facilities[f][3]
facility_pay_schedule[f][0] = facility_opening_cost[f] - current_client_pay
open_facilities[f][2] = current_time
if not len(open_facilities[f][1]) == 0:
if j in open_facilities[f][1]:
if len(open_facilities[f][1]) > 1:
remaining_time = (facility_opening_cost[f] - current_client_pay) / (len(open_facilities[f][1]) - 1)
else:
remaining_time = facility_opening_cost[f] - current_time + 1
else:
remaining_time = (facility_opening_cost[f] - current_client_pay) / len(open_facilities[f][1])
else:
# No currently contributing clients, so it should not get paid for anytime soon
remaining_time = facility_opening_cost[f] - current_time + 1
# print("new pay schedule is", facility_pay_schedule[f][0])
facility_pay_schedule[f][0] = current_time + remaining_time
next_paid_facility.add([facility_pay_schedule[f][0], f])
# print("new val", facility_pay_schedule[f][0])
# remove all copies of client j from facility f
open_facilities[f][1].remove(j)
if num_open_clients == len(facility_opening_cost):
print("Time of first facility opening", current_time)
client_v[j] = [current_time, i] # record time and facility when client j is connected
clients_copy[j] = -1
num_open_clients -= 1
# Remove facility i from the list of facilities waiting to be paid for
if not facility_pay_schedule[i] == -1:
next_paid_facility.remove([facility_pay_schedule[i][0], i])
facility_pay_schedule[i] = -1
# print("Facility",i,"is now open")
return num_open_clients
def cluster_assignment_via_pruning(clients, facilities, client_w, c, facility_pay_schedule, client_v, opening_cost):
"""
Implementing the pruning procedure
"""
primal_obj = 0
assigned_clients = []
# Candidate facilities to remain open
temp_open_facilities = [i for i in range(len(facility_pay_schedule)) if facility_pay_schedule[i] == -1]
opened_facilities = []
while temp_open_facilities != []:
i = temp_open_facilities[0]
temp_open_facilities = temp_open_facilities[1:] # slice off current facility
facility_assignment = [i, []]
for j in range(len(client_v)):
# How do I tell if client j ever contributed to facility i? Should I check if client_w[i][j] != -1?
witness = client_v[j][1]
# Hopefully client_w is up-to-date
witness_cost = (client_v[j][0] - client_w[witness][j]) + c[witness][j]
current_cost = (client_v[j][0] - client_w[i][j]) + c[i][j]
if current_cost <= witness_cost and client_w[i][j] != -1:
primal_obj += c[i][j]
assigned_clients += [j]
facility_assignment[1] += [clients[j]]
temp_open_facilities = [h for h in temp_open_facilities if not (client_w[i][j] > 0 and client_w[h][j] > 0)]
if not len(facility_assignment[1]) == 0:
opened_facilities += [facility_assignment]
# Continuing Method 2: assigned clients not yet connected to their nearest open facility
extra_cost = 0
for j in range(len(client_v)):
if not j in assigned_clients:
nearest_facility_list = [[c[opened_facilities[i][0]][j], i, opened_facilities[i][0]] for i in range(len(opened_facilities))]
nearest_facility_sublist = min(nearest_facility_list)
nearest_facility_index = nearest_facility_sublist[1]
primal_obj += c[nearest_facility_sublist[2]][j]
extra_cost += c[nearest_facility_sublist[2]][j]
opened_facilities[nearest_facility_index][1] += [clients[j]]
primal_obj += len(opened_facilities) * opening_cost
print("Pruning Assigned Clients", len(assigned_clients))
print("Pruning Open Facilities", len(opened_facilities))
return primal_obj, opened_facilities, extra_cost
def plot_final_points(clients, client_assignments):
"""
Here the clients are grouped by open facility
"""
facility_centers = []
for j in range(len(client_assignments)):
plt.plot([clients[j][0]], [clients[j][1]], marker='o', color='red')
# iterate over clients assigned to this facility
if j == 0:
for k in client_assignments[j]:
facility_centers += [[clients[k[0]][0], clients[k[0]][1]]]
# plot the facilities
for i in facility_centers:
plt.plot([i[0]], [i[1]], marker='s', color='black')
plt.show()
def plot_final_clusters(open_facilities):
"""
Here the clients are grouped by open facility
"""
ax = plt.gca()
facility_centers = []
for i in open_facilities:
color = next(ax._get_lines.prop_cycler)['color']
# iterate over clients assigned to this facility
for j in i[1]:
plt.plot([j[0]], [j[1]], marker='o', color=color)
facility_centers += [[i[0][0], i[0][1]]]
# plot the facilities
for i in facility_centers:
plt.plot([i[0]], [i[1]], marker='s', color='black')
plt.show()容量制約なし施設配置問題に対する主双対アルゴリズム
import numpy as np
from numpy import linalg
import matplotlib.pyplot as plt
from sortedcontainers import *
from timeit import default_timer as timer
from scipy.spatial import KDTree
import csv
# Import that came from the point clusters code
# %matplotlib inline
import seaborn as sns; sns.set() # for plot styling
from sklearn.datasets import make_blobs
from sklearn.datasets import make_moons
from sklearn.datasets import make_circles
# Methods to generate starting data
def compute_distances(clients):
"""
Compute the distance between each facility i and client j
Return a sorted list of where each entry has the form
[distance, facility_index, client_index]
"""
print("Starting to compute distances")
start = timer()
client_distances = []
counter = [0 for x in range(len(clients))]
biggest_edge = [0 for x in range(len(clients))]
c = [[0 for j in range(len(clients))] for i in range(len(clients))]
for i in range(len(clients)):
for j in range(i, len(clients)):
current_distance = np.linalg.norm(clients[i]-clients[j])
client_distances += [[current_distance, i, j]]
counter[i] += current_distance
c[i][j] = current_distance
if biggest_edge[i] < current_distance:
biggest_edge[i] = current_distance
# Distance is symmetric
if not i == j:
client_distances += [[current_distance, j, i]]
counter[j] += current_distance
c[j][i] = current_distance
if biggest_edge[j] < current_distance:
biggest_edge[j] = current_distance
client_distances = SortedList(client_distances)
end = timer()
results = [len(clients) * biggest_edge[i] - counter[i] for i in range(len(clients))]
print("Finished computing distances. Time:", end-start)
return client_distances, c
def compute_distances2(facilities, clients):
"""
Compute the distance between each facility i and client j
Return a sorted list of where each entry has the form
[distance, facility_index, client_index]
"""
print("Starting to compute distances")
c = [[0 for i in range(len(clients))] for i in range(len(facilities))]
start = timer()
client_distances = []
counter = [0 for x in range(len(clients))]
biggest_edge = [0 for x in range(len(clients))]
for i in range(len(facilities)):
for j in range(len(clients)):
current_distance = np.linalg.norm(facilities[i]-clients[j])
client_distances += [[current_distance, i, j]]
counter[i] += current_distance
c[i][j] = current_distance
if biggest_edge[i] < current_distance:
biggest_edge[i] = current_distance
# Distance is symmetric
# if not i == j:
# client_distances += [[current_distance, j, i]]
# counter[j] += current_distance
# #c[j][i] = current_distance
# if biggest_edge[j] < current_distance:
# biggest_edge[j] = current_distance
client_distances = SortedList(client_distances)
end = timer()
results = [len(clients) * biggest_edge[i] - counter[i] for i in range(len(clients))]
print("Finished computing distances. Time:", end-start)
return client_distances, c
def compute_distances3(facilities, clients, theta):
"""
Compute the distance between each facility i and client j
Return a sorted list of where each entry has the form
[distance, facility_index, client_index]
"""
print("Starting to compute distances")
c = [[0 for i in range(len(clients))] for i in range(len(facilities))]
start = timer()
client_distances = []
counter = [0 for x in range(len(clients))]
biggest_edge = [0 for x in range(len(clients))]
tree = KDTree(facilities)
dis, idx = tree.query(clients, k=theta)#各クライアントに近い100の施設
for j in range(len(clients)):
for _i in range(theta):
i = idx[j][_i]
current_distance = dis[j][_i]
client_distances += [[current_distance, i, j]]
counter[i] += current_distance
c[i][j] = current_distance
if biggest_edge[i] < current_distance:
biggest_edge[i] = current_distance
client_distances = SortedList(client_distances)
end = timer()
results = [len(clients) * biggest_edge[i] - counter[i] for i in range(len(clients))]
print("Finished computing distances. Time:", end-start)
return client_distances, c
def plot_initial_points(clients):
"""
Plot facilities and clients on the same graph
"""
client_coords = [[],[]]
for j in clients:
client_coords[0] += [j[0]]
client_coords[1] += [j[1]]
# Clients are red circles
plt.plot(client_coords[0], client_coords[1], 'ro')
plt.show()
def make_data_set(num_samples, num_facilities, X, theta=100):
"""
Prepare data set X for primal-dual algorithm
theta: number of near neighbor facilities for each client
"""
clients = [np.array((X[i][0], X[i][1])) for i in range(num_samples)]
facilities = np.random.permutation(clients)[:num_facilities]
#distances, c = compute_distances2(facilities, clients)
#compute_distances2(facilities, clients, theta)
distances, c = compute_distances(clients)
# plot_initial_points(clients)
return clients, facilities, distances, c
def cluster_data(num_clusters, num_samples, var=0.40, theta=100):
X, y_true = make_blobs(n_samples=num_samples, centers=num_clusters, cluster_std=var, random_state=0)
return make_data_set(num_samples, num_facilities, X, theta)
def crescent_moons_data(num_samples, theta):
X, y_true = make_moons(num_samples, noise=.05, random_state=0)
return (num_samples, num_facilities, X, theta)
def noisy_circles_data(num_samples, theta):
X, y_true = make_circles(num_samples, factor=.5, noise=.05, random_state=0)
return (num_samples, num_facilities, X, theta)
def anisotropic_data(num_samples, theta):
X, y = make_blobs(n_samples=num_samples, random_state=170)
transformation = [[0.6, -0.6], [-0.4, 0.8]]
X_aniso = np.dot(X, transformation)
return make_data_set(num_samples, num_facilities, X_aniso, theta)
def varied_variance_data(num_samples, theta):
X, y = make_blobs(n_samples=num_samples, cluster_std=[0.4, 1.2, 0.8], random_state=170)
return (num_samples, num_facilities, X, theta)
def read_data_set(file_name):
"""
Prepare data set X for primal-dual algorithm
"""
# read dataset from file
f = open(file_name, "r")
lines = f.readlines()
points_list = []
for line in lines:
points = line.split(",")
points_list += [[float(points[0]), float(points[1][:-1])]]
f.close()
return make_data_set(len(points_list), points_list)
def write_sorted_distances(distances, file_name):
"""
Write sorted dataset to file
"""
num_samples = len(distances)
# Write dataset to file
f = open(file_name, "w+")
for i in range(num_samples - 1):
f.write(str(distances[i][0]) + "," + str(distances[i][1]) + "," + str(distances[i][2]) + "\n")
f.write(str(distances[i][0]) + "," + str(distances[i][1]) + "," + str(distances[i][2]))
f.close()num_points = 150
num_facilities =100
theta =100
cost= 10
data_name = "cluster"
# data_name = "moon"
# data_name = "circle"
# data_name = "aniso"
# data_name = "variedvar"
if data_name == "cluster":
clients, facilities, distances, c = cluster_data(4, num_points, num_facilities, theta)
elif data_name == "moon":
clients, facilities,distances, c = crescent_moons_data(num_points)
elif data_name == "circle":
clients, facilities,distances, c = noisy_circles_data(num_points)
elif data_name == "aniso":
clients, facilities,distances, c = anisotropic_data(num_points)
elif data_name == "variedvar":
clients, facilities, distances, c = varied_variance_data(num_points)
facility_location(clients, facilities, distances, c, cost, draw=True)Starting to compute distances
Finished computing distances. Time: 0.0731076199990639
Starting primal-dual algorithm--------------------------------------------------------------------------- IndexError Traceback (most recent call last) Cell In[34], line 22 19 elif data_name == "variedvar": 20 clients, facilities, distances, c = varied_variance_data(num_points) ---> 22 facility_location(clients, facilities, distances, c, cost, draw=True) Cell In[32], line 96, in facility_location(clients, facilities, distances_list, c, opening_cost, draw) 92 next_facility = [next_edge + 1, 0] 93 if next_edge <= next_facility[0]: 94 # An edge goes tight 95 # if not clients_copy[distance[2]] == -1: ---> 96 num_open_clients = service_tight_edge(open_facilities, facilities, clients_copy, distance, client_w, facility_opening_cost, facility_pay_schedule, next_paid_facility, num_open_clients, counter, client_v) 97 current_time = distance[0] 98 current_edge_index += 1 Cell In[32], line 153, in service_tight_edge(open_facilities, facilities, clients_copy, distance, client_w, facility_opening_cost, facility_pay_schedule, next_paid_facility, num_open_clients, counter, client_v) 149 # print("Checking connection", i, j, "at time", distance[0]) 150 151 # Set the time for when client j starts to contribute 152 if not clients_copy[j] == -1: --> 153 client_w[i][j] = distance[0] 155 # check to see if client j neighbors facility i 156 # Assign and remove it if it does, provided the facility is paid for 157 if facility_pay_schedule[i] == -1: IndexError: list index out of range
#KD-tree
# tree = KDTree(facilities)
# dis, idx = tree.query(clients, k=100)#各クライアントに近い100の施設
# dis[0], idx[0]len(clients)1500実験結果
| クライアント数 | 施設数 | 距離計算 (s) | 主双対法 (s) |
|---|---|---|---|
| 10_000 | 1_000 | 66 | 16 |
| 100_000 | 1_000 | 900 | 239 |
| 150_000 | 1_000 | 1785 | 952 |
// a C++ implementation of the primal-dual clustering algorithm
#include "algorithm"
#include "array"
#include "chrono"
#include "cmath"
#include "fstream"
#include "iostream"
#include "string"
#include "vector"
using namespace std;
vector<array<double, 2> > parseCSV(const char* filename) {
// Reads a text file, inserts data points into a vector
ifstream file(filename);
string value_1;
string value_2;
vector<array<double, 2> > data_points;
while (file.good()) {
getline(file, value_1, ',');
getline(file, value_2, '\n');
if (!value_1.empty()) {
array<double, 2> data_pt = {stod(value_1), stod(value_2)};
data_points.push_back(data_pt);
}
}
return data_points;
}
vector<array<double, 3> > parseCSVSorted(const char* filename) {
// Reads a text file, inserts data points into a vector
ifstream file(filename);
string value_1;
string value_2;
string value_3;
vector<array<double, 3> > data_points;
while (file.good()) {
getline(file, value_1, ',');
getline(file, value_2, ',');
getline(file, value_3, '\n');
// This seems to be necessary here...
if (!value_1.empty()) {
array<double, 3> data_pt = {stod(value_1), stod(value_2), stod(value_3)};
data_points.push_back(data_pt);
}
}
return data_points;
}
vector<vector<double> > getDistanceGrid(vector<array<double, 3> > &distanceList, int numClients) {
// put distances from sorted list into grid form
vector<vector<double> > distanceGrid;
vector<double> current_vector;
current_vector.resize(numClients, 0);
distanceGrid.resize(numClients, current_vector);
int x; int y; double distance;
// cout << "Distance List Size" << distanceList.size() << endl;
for (int i = 0; i < distanceList.size(); i++) {
x = (int) distanceList.at(i)[1];
y = (int) distanceList.at(i)[2];
distance = distanceList.at(i)[0];
distanceGrid.at(x).at(y) = distance;
}
return distanceGrid;
}
// -------- Methods for computing and sorting distances ---------
// http://www.cplusplus.com/forum/beginner/178293/
double euclideanDistance(double x1, double y1, double x2, double y2) {
double x = x1 - x2;
double y = y1 - y2;
double dist;
dist = pow(x, 2) + pow(y, 2);
dist = sqrt(dist);
return dist;
}
// Used for sorting distance vector
bool compareArrays(array<double, 3> a1, array<double, 3> a2) {
return (a1[0] < a2[0]);
}
// I hope this will create a "reverse" sort, to be able to use the vector pop() operation
bool compareArraysTwo(array<double, 2> a1, array<double, 2> a2) {
return (a1[0] > a2[0]);
}
vector<array<double, 3> > computeDistances(vector<array<double, 2> > clients){
/*
* Compute the distance between each facility i and client j
* Return a sorted list of where each entry has the form
* {distance, facility_index, client_index}
*/
cout << "Starting to compute distances" << endl;
vector<array<double, 3> > client_distances;
vector<array<double, 2> >::const_iterator i;
vector<array<double, 2> >::const_iterator j;
double current_distance;
double x1; double x2; double y1; double y2;
array<double, 2> entry_1;
array<double, 2> entry_2;
double facility_index = 0;
double client_index = 0;
for(i = clients.begin(); i != clients.end(); ++i) {
for (j = i; j != clients.end(); ++j) {
entry_1 = *i; entry_2 = *j;
x1 = entry_1[0]; x2 = entry_1[1];
y1 = entry_2[0]; y2 = entry_2[1];
current_distance = euclideanDistance(x1, x2, y1, y2);
array<double, 3> current_entry = {current_distance, facility_index, client_index};
client_distances.push_back(current_entry);
// Distance is symmetric
if (facility_index != client_index) {
current_entry = {current_distance, client_index, facility_index};
client_distances.push_back(current_entry);
}
client_index++;
}
facility_index++;
client_index = facility_index;
}
sort(client_distances.begin(), client_distances.end(), compareArrays);
cout << "Finished computing distances" << endl;
return client_distances;
}
// https://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c
bool file_exists (const std::string& name) {
ifstream f(name.c_str());
return f.good();
}
vector<array<double, 3> > retrieveDistances (const char* filename, vector<array<double, 2> > &data_points) {
vector<array<double, 3> > distances;
string sorted_filename(filename);
sorted_filename.insert(sorted_filename.size()-4,"_sorted");
// Check to see if data set has already been computed
if (file_exists(sorted_filename)) {
// File already exists
distances = parseCSVSorted(sorted_filename.c_str());
}
else {
// compute distances and store results
auto start = chrono::high_resolution_clock::now();
distances = computeDistances(data_points);
auto finish = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed = finish - start;
cout << "Elapsed time: " << elapsed.count() << " s\n";
// Write distances to file
ofstream output_file;
output_file.open(sorted_filename.c_str());
array<double, 3> assignment;
vector<array<double, 3> >::const_iterator results;
for (results = distances.begin(); results < distances.end(); results++) {
assignment = *results;
output_file << assignment[0] << "," << assignment[1] << "," << assignment[2] << "\n";
}
output_file.close();
}
return distances;
}
// -------- Methods for the primal dual algorithm ----------
bool openClient (vector<int> client) {
// Need to avoid looking up elements in empty vectors
if (!client.empty()) {
// if (client.back() == -1) {
if (find(client.begin(), client.end(), -1) != client.end()) {
return false;
}
}
return true;
}
int update_facilities(int i, vector<int> &update_clients, vector<vector<int> > &open_facilities, vector<array<double, 2> > &facility_pay_schedule, vector<array<double, 2> > &next_paid_facility, vector<vector<int> > &clients_copy, vector<vector<double> > &client_w, int num_open_clients, double current_time, vector<double> &client_contribute_times, vector<double> &facility_contributions, vector<double> &facility_opening_cost, vector<array<double, 2> > &client_v) {
/*
* Facility i is now paid for, so we now prevent all clients that are contributing
* to facility i from contributing to any other facilities.
*
* We only update clients specified in update_clients
*/
// cout << "update facilities" << endl;
// iterate over clients contributing to facility i
vector<int> clients_to_add;
vector<int>::const_iterator client_index;
vector<int> assigned_facilities;
vector<int>::const_iterator facilities_index;
vector<int> assigned_clients;
vector<int>::iterator client_index_2;
array<double, 2> old_value;
array<double, 2> new_value;
double num_clients_facility_f;
double update_client_pay;
double current_client_pay;
double remaining_time;
int j; int f; int h;
if (!openClient(update_clients)) {
return num_open_clients;
}
for (client_index = update_clients.begin(); client_index < update_clients.end(); client_index++) {
j = *client_index;
// cout << "client index is: " << j << endl;
if (j == -1) {
break;
}
assigned_facilities = clients_copy.at(j);
// iterate over facilities that client j contributes to
if (openClient(assigned_facilities)) {
for (facilities_index = assigned_facilities.begin(); facilities_index < assigned_facilities.end(); facilities_index++) {
f = *facilities_index;
//assigned_clients = open_facilities.at(f);
if (openClient(open_facilities.at(f)) && f != i) {
if (facility_pay_schedule.at(f)[0] != -1) {
if (client_w.at(f).at(j) != -1) {
num_clients_facility_f = open_facilities.at(f).size();
update_client_pay = num_clients_facility_f * (current_time - client_contribute_times.at(f));
old_value = {facility_pay_schedule.at(f)[0], (double) f};
facility_contributions.at(f) += update_client_pay;
current_client_pay = facility_contributions.at(f);
facility_pay_schedule.at(f)[0] -= facility_contributions.at(f);
client_contribute_times.at(f) = current_time;
if (open_facilities.at(f).size() != 0) {
if(find(open_facilities.at(f).begin(), open_facilities.at(f).end(), j) != open_facilities.at(f).end()) {
if (open_facilities.at(f).size() > 1) {
remaining_time = (facility_opening_cost.at(f) - current_client_pay) / (open_facilities.at(f).size() - 1);
}
else {
remaining_time = facility_opening_cost.at(f) - current_time + 1;
}
}
else {
remaining_time = (facility_opening_cost.at(f) - current_client_pay) / (open_facilities.at(f).size());
}
}
else {
remaining_time = facility_opening_cost.at(f) - current_time + 1;
}
facility_pay_schedule.at(f)[0] = current_time + remaining_time;
new_value = {facility_pay_schedule.at(f)[0], (double) f};
// TODO: make this replace & sort more efficient.
replace(next_paid_facility.begin(), next_paid_facility.end(), old_value, new_value);
sort(next_paid_facility.begin(), next_paid_facility.end(), compareArraysTwo);
}
// remove all copies of client j from facility f
// TODO: make this a log(n) operation
open_facilities.at(f).erase(remove(open_facilities.at(f).begin(), open_facilities.at(f).end(), j), open_facilities.at(f).end());
}
// set contribution of client j to facility f equal to 0
// client_w.at(f).at(j) = 0;
}
}
clients_copy.at(j).push_back(-1);
client_v.at(j) = {current_time, (double) i};
num_open_clients -= 1;
// cout << "Current Number of Open Clients: " << num_open_clients << endl;
}
}
// Remove facility i from the list of facilities waiting to be paid for
if (facility_pay_schedule.at(i)[0] != -1) {
old_value = {facility_pay_schedule.at(i)[0], (double) i};
// TODO: make this a log(n) operation
next_paid_facility.erase(remove(next_paid_facility.begin(), next_paid_facility.end(), old_value), next_paid_facility.end());
facility_pay_schedule.at(i)[0] = -1;
}
return num_open_clients;
}
int service_tight_edge(vector<vector<int> > &open_facilities, vector<array<double, 2> > &facilities, vector<vector<int> > &clients_copy, array<double, 3> distance, vector<vector<double> > &client_w, vector<double> &facility_opening_cost, vector<array<double, 2> > &facility_pay_schedule, vector<array<double, 2> > &next_paid_facility, int num_open_clients, vector<double> &client_contribute_times, vector<double> &facility_contributions, vector<array<double, 2> > &client_v){
/*
* Subroutine for the facility location algorithm for the event that
* an edge i,j goes tight
*/
// cout << "service tight edge" << endl;
vector<int> update_clients;
array<double, 2> old_value;
array<double, 2> new_value;
double update_client_pay;
double current_client_pay;
int num_clients_facility_i;
int i = (int) distance[1];
int j = (int) distance[2];
// set the time for when client j starts to contribute
if (openClient(clients_copy.at(j))) {
client_w.at(i).at(j) = distance[0];
}
// check to see if client j neighbors facility i
// Assign and remove it if it does, provided the facility is paid for
if (facility_pay_schedule.at(i)[0] == -1) {
open_facilities.at(i).push_back(j);
update_clients.push_back(j);
num_open_clients = update_facilities(i, update_clients, open_facilities, facility_pay_schedule, next_paid_facility, clients_copy, client_w, num_open_clients, distance[0], client_contribute_times, facility_contributions, facility_opening_cost, client_v);
}
// Otherwise, client j contributes to the cost of facility i
// Update contributions of other clients too
else {
if (openClient(clients_copy.at(j))) {
// Add facility i to client j
clients_copy.at(j).push_back(i);
}
// If facility i currently has no contributing clients
if (!openClient(open_facilities.at(i))) {
// Add client j to facility i. Find and replace the -1
replace(open_facilities.at(i).begin(), open_facilities.at(i).end(), -1, j);
// Track time last contribution update to facility i
client_contribute_times.at(i) = distance[0];
// Initial facility contribution is already set to 0
// Update when facility i will be paid for
old_value = {facility_pay_schedule.at(i)[0], (double) i};
facility_pay_schedule.at(i)[0] = facility_opening_cost.at(i);
new_value = {facility_pay_schedule.at(i)[0], (double) i};
// TODO: make this replace & sort more efficient.
replace(next_paid_facility.begin(), next_paid_facility.end(), old_value, new_value);
sort(next_paid_facility.begin(), next_paid_facility.end(), compareArraysTwo);
}
// Facility i already has some clients
else {
// Update client contribution to facility i
num_clients_facility_i = open_facilities.at(i).size();
update_client_pay = num_clients_facility_i * (distance[0] - client_contribute_times.at(i));
facility_contributions.at(i) += update_client_pay;
current_client_pay = facility_contributions.at(i);
// Add client j to facility i
if (openClient(clients_copy.at(j))) {
open_facilities.at(i).push_back(j);
num_clients_facility_i++;
}
client_contribute_times.at(i) = distance[0];
// sort(client_contribute_times.at(i).begin(), client_contribute_times.at(i).end());
// Check to see if facility is paid for
if (current_client_pay >= facility_opening_cost.at(i)) {
// Clients connected to facility i are removed from contributing to other facilities
num_open_clients = update_facilities(i, open_facilities.at(i), open_facilities, facility_pay_schedule, next_paid_facility, clients_copy, client_w, num_open_clients, distance[0], client_contribute_times, facility_contributions, facility_opening_cost, client_v);
}
else {
// Update when facility i will be paid for
old_value = {facility_pay_schedule.at(i)[0], (double) i};
double remaining_time;
if (open_facilities.at(i).size() != 0) {
remaining_time = (facility_opening_cost.at(i) - current_client_pay) / num_clients_facility_i;
}
else {
remaining_time = facility_opening_cost.at(i) - distance[0] + 1;
}
facility_pay_schedule.at(i)[0] = distance[0] + remaining_time;
new_value = {facility_pay_schedule.at(i)[0], (double) i};
// TODO: make this replace & sort more efficient.
replace(next_paid_facility.begin(), next_paid_facility.end(), old_value, new_value);
sort(next_paid_facility.begin(), next_paid_facility.end(), compareArraysTwo);
}
}
}
return num_open_clients;
}
vector<vector<int> > output_cluster_results(vector<array<double, 2> > &clients, vector<vector<double> > &client_w, vector<vector<double> > &distanceGrid, vector<array<double, 2> > &facility_pay_schedule, vector<vector<int> > &open_facilities, vector<array<double, 2> > &client_v, double current_time) {
// I'm going to check contributions based off of current_time - client_w (e.g. end time minus start time)
vector<int> assigned_clients;
vector<int> temp_open_facilities;
for (int i = 0; i < facility_pay_schedule.size(); i++) {
if (facility_pay_schedule.at(i)[0] == -1)
temp_open_facilities.push_back(i);
}
vector<vector<int> > opened_facilities;
int i;
while (!temp_open_facilities.empty()) {
// get and remove a facility
i = temp_open_facilities.back();
temp_open_facilities.pop_back();
// first coordinate is facility, rest are indices of assigned clients
vector<int> facility_assignment;
facility_assignment.push_back(i);
for (int j = 0; j < client_v.size(); j++) {
//int current_client = open_facilities.at(i).at(j);
// double current_contribution = current_time - client_w.at(i).at(current_client);
int current_client = j;
// need to use time between when client starts contributing to facility
// and when client has its contributions capped
int witness = (int) client_v.at(current_client)[1];
double witness_cost = (client_v.at(current_client)[0] - client_w.at(witness).at(current_client)) + distanceGrid.at(witness).at(current_client);
double current_cost = (client_v.at(current_client)[0] - client_w.at(i).at(current_client)) + distanceGrid.at(i).at(current_client);
if (current_cost <= witness_cost && client_w.at(i).at(current_client) != -1) {
assigned_clients.push_back(current_client);
facility_assignment.push_back(current_client);
vector<int> facilities_to_remove;
// track other potential open facilities that are counting on contribution from j
for (int h = 0; h < temp_open_facilities.size(); h++) {
int current_facility = temp_open_facilities.at(h);
vector<int> current_client_list = open_facilities.at(current_facility);
if (client_w.at(i).at(current_client) > 0 && client_w.at(current_facility).at(current_client) > 0) {
// blacklist the current facility
facilities_to_remove.push_back(current_facility);
}
}
// erase the offending facilities
for (int h = 0; h < facilities_to_remove.size(); h++) {
temp_open_facilities.erase(remove(temp_open_facilities.begin(), temp_open_facilities.end(), facilities_to_remove.at(h)), temp_open_facilities.end());
}
}
}
if (facility_assignment.size() > 1) {
opened_facilities.push_back(facility_assignment);
}
}
for (int j = 0; j < clients.size(); j++) {
if (find(assigned_clients.begin(), assigned_clients.end(), j) == assigned_clients.end()) {
// cout << "Assigning client " << j << endl;
int current_facility;
double current_distance;
double min_distance = -1;
int min_facility_index;
for (int h = 0; h < opened_facilities.size(); h++) {
current_facility = opened_facilities.at(h).at(0);
current_distance = distanceGrid.at(current_facility).at(j);
if (min_distance == -1) {
min_distance = current_distance;
min_facility_index = h;
}
else if (current_distance < min_distance) {
min_distance = current_distance;
min_facility_index = h;
}
}
opened_facilities.at(min_facility_index).push_back(j);
}
}
return opened_facilities;
}
int facility_location(vector<array<double, 2> > clients, vector<array<double, 3> > client_distances, double opening_cost, string file_name){
/*
* Implementation of algorithm outline found on page 183 of
* "The Design of Approximation Algorithms" by Williamson and Shmoys,
* and page 282 of "Approximation Algorithms for Metric Facility
* Location and k-Median Problems Using the Primal-Dual Schema and
* Lagrangian Relaxation" by Jain and Vazirani
*
* Takes as input a lists representing the facilites and clients,
* a (sorted) list of distances, and an integer for the facility opening cost.
* Executes the primal-dual uncapacitated facility location algorithm
* and finds a good approximate assignment for clients to facilities.
*/
array<double, 3> distance;
array<double, 2> next_facility;
double next_edge;
cout << "Starting primal-dual algorithm" << endl;
auto start_2 = chrono::high_resolution_clock::now();
// note that facilities and clients are the same here
int num_client_distances = client_distances.size();
int num_clients = clients.size();
// obtain the grid of distances c[i][j] to use for clustering assignments
vector<vector<double> > distanceGrid = getDistanceGrid(client_distances, num_clients);
// Variable w -- how much a client contributes to a particular facility
// client_w[i][j] stores the time when client j starts to contribute to facility i
vector<vector<double> > client_w;
vector<double> current_vector;
current_vector.resize(num_clients, -1);
client_w.resize(num_clients, current_vector);
// Record time and facility when clients are declared connected
vector<array<double, 2> > client_v;
client_v.resize(num_clients, {-1, -1});
// We assume all facilitiy costs are equal
vector<double> facility_opening_cost;
facility_opening_cost.resize(num_clients, opening_cost);
// Variable S -- gets a copy of the list of clients
// Each client stores a list of facilities that it's assigned to
vector<vector<int> > clients_copy;
clients_copy.resize(num_clients);
int num_open_clients = clients_copy.size();
// Variable T -- starts empty, but will accept facilities
vector<vector<int> > open_facilities;
vector<int> initial_list_open_facilities;
initial_list_open_facilities.push_back(-1);
open_facilities.resize(num_clients, initial_list_open_facilities);
// Want to also track when clients start to contribute to a given facility
vector<double> client_contribute_times;
client_contribute_times.resize(num_clients, 0);
// We will also track how much is currently being contributed to a given facility
vector<double> facility_contributions;
facility_contributions.resize(num_clients, 0);
// Initial time each facility will be paid for is one more than the max edge length
// Second coordinate keeps track of the facility index
// double maximum_distance = client_distances.at(num_client_distances - 1)[0];
vector<array<double, 2> > facility_pay_schedule;
facility_pay_schedule.resize(num_clients);
for (int i = 0; i < num_clients; i++){
array<double, 2> initial_facility_pay = {opening_cost + 1, (double) i};
facility_pay_schedule.at(i) = initial_facility_pay;
}
vector<array<double, 2> > next_paid_facility = facility_pay_schedule;
// Keep track of the index of the edge we service next
int current_edge_index = 0;
// Keep track of how much time has passed
double current_time = 0;
// Here we flag "closed" clients and facilities with the value -1
while (num_open_clients > 0) {
if (current_edge_index == num_client_distances) {
cout << "All edges have been evaluated. Continuing algorithm" << endl;
while (!next_paid_facility.empty()) {
next_facility = next_paid_facility.back();
next_paid_facility.pop_back();
if (next_facility[0] >= 0) {
current_time = next_facility[0];
}
num_open_clients = update_facilities(next_facility[1], open_facilities.at(next_facility[1]), open_facilities, facility_pay_schedule, next_paid_facility, clients_copy, client_w, num_open_clients, current_time, client_contribute_times, facility_contributions, facility_opening_cost, client_v);
if (num_open_clients <= 0) {
break;
}
}
break;
}
// cout << "Looking at edge " << current_edge_index << endl;
distance = client_distances.at(current_edge_index);
next_edge = distance[0];
// Decide which event happens next
if (!next_paid_facility.empty()) {
next_facility = next_paid_facility.back();
next_paid_facility.pop_back();
}
else {
// Don't go here!
next_facility = {next_edge + 1, 0};
}
if (next_edge <= next_facility[0]) {
// An edge goes tight
// Save some time by checking this first
// if (openClient(clients_copy.at(distance[2]))) {
num_open_clients = service_tight_edge(open_facilities, clients, clients_copy, distance, client_w, facility_opening_cost, facility_pay_schedule, next_paid_facility, num_open_clients, client_contribute_times, facility_contributions, client_v);
// }
current_time = distance[0];
current_edge_index += 1;
}
else {
// A facility is now paid for
num_open_clients = update_facilities(next_facility[1], open_facilities.at(next_facility[1]), open_facilities, facility_pay_schedule, next_paid_facility, clients_copy, client_w, num_open_clients, current_time, client_contribute_times, facility_contributions, facility_opening_cost, client_v);
current_time = next_facility[0];
}
// cout << "Current number of open clients: " << num_open_clients << endl;
}
auto finish_2 = chrono::high_resolution_clock::now();
chrono::duration<double> elapsed_2 = finish_2 - start_2;
cout << "Finished primal-dual algorithm" << endl;
cout << "Elapsed time: " << elapsed_2.count() << " s\n";
// Get results of pruning step
vector<vector<int> > opened_facilities = output_cluster_results(clients, client_w, distanceGrid, facility_pay_schedule, open_facilities, client_v, current_time);
vector<vector<int> >::const_iterator results;
vector<int> assignment;
vector<int>::const_iterator results_2; // list of clients assigned to facility
int client_assigned_to_facility;
int num_facilities_opened = opened_facilities.size();
ofstream output_file;
output_file.open("pd_result_"+file_name);
// output algorithm results
for (results = opened_facilities.begin(); results < opened_facilities.end(); results++) {
assignment = *results;
output_file << clients.at(*assignment.begin())[0] << "," << clients.at(*assignment.begin())[1] << endl;
// both facility and assigned clients are written here
for (results_2 = assignment.begin(); results_2 < assignment.end(); results_2++) {
client_assigned_to_facility = *results_2;
output_file << clients.at(client_assigned_to_facility)[0] << "," << clients.at(client_assigned_to_facility)[1] << endl;
}
output_file << endl;
}
output_file.close();
// output_file.close();
cout << "Number of facilities: " << num_facilities_opened << endl;
return 0;
}
// ---------------------------------------------------------
int main(int argc, char** argv) {
vector<array<double, 2> > data_points = parseCSV(argv[1]);
vector<array<double, 3> > distances = retrieveDistances(argv[1], data_points);
double lambda = atof(argv[2]);
// Run primal-dual algorithm
facility_location(data_points, distances, lambda, argv[1]);
/*
while (true){
cout << "Enter lambda: ";
cin >> lambda;
facility_location(data_points, distances, lambda, argv[1]);
}
*/
return 0;
}