IMGR.py

import VRPLibReader
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import TABU

distmat = VRPLibReader.distmat
dimension = VRPLibReader.n
depot = VRPLibReader.site[0]
sites = VRPLibReader.site
cap = VRPLibReader.capacity
visited = np.zeros((dimension, dimension))
flag = np.zeros((dimension, dimension))
supply = VRPLibReader.things
finished = np.zeros(dimension)
finished[0] = 1


def value(a, distmat):
    s = 0
    for i in range(len(a) - 1):
        s += distmat[a[i]][a[i + 1]]
    return s


def find(a, b):
    for i in range(a.shape[0]):
        if np.sum(a[i, :] == b) == b.size:
            return i


def createKDTree(dataSet, depth):
    n = np.shape(dataSet)[0]
    treeNode = {}
    if n == 0:
        return None
    else:
        n, m = np.shape(dataSet)
        split_axis = depth % m
        depth += 1
        treeNode['split'] = split_axis
        dataSet = np.array(sorted(dataSet, key=lambda a: a[split_axis]))
        num = n // 2
        treeNode['median'] = dataSet[num]
        treeNode['left'] = createKDTree(dataSet[:num], depth)
        treeNode['right'] = createKDTree(dataSet[num + 1:], depth)
        # treeNode['number'] = int(np.where((VRPLibReader.site == dataSet[num]).all(1))[0])
        treeNode['number'] = find(VRPLibReader.site, dataSet[num])
        return treeNode


def searchTree(tree, i):
    data = sites[i]
    k = len(data)
    if tree is None:
        return [0] * k, float('inf'), 0
    split_axis = tree['split']
    median_point = tree['median']
    num = tree['number']
    if data[split_axis] <= median_point[split_axis]:
        nearestPoint, nearestDistance, nearestnum = searchTree(tree['left'], i)
    else:
        nearestPoint, nearestDistance, nearestnum = searchTree(tree['right'], i)
    nowDistance = np.linalg.norm(data - median_point)
    if nearestDistance > nowDistance > 0 and flag[i][num] == 0 and flag[num][i] == 0 and finished[num] == 0:
        nearestnum = num
        nearestDistance = nowDistance
        nearestPoint = median_point.copy()
    splitDistance = abs(data[split_axis] - median_point[split_axis])
    if splitDistance > nearestDistance > 0:
        return nearestPoint, nearestDistance, nearestnum
    else:
        if data[split_axis] <= median_point[split_axis]:
            nextTree = tree['right']
        else:
            nextTree = tree['left']
        nearPoint, nearDistance, nearnum = searchTree(nextTree, i)
        if nearestDistance > nearDistance > 0 and flag[i][nearnum] == 0 and flag[nearnum][i] == 0 and finished[
            num] == 0:
            nearestnum = nearnum
            nearestDistance = nearDistance
            nearestPoint = nearPoint.copy()
        return nearestPoint, nearestDistance, nearestnum


class ZHeap:
    def __init__(self, item=[]):
        # 初始化。item为数组
        self.items = item
        self.heapsize = len(self.items)

    def LEFT(self, i):
        return 2 * i + 1

    def RIGHT(self, i):
        return 2 * i + 2

    def PARENT(self, i):
        return (i - 1) / 2

    def MIN_HEAPIFY(self, i):
        # 最小堆化：使以i为根的子树成为最小堆
        l = self.LEFT(i)
        r = self.RIGHT(i)
        if l < self.heapsize and self.items[l][2] < self.items[i][2]:
            smallest = l
        else:
            smallest = i

        if r < self.heapsize and self.items[r][2] < self.items[smallest][2]:
            smallest = r

        if smallest != i:
            self.items[i], self.items[smallest] = self.items[smallest], self.items[i]
            self.MIN_HEAPIFY(smallest)

    def INSERT(self, val):
        # 插入一个值val，并且调整使满足堆结构
        self.items.append(val)
        idx = len(self.items) - 1
        parIdx = int(self.PARENT(idx))
        while parIdx >= 0:
            if self.items[parIdx][2] > self.items[idx][2]:
                self.items[parIdx], self.items[idx] = self.items[idx], self.items[parIdx]
                idx = parIdx
                parIdx = int(self.PARENT(parIdx))
            else:
                break
        self.heapsize += 1

    def DELETE(self, i):
        last = len(self.items) - 1
        if last < 0:
            return None
        self.items[i], self.items[last] = self.items[last], self.items[i]
        val = self.items.pop()
        self.heapsize -= 1
        self.MIN_HEAPIFY(i)
        return val

    def BUILD_MIN_HEAP(self):
        # 建立最小堆, O(nlog(n))
        i = self.PARENT(len(self.items) - 1)
        while i >= 0:
            self.MIN_HEAPIFY(i)
            i -= 1

    def SHOW(self):
        print(self.items)


class ZPriorityQ(ZHeap):
    def __init__(self, item=[]):
        ZHeap.__init__(self, item)

    def insert(self, val):
        ZHeap.INSERT(self, val)

    def delete(self, i):
        val = ZHeap.DELETE(self, i)
        return val


def suit(edge):
    if finished[edge[1]] == 1 or finished[edge[0]] == 1:
        return False
    if visited[edge[0]].sum() < 2 and visited[edge[1]].sum() < 2:
        if visited[edge[0]].sum() + visited[edge[1]].sum() <= 1:
            return True
    return False


start_time = datetime.datetime.now()
tree_depth = int(math.log10(dimension) - 2)
KD_Tree = createKDTree(sites[1:], tree_depth)
PAI = []
for i in range(len(distmat)):
    PAI.append(0 * distmat[0][i]*(1-supply[i]/supply.sum()))
    distmat[i][i] = float('inf')
for i in range(len(distmat)):
    for j in range(len(distmat)):
        distmat[i][j] += (PAI[i] + PAI[j])
temp = []
for i in range(1, len(distmat)):
    a = searchTree(KD_Tree, i)[2]
    flag[i][a] = flag[a][i] = 1
    temp.append([i, a, distmat[i][a]])
Min_Priority = ZPriorityQ()
for i in temp:
    Min_Priority.insert(i)
paths = []
while len(Min_Priority.items) != 0:
    weight = 0
    path = []
    i = 0
    while weight < cap:
        if i >= len(Min_Priority.items):
            break
        edge = Min_Priority.items[i]
        if edge:
            if suit(edge):
                if weight + supply[edge[0]] * (1 - visited[edge[0]].sum()) + supply[edge[1]] * (
                        1 - visited[edge[1]].sum()) <= cap:
                    if visited[edge[0]].sum() == 0:
                        weight += supply[edge[0]]
                    if visited[edge[1]].sum() == 0:
                        weight += supply[edge[1]]
                    path.append(edge[:2])
                    visited[edge[0]][edge[1]] += 1
                    visited[edge[1]][edge[0]] += 1
                    Min_Priority.delete(i)
                    if visited[edge[0]].sum() < 2 and flag[edge[0]].sum() < dimension - 2:
                        a = searchTree(KD_Tree, edge[0])[2]
                        flag[edge[0]][a] = flag[a][edge[0]] = 1
                        Min_Priority.insert([edge[0], a, distmat[edge[0]][a]])
                else:
                    i += 1
            else:
                Min_Priority.delete(i)
                if visited[edge[0]].sum() < 2 and flag[edge[0]].sum() < dimension - 2:
                    a = searchTree(KD_Tree, edge[0])[2]
                    if a == 0:
                        break
                    flag[edge[0]][a] = flag[a][edge[0]] = 1
                    Min_Priority.insert([edge[0], a, distmat[edge[0]][a]])
    for i in path:
        visited[i[0]][0] = 1
        visited[i[1]][0] = 1
        finished[i[0]] = finished[i[1]] = 1
        flag[i[0]] = 1
    if path:
        paths.append(path)
c = 0
path = []
while finished.all != 1:
    while finished[c] == 1:
        c += 1
        if c==len(finished):
            break
    if c != len(finished):
        finished[c]=1
        path.append([c, 0])
    if c == len(finished):
        paths.append(path)
        break
end_time = datetime.datetime.now()


class Customer:

    def __init__(self, num, demand):
        self.num = num
        self.demand = demand
        self.is_visited = False


class UAV:

    def __init__(self, cap):
        self.cap = cap
        self.load = 0
        self.routes = []
        self.current_location = 0

    def check_if_fit(self, demand):
        return self.load + demand <= self.cap

    def transport(self, cus):
        self.routes.append(cus)
        self.load += cus.demand
        self.current_location = cus.num


realpath = []
temp = np.zeros(dimension)
for i in paths:
    tem = [0, ]
    for j in i:
        if not temp[j[0]]:
            tem.append(j[0])
            temp[j[0]] = 1
        if not temp[j[1]]:
            temp[j[1]] = 1
            tem.append(j[1])
    realpath.extend(tem)
print(realpath)
print('算法时间:', end_time - start_time)
MAXCOUNT = 1000
tem = []
temps = []
while realpath[0] == 0:
    del realpath[0]
while realpath[len(realpath) - 1] == 0:
    del realpath[len(realpath) - 1]
for z in realpath:
    if (z != 0):
        tem.append(z)
    else:
        a = []
        for j in tem:
            a.append(j)
        temps.append(a)
        tem.clear()
temps.append(tem)
for z in temps:
    z.insert(0, 0)
newtabu = []
for i in temps:
    city = []
    for z in i:
        city.append(VRPLibReader.site[z])
    cities = np.array(city)


    def calDist(xindex, yindex):
        return (np.sum(np.power(cities[xindex] - cities[yindex], 2))) ** 0.5


    def calPathDist(indexList):
        sum = 0.0
        for i in range(1, len(indexList)):
            sum += calDist(indexList[i], indexList[i - 1])
        return sum


    # path1长度比path2短则返回true
    def pathCompare(path1, path2):
        if calPathDist(path1) <= calPathDist(path2):
            return True
        return False


    def generateRandomPath(bestPath):
        a = np.random.randint(len(bestPath))
        while True:
            b = np.random.randint(len(bestPath))
            if np.abs(a - b) > 1:
                break
        if a > b:
            return b, a, bestPath[b:a + 1]
        else:
            return a, b, bestPath[a:b + 1]


    def reversePath(path):
        rePath = path.copy()
        rePath[1:-1] = rePath[-2:0:-1]
        return rePath


    def updateBestPath(bestPath):
        count = 0
        while count < MAXCOUNT and len(bestPath) > 3:
            start, end, path = generateRandomPath(bestPath)
            rePath = reversePath(path)

            if pathCompare(path, rePath):
                count += 1
                continue

            else:
                count = 0
                bestPath[start:end + 1] = rePath
        return bestPath


    def opt2():
        # 随便选择一条可行路径
        bestPath = np.arange(0, len(cities))
        bestPath = np.append(bestPath, 0)
        bestPath = updateBestPath(bestPath)
        return bestPath


    ttt = opt2().tolist()
    aaa = []
    for j in range(1, len(ttt)):
        aaa.append(i[ttt[j]])
    bbb = aaa * 2
    a = 0
    while bbb[a] != 0:
        a += 1
    newtabu.append(0)
    a += 1
    while bbb[a] != 0:
        newtabu.append(bbb[a])
        a += 1
        if a == len(bbb):
            break
newtabu.append(0)
newlen = value(newtabu, VRPLibReader.getdistmat(VRPLibReader.site))
print("newpath", newtabu)
print("newlen", newlen)
num = newtabu
total_cost = newlen
while num[0] == 0:
    del num[0]
while num[len(num) - 1] == 0:
    del num[len(num) - 1]
t1 = []
temp_arr = []
for i in num:
    if i != 0:
        t1.append(i)
    else:
        a = []
        for j in t1:
            a.append(j)
        temp_arr.append(a)
        t1.clear()
temp_arr.append(t1)
for i in temp_arr:
    i.append(0)
    i.insert(0, 0)
uavs = []
for i in range(len(temp_arr)):
    uavs.append(UAV(VRPLibReader.capacity))
for i in range(len(temp_arr)):
    for j in temp_arr[i]:
        aaa = Customer(j, VRPLibReader.things[j])
        uavs[i].transport(aaa)
TABU.tabu_solver(uavs,total_cost)