aco3.py

import random
import math
import numpy as np

seed = 1
boost = 20
maxIter = 50
numCities = 15
mu, sigma = 500, 200
demand = np.random.normal(mu,sigma,numCities)
demand[0] = 0
maxDistance = 40000
RP = 0.8 # Restriction percentage set as 80%
upperX = 24000
upperY = 32000
T = 8
t_u = 25/60
t_l = 25/60
upperDemand = 1200
Q = 2000
numBH = numCities

def randomMatrix(n, upperBound, seed):
    random.seed(seed)
    cities = []
    for r in range(n):
        sm = []
        cities.append(sm)
        for c in range(n):
            sm.append(upperBound*random.random())
    return cities

def sampleMatrix():
    return [[0,10,15,100],[100,0,100,5],[15,100,0,5],[100,5,5,0]]

def dis(a,b):
    return math.sqrt(math.pow((a[0]-b[0]),2) + math.pow((a[1]-b[1]),2))

def cityCoord(n,seed):
    random.seed(seed)
    cityC = []
    for r in range(n):
        coordinate = []
        coordinate.append(upperX*random.random())
        coordinate.append(upperY*random.random())
        cityC.append(coordinate)
    return cityC

def cityMatrix(n, cityC):
    cities = [[None for _ in range(n)] for _ in range(n)]
    for r in range(n):
        for c in range(n):
            cities[r][c] = dis(cityC[r],cityC[c])
    return cities

def Demand(n,upperDemand):
    demand = []
    for _ in range(n):
        demand.append(upperDemand*random.random())
        demand[0] = 0
    return demand

def pathLength(cities,path):
    return sum([cities[r][c] for (r,c) in zip(path,path[1:]+[path[0]])])

def updatePher(pher,path,boost):
    for (r,c) in zip(path,path[1:]+[path[0]]):
        pher[r][c] = pher[r][c] + boost

def evaporatePher(pher,maxIter,boost):
    decr = boost/float(maxIter)
    for r in range(len(pher)):
        for c in range(len(pher[r])):
            if pher[r][c] > decr:
                pher[r][c] = pher[r][c] - decr
            else:
                pher[r][c] = 0.0

def doSumWeight(cities, pher, used, current):
    runningTotal = 0.0
    for city in range(len(cities)):
        if not used.has_key(city):
            runningTotal = (runningTotal +
                            cities[current][city]*(1.0+pher[current][city]))
    return runningTotal

def findSumWeight(cities, pher, used, current, soughtTotal):
    runningTotal = 0.0
    next = 0
    for city in range(len(cities)):
        if runningTotal >= soughtTotal:
            break
        if not used.has_key(city):
            runningTotal = (runningTotal + 
                            cities[current][city]*(1.0+pher[current][city]))
            next = city
    return next

def genPath(cities, pher):
    #current = random.randint(0,len(cities)-1)
    current = 0
    path = [current]
    used = {current:1}
    while len(used) < len(cities):
        sumWeight = doSumWeight(cities,pher,used,current)
        rndValue = random.random()*sumWeight
        current = findSumWeight(cities,pher,used,current,rndValue)
        path.append(current)
        used[current] = 1
    return path

def bestPath(cities,seed,boost):
    pher = randomMatrix(len(cities),0,0)
    random.seed(seed)
    bestLen = numCities*maxDistance
    bestPath = []
    for iter in range(maxIter):
        #print (iter,bestLen)
        path = genPath(cities,pher)
        pathLen = pathLength(cities,path)
        if pathLen < bestLen:
            updatePher(pher,path,boost)
            bestLen = pathLen
            bestPath = path
        evaporatePher(pher,maxIter,boost)
    return bestPath

def TSPtoVRP(demand,path):
    sumDemand = 0
    sumTime = 0
    count = 0
    while count < len(path):
        sumDemand += demand[path[count]]
        if sumDemand > Q or sumTime > T:
            sumDemand = 0
            sumTime = 0
            path.insert(count,0)
        elif sumDemand == Q:
            sumDemand = 0
            sumTime = 0
            path.insert(count,0)
        #print "count = %s, len(path) = %s, sumDemand = %s, Q = %s"%(count,len(path),sumDemand,Q)
        sumTime += t_u
        count += 1
    return path

def main():
    print "starting"
    #cities = randomMatrix(numCities,maxDistance,seed)
    cityC = cityCoord(numCities,seed)
    BHC = cityCoord(numBH,seed)
    cities = cityMatrix(numCities,cityC)
    path = bestPath(cities,seed,boost)
    print path
    print "len = ",pathLength(cities,path)
  #  demand = Demand(numCities, upperDemand)
    path = TSPtoVRP(demand,path)
    print path
    countZero = 0
    for i in path:
        if i == 0:
            countZero += 1
    #print demand
    for i in range(len(path)):
        print (path[i],demand[path[i]])
    iter = 0
    lastZero = 0
    currentZero = 1
    fullPath = []
    for i in range(len(path)):
        p = []
        p.append(0)
        p.append(path[i])
        fullPath.append(p)
    print "cityC"
    print "len = ",len(cityC)
    while iter < countZero:
        while currentZero < len(path) and path[currentZero] != 0:
            currentZero += 1
        L_v = 0
        for i in range(lastZero,currentZero):
            L_v += demand[path[i]]
        i = lastZero
        sumDem = 0
        #print "within while loop for RP"
        while (Q-L_v+sumDem) < RP*Q:
            sumDem += demand[path[i]]
            #print "i = %s, demand[i] = %s, L_v = %s, Q-L_v+sumDem = %s"%(i,demand[path[i]],L_v,Q-L_v+sumDem)
            i += 1
        #print ("iter = %s, L_v = %s, i = %s, path[i] = %s, lastZero = %s, currentZero = %s"%(iter,L_v,i-1,path[i-1],lastZero,currentZero))
        print (iter,i-2)
        iter += 1
        centroid = []
        centroid.append(0)
        centroid.append(0)
        for j in range(i-2,currentZero):
            centroid[0] += cityC[path[j]][0]
            centroid[1] += cityC[path[j]][1]
        print centroid[0]/(currentZero-i+2),centroid[1]/(currentZero-i+2)
        lastZero = currentZero
        currentZero += 1
        
        minI = 0
        for i in range(numBH):
            if dis(BHC[i],centroid) < dis(BHC[minI],centroid):
                minI = i    

        a = []
        a.append(1)
        a.append(i)
        fullPath.insert(random.randint(0,len(cityC)-1),a)
    
    print fullPath
'''
        mI = i
        for j in range(i-1,currentZero):
            for k in range(j,currentZero):
                if dis(cityC[path[j]],cityC[path[k]]) < dis(cityC[path[mI]],cityC[path[k]]):
                    a = []
                    a.append(1)
                    a.append(mI)
                    fullPath.insert(j,a)
'''

        

if __name__ == "__main__":
    main()