From 138696b4f7b209d185c550f6a0d26606e690679a Mon Sep 17 00:00:00 2001 From: rrooman69 Date: Mon, 4 Apr 2016 20:56:16 +0200 Subject: [PATCH] Create cvrp.py added calculation of distances from GEO coordinates and reading of EXPLICIT distances --- cvrp.py | 194 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 194 insertions(+) create mode 100644 cvrp.py diff --git a/cvrp.py b/cvrp.py new file mode 100644 index 0000000..fb1195b --- /dev/null +++ b/cvrp.py @@ -0,0 +1,194 @@ +from math import radians, cos, sin, asin, sqrt,pi,floor,atan2 +from math import acos + +import random +import copy +import union_find as uf +import itertools +from collections import deque +import pp + +RRR = 6378388 # r van aarde in m + +class ProblemInstance: + def __str__(self): + return """Name: %s +Type : %s +Comment: %s +Nodes: %s +Dimension: %s +Capacity: %s +Edge weight: %s +Node demand: %s +Node coordinate type: %s +""" % (self.name, + self.type, + self.comment, + self.node_coord, + self.dimension, + self.capacity, + self.node_distances, + self.node_demand, + self.ewt) + + def prepare(self): + if self.ewt == "EUC_2D": + def metric((x, y), (a, b)): + return sqrt((x - a) ** 2 + (y - b) ** 2) + if self.ewt == "GEO": + def metric((x, y), (a, b)): + x,y,a,b = map(radians,[x,y,a,b]) + q1= cos(y-b) + q2= cos(x-a) + q3= cos(x+a) + return int(RRR* acos(0.5*((1.0+q1)*q2-(1.0-q1)*q3))+1.0) + else: + def metric(a, b): + return 0 + if self.ewt != "EXPLICIT": + self.node_distances = [ [ metric(self.node_coord[i], self.node_coord[j]) for i in range(self.dimension) ] for j in range(self.dimension) ] + return self + + def __init__(self): + self.name = "" + self.type = "" + self.comment = "" + self.dimension = 0 + self.capacity = 0 + self.nct = "NO_COORDS" + self.ewt = "" + self.ewf = "" + self.edf = "" + self.ddt = "NO_DISPLAY" + + self.node_coord = {} + self.node_distances = {} + self.node_demand = {} + + +class Solution: + def __init__(self, p): + self.problem = p.name + self.capacity = p.capacity + print 'Als je', self.capacity,'posten kan dragen vinden we volgende routes:' + self.paths = [uf.Path(p.node_demand[i], i) for i in xrange(p.dimension)] + self.used_savings = deque() + self.solution = sum([2 * p.node_distances[0][i] for i in range(1, p.dimension)]) + + def feasibleSaving(self, saving): + (_, i, j, _) = saving + return uf.Joinable(self.capacity, self.paths[i], self.paths[j]) + + def processSaving(self, saving): + (k, i, j, s) = saving + x = self.paths[i] + y = self.paths[j] + if uf.Joinable(self.capacity, x, y): + uf.Union(x, y) + self.used_savings.append(k) + self.solution -= s + + +def randomSampling(solution, savings, current_best): + p = random.uniform(0.05, 0.2) + for s in savings: + if random.random() < p: + if solution.feasibleSaving(s): + solution.processSaving(s) + return solution + + +def saving(p, i, j): + c = p.node_distances + return c[0][i] + c[0][j] - c[i][j] + +def buildSavingsList(p): + savings = sorted([(i, j, saving(p, i, j)) for i in xrange(1, p.dimension) for j in xrange(i, p.dimension) if i != j], key=lambda tup: tup[2], reverse=True) + return [(i,) + savings[i] for i in xrange(len(savings))] + +def remote_job(chunk_size, solution, deque_savings, best): + local_best_cost = best + local_best_path = None + + local_solution = copy.deepcopy(solution) + probe = randomSampling(local_solution, deque_savings, local_best_cost) + local_sum = probe.solution + + for _ in xrange(1, chunk_size): + local_solution = copy.deepcopy(solution) + probe = randomSampling(local_solution, deque_savings, local_best_cost) + local_sum += probe.solution + if probe.solution < local_best_cost: + local_best_cost = probe.solution + local_best_path = probe.used_savings + return (local_best_cost, local_best_path, local_sum) + + +def MonteCarloSimulation(problem, ncpus,nprobes,servers): + ppservers = tuple(servers) + job_server = pp.Server(ppservers=ppservers) + + sol = Solution(problem) + best_cost = sol.solution + best_path = sol.used_savings + + savings = buildSavingsList(problem) + deq_savings = deque(savings) + if ncpus == None or ncpus < 1: + ncpus = job_server.get_ncpus() + chunk_size = nprobes / ncpus + print "There are %d cpus available." % ncpus + for _ in xrange(len(savings) - 10): + saving = deq_savings.popleft() + if sol.feasibleSaving(saving): + + unprocessed = sol + no = 0 + jobs = [job_server.submit(remote_job, (chunk_size, unprocessed, deq_savings, best_cost), (randomSampling,), ("random", "copy", "itertools")) for _ in xrange(ncpus)] + results = [job() for job in jobs] + for (lbest_cost, lbest_path, lsum) in results: + if lbest_cost < best_cost: + best_cost = lbest_cost + best_path = lbest_path + no += lsum + + processed = copy.deepcopy(sol) + processed.processSaving(saving) + yes = 0 + jobs = [job_server.submit(remote_job, (chunk_size, processed, deq_savings, best_cost), (randomSampling,), ("random", "copy", "itertools")) for _ in xrange(ncpus)] + results = [job() for job in jobs] + for (lbest_cost, lbest_path, lsum) in results: + if lbest_cost < best_cost: + best_cost = lbest_cost + best_path = lbest_path + yes += lsum + + if yes < no: + sol = processed + + l = savings[len(savings) - 10:] + subsets = list(itertools.chain(*[itertools.combinations(l, ni) for ni in xrange(1, len(l) + 1)])) + for subset in subsets: + solution_copy = copy.deepcopy(sol) + for s in subset: + if solution_copy.feasibleSaving(s): + solution_copy.processSaving(s) + if solution_copy.solution < best_cost: + best_cost = solution_copy.solution + best_path = solution_copy.used_savings + paths = [[i] for i in xrange(problem.dimension)] + for (i, j) in [savings[n][1:3] for n in best_path]: + if paths[i][0] == i and paths[j][0] == j: + newpath = paths[i][::-1] + paths[j] + elif paths[i][0] == i and paths[j][-1] == j: + newpath = paths[j] + paths[i] + elif paths[i][-1] == i and paths[j][0] == j: + newpath = paths[i] + paths[j] + elif paths[i][-1] == i and paths[j][-1] == j: + newpath = paths[i] + paths[j][::-1] + for k in newpath: + paths[k] = newpath + + job_server.print_stats() + return best_cost, [k for k, _ in itertools.groupby(sorted(paths[1:]))] +