from time import perf_counter           #Performance counter - clock with highest available resolution
from math import sqrt                   #Square root
from bisect import bisect_left          #Provides would-be index of element to insert
from collections import defaultdict     #Dictionary that does not throw KeyErrors
import numpy as np                      #http://www.numpy.org/
import random                           #RNG
import sys                              #Command-line arguments
import os                               #Directory check and creation

def DRTMultiAVRgen(M):

    #If no directory exists...
    if not os.path.exists('DRTMultiAVROutputs/'):
        #Create output directory
        os.makedirs('DRTMultiAVROutputs/')

    #Units
    #a_max, a_min : revolutions / min^2
    #speed, peak_speed, speed_new : revolutions / minute (RPM)
    #Acceleration equal in magnitude to deceleration per Bijinemula et al. Sec. III.A Para. 4 
    a_max = 600000
    a_min = -a_max

    if M>0:

        # Multiple AVR Tasks generation
        fileMname = 'DRTMultiAVROutputs/DRTAlg_MultiAVR_'+str(M)+'.txt'

        # start1 = perf_counter()
        random.seed(100)
        Ustar = 0.25
        Udict = dict()
        for i in range(M):
            Udict[i] = random.uniform(0.85*Ustar, Ustar)
        Udict[M-1] = Ustar

        while True:
            while True:
                boundarySpeeds = sorted(random.sample(range(500,6500),M-1))
                boundarySpeeds = [500]+boundarySpeeds+[6500]
                for i in range(M):
                    k = i
                    if boundarySpeeds[i+1] - boundarySpeeds[i] < 3000/M:
                        continue
                        
                if k==M-1:
                    break
            br = 0
            executionTimes = []
            for i in range(M):
                executionTimes.append(1000**2*Udict[i]*60/boundarySpeeds[i])   #in micro seconds
            for c1,c2 in zip(executionTimes[:-1],executionTimes[1:]):
                if c1<c2:
                    br=1
                    break
            if br==0:
                break

        # print(executionTimes)
        # print()
        # print(boundarySpeeds)
        # print()
        print('Modes = ',M)
        # end1 = perf_counter()
        # print('Time Taken to compute Random values is ', (end1-start1)*1000)

    else:
        #if M<0, run the algorithm on a predefined taskset
        import json
        taskSetFileName = 'taskSet'+str(abs(M))+'.json'
        fileMname = 'DRTMultiAVROutputs/DRTAlg_'+str(abs(M))+'.txt'
        
        with open(taskSetFileName) as f:
            taskset = json.load(f)

        boundarySpeeds = sorted(taskset['boundarySpeeds'])
        executionTimes = sorted(taskset['executionTimes'], reverse = True)

    #Start timer
    start = perf_counter()

    #Boundary Speeds
        #Squares of right boundary speeds used to avoid repeated sqrt computation later
        #Speeds in the first step are not counted per Lemma 2 - start from first right boundary speed
    boundarySpeeds = [speed**2 for speed in boundarySpeeds]
    maxSpeed = boundarySpeeds[-1]

    #Setup nodeSpeeds
    nodeSpeeds = set()

    #Iterate through all boundary speeds, creating nodes for all reachable speeds
    for i in range(len(boundarySpeeds)):
        speed = boundarySpeeds[i]           #Select boundary speed
        while speed <= maxSpeed:            #While speed has not exceed maxspeed
            nodeSpeeds.add(speed)               #Create node for speed
            speed = speed + 2*a_max             #Add subsequent, increasing speed via a_max
        
        speed = boundarySpeeds[i]           #Select boundary speed
        while speed>=boundarySpeeds[0]:     #While speed has not fallen below the lowest boundary speed
            nodeSpeeds.add(speed)               #Create node for speed
            speed = speed + 2*a_min             #Add subsequent, decreasing speed via a_min
    nodeSpeeds = sorted(nodeSpeeds)         #Nodes sorted by speed increasing order

    #Setup nodes
    nodes = []

    for i,j in zip(nodeSpeeds[:-1],nodeSpeeds[1:]): #For each tuple in a list of created tuples
        nodes.append((i,j))                             #Add each tuple to nodes

    adjMatrix = dict()  #Empty adjacency matrix for creating DRT graph

    for i in range(len(nodes)):                 #For every node tuple
        currentRB = nodes[i][1]                     #Select right boundary
        for j in range(i,len(nodes)):               #Iterate through all possible next right boundaries
            nextRB = nodes[j][1]                        #Select next right boundary
            maxReach = currentRB + 2*a_max              #Calculate maximum next speed

            #If next right boundary is not reachable, break
            if nextRB > maxReach:
                break

            #If next right boundary is reachable via constant a_max...
            if nextRB == maxReach: 
                #Calculate minimum interarrival time - Mohaqeqi et al. Sec. 3.2 Case 2
                adjMatrix[(i,j)] = 60*(sqrt(nextRB)-sqrt(currentRB))/a_max 

            #...otherwise, the next right boundary is reachable via variable acceleration
            else:
                #Calculate the prospective peak speed
                midSpeed = (2*a_min*a_max+a_min*currentRB-a_max*nextRB)/(a_min-a_max)

                #If the prospective peak speed does not exceed maximum speed...
                if midSpeed <= maxSpeed:
                    #Calculate3 minimum interarrival time in seconds - Mohaqeqi et al. Sec 3.2 Eqn 20
                    adjMatrix[(i,j)] = 60*((sqrt(midSpeed)-sqrt(currentRB))/a_max +(sqrt(nextRB)-sqrt(midSpeed))/a_min)

                #...otherwise, the prospective peak speed exceeds maximum speed
                else:
                    #Calculate the time to reach maximum speed - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_1^*
                    t1 = (sqrt(maxSpeed) - sqrt(currentRB))/a_max
                    #Calculate the time over which maximum speed is maintained - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_2^*
                    t2 = (1-((maxSpeed-currentRB)/(2*a_max))-((nextRB-maxSpeed)/(2*a_min)))/sqrt(maxSpeed)
                    #Calculate the time to descend from maximum speed to final speed - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_3^*
                    t3 = (sqrt(nextRB)-sqrt(maxSpeed))/a_min

                    #Sum individual time segments and conver to seconds
                    adjMatrix[(i,j)] = 60*(t1 + t2 + t3)

            #Mirror adjacency matrix
            adjMatrix[(j,i)] = adjMatrix[(i,j)]
            

    #Dictionary for logging speeds, completion times - Supports dynamic programming.
    hashTable = defaultdict(dict)

    #Function for calculating maximum demand given a starting node over a set duration of time
    def calc_demand(i,time):
        
        #Initialization
        demand_max = 0
        time = round(time,5)
        
        #Stored demand check
        if i in hashTable.keys():
            if time in hashTable[i].keys():
                return hashTable[i][time]
        
        #If node is not the first right bounday speed...
        if nodes[i][1]!=boundarySpeeds[0]:
            #Assign execution time based on index
            demNode = executionTimes[bisect_left(boundarySpeeds,nodes[i][1])-1]
        #...otherwise, assign first execution time
        else:
            demNode = executionTimes[0]
        
        #Iterate through all nodes
        for j in range(len(nodes)):

            #If the node pair selected is reachable from each other...
            if (i,j) in adjMatrix.keys():

                #Extract the time required to reach the next node
                timeNextNode = adjMatrix[(i,j)]
                
                #If insufficient time remains to reach the next node...
                if time - timeNextNode < 0:
                    #Do not count towards demand, try next node
                    continue
                
                #If time remaining is exactly the required time...
                if time - timeNextNode == 0:
                    #Update demand if necessary and continue
                    if demNode > demand_max:
                        demand_max = demNode
                        continue
                
                #Calculate new demand
                demand = demNode + calc_demand(j,time-timeNextNode)
            
                #Update demand if necessary
                if demand > demand_max:
                    demand_max = demand
        
        #Update hash table with max demand and return
        hashTable[i][time] = demand_max
        return demand_max

    #f = open('DRTAlgoOutput.txt','w')
    for tot_time in np.arange(0.01,1.01,0.01):
        max_demand = 0
        for i in range(len(nodes)):
            demand = calc_demand(i,tot_time)
            if demand > max_demand:
                max_demand = demand
        #print('Max demand for time: {} = {}'.format(tot_time,max_demand))
        #f.write('Max demand for time: {} = {}\n'.format(tot_time,max_demand))
    #f.close()

    #End performance counting, calculate and log total time
    end = perf_counter()
    fileM = open(fileMname,'a')
    total_time = end - start
    fileM.write(str((end-start))+'\n')
    #print('Total time is ',total_time)
    fileM.close()

if __name__ == '__main__':
    M = int(sys.argv[1])
    DRTMultiAVRgen(M)
from time import perf_counter           #Performance counter - clock with highest available resolution
from math import sqrt                   #Square root
from bisect import bisect_left          #Provides would-be index of element to insert
from collections import defaultdict     #Dictionary that does not throw KeyErrors
import numpy as np                      #http://www.numpy.org/
import random                           #RNG
import sys                              #Command-line arguments
import os                               #Directory check and creation

def DRTMultiAVRgen(M):

    #If no directory exists...
    if not os.path.exists('DRTMultiAVROutputs/'):
        #Create output directory
        os.makedirs('DRTMultiAVROutputs/')

    #Units
    #a_max, a_min : revolutions / min^2
    #speed, peak_speed, speed_new : revolutions / minute (RPM)
    #Acceleration equal in magnitude to deceleration per Bijinemula et al. Sec. III.A Para. 4 
    a_max = 600000
    a_min = -a_max

    if M>0:

        # Multiple AVR Tasks generation
        fileMname = 'DRTMultiAVROutputs/DRTAlg_MultiAVR_'+str(M)+'.txt'

        # start1 = perf_counter()
        random.seed(100)
        Ustar = 0.25
        Udict = dict()
        for i in range(M):
            Udict[i] = random.uniform(0.85*Ustar, Ustar)
        Udict[M-1] = Ustar

        while True:
            while True:
                boundarySpeeds = sorted(random.sample(range(500,6500),M-1))
                boundarySpeeds = [500]+boundarySpeeds+[6500]
                for i in range(M):
                    k = i
                    if boundarySpeeds[i+1] - boundarySpeeds[i] < 3000/M:
                        continue
                        
                if k==M-1:
                    break
            br = 0
            executionTimes = []
            for i in range(M):
                executionTimes.append(1000**2*Udict[i]*60/boundarySpeeds[i])   #in micro seconds
            for c1,c2 in zip(executionTimes[:-1],executionTimes[1:]):
                if c1<c2:
                    br=1
                    break
            if br==0:
                break

        # print(executionTimes)
        # print()
        # print(boundarySpeeds)
        # print()
        print('Modes = ',M)
        # end1 = perf_counter()
        # print('Time Taken to compute Random values is ', (end1-start1)*1000)

    else:
        #if M<0, run the algorithm on a predefined taskset
        import json
        taskSetFileName = 'taskSet'+str(abs(M))+'.json'
        fileMname = 'DRTMultiAVROutputs/DRTAlg_'+str(abs(M))+'.txt'
        
        with open(taskSetFileName) as f:
            taskset = json.load(f)

        boundarySpeeds = sorted(taskset['boundarySpeeds'])
        executionTimes = sorted(taskset['executionTimes'], reverse = True)

    #Start timer
    start = perf_counter()

    #Boundary Speeds
        #Squares of right boundary speeds used to avoid repeated sqrt computation later
        #Speeds in the first step are not counted per Lemma 2 - start from first right boundary speed
    boundarySpeeds = [speed**2 for speed in boundarySpeeds]
    maxSpeed = boundarySpeeds[-1]

    #Setup nodeSpeeds
    nodeSpeeds = set()

    #Iterate through all boundary speeds, creating nodes for all reachable speeds
    for i in range(len(boundarySpeeds)):
        speed = boundarySpeeds[i]           #Select boundary speed
        while speed <= maxSpeed:            #While speed has not exceed maxspeed
            nodeSpeeds.add(speed)               #Create node for speed
            speed = speed + 2*a_max             #Add subsequent, increasing speed via a_max
        
        speed = boundarySpeeds[i]           #Select boundary speed
        while speed>=boundarySpeeds[0]:     #While speed has not fallen below the lowest boundary speed
            nodeSpeeds.add(speed)               #Create node for speed
            speed = speed + 2*a_min             #Add subsequent, decreasing speed via a_min
    nodeSpeeds = sorted(nodeSpeeds)         #Nodes sorted by speed increasing order

    #Setup nodes
    nodes = []

    for i,j in zip(nodeSpeeds[:-1],nodeSpeeds[1:]): #For each tuple in a list of created tuples
        nodes.append((i,j))                             #Add each tuple to nodes

    adjMatrix = dict()  #Empty adjacency matrix for creating DRT graph

    for i in range(len(nodes)):                 #For every node tuple
        currentRB = nodes[i][1]                     #Select right boundary
        for j in range(i,len(nodes)):               #Iterate through all possible next right boundaries
            nextRB = nodes[j][1]                        #Select next right boundary
            maxReach = currentRB + 2*a_max              #Calculate maximum next speed

            #If next right boundary is not reachable, break
            if nextRB > maxReach:
                break

            #If next right boundary is reachable via constant a_max...
            if nextRB == maxReach: 
                #Calculate minimum interarrival time - Mohaqeqi et al. Sec. 3.2 Case 2
                adjMatrix[(i,j)] = 60*(sqrt(nextRB)-sqrt(currentRB))/a_max 

            #...otherwise, the next right boundary is reachable via variable acceleration
            else:
                #Calculate the prospective peak speed
                midSpeed = (2*a_min*a_max+a_min*currentRB-a_max*nextRB)/(a_min-a_max)

                #If the prospective peak speed does not exceed maximum speed...
                if midSpeed <= maxSpeed:
                    #Calculate3 minimum interarrival time in seconds - Mohaqeqi et al. Sec 3.2 Eqn 20
                    adjMatrix[(i,j)] = 60*((sqrt(midSpeed)-sqrt(currentRB))/a_max +(sqrt(nextRB)-sqrt(midSpeed))/a_min)

                #...otherwise, the prospective peak speed exceeds maximum speed
                else:
                    #Calculate the time to reach maximum speed - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_1^*
                    t1 = (sqrt(maxSpeed) - sqrt(currentRB))/a_max
                    #Calculate the time over which maximum speed is maintained - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_2^*
                    t2 = (1-((maxSpeed-currentRB)/(2*a_max))-((nextRB-maxSpeed)/(2*a_min)))/sqrt(maxSpeed)
                    #Calculate the time to descend from maximum speed to final speed - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_3^*
                    t3 = (sqrt(nextRB)-sqrt(maxSpeed))/a_min

                    #Sum individual time segments and conver to seconds
                    adjMatrix[(i,j)] = 60*(t1 + t2 + t3)

            #Mirror adjacency matrix
            adjMatrix[(j,i)] = adjMatrix[(i,j)]
            

    #Dictionary for logging speeds, completion times - Supports dynamic programming.
    hashTable = defaultdict(dict)

    #Function for calculating maximum demand given a starting node over a set duration of time
    def calc_demand(i,time):
        
        #Initialization
        demand_max = 0
        time = round(time,5)
        
        #Stored demand check
        if i in hashTable.keys():
            if time in hashTable[i].keys():
                return hashTable[i][time]
        
        #If node is not the first right bounday speed...
        if nodes[i][1]!=boundarySpeeds[0]:
            #Assign execution time based on index
            demNode = executionTimes[bisect_left(boundarySpeeds,nodes[i][1])-1]
        #...otherwise, assign first execution time
        else:
            demNode = executionTimes[0]
        
        #Iterate through all nodes
        for j in range(len(nodes)):

            #If the node pair selected is reachable from each other...
            if (i,j) in adjMatrix.keys():

                #Extract the time required to reach the next node
                timeNextNode = adjMatrix[(i,j)]
                
                #If insufficient time remains to reach the next node...
                if time - timeNextNode < 0:
                    #Do not count towards demand, try next node
                    continue
                
                #If time remaining is exactly the required time...
                if time - timeNextNode == 0:
                    #Update demand if necessary and continue
                    if demNode > demand_max:
                        demand_max = demNode
                        continue
                
                #Calculate new demand
                demand = demNode + calc_demand(j,time-timeNextNode)
            
                #Update demand if necessary
                if demand > demand_max:
                    demand_max = demand
        
        #Update hash table with max demand and return
        hashTable[i][time] = demand_max
        return demand_max

    #f = open('DRTAlgoOutput.txt','w')
    for tot_time in np.arange(0.01,1.01,0.01):
        max_demand = 0
        for i in range(len(nodes)):
            demand = calc_demand(i,tot_time)
            if demand > max_demand:
                max_demand = demand
        #print('Max demand for time: {} = {}'.format(tot_time,max_demand))
        #f.write('Max demand for time: {} = {}\n'.format(tot_time,max_demand))
    #f.close()

    #End performance counting, calculate and log total time
    end = perf_counter()
    fileM = open(fileMname,'a')
    total_time = end - start
    fileM.write(str((end-start))+'\n')
    #print('Total time is ',total_time)
    fileM.close()

if __name__ == '__main__':
    M = int(sys.argv[1])
    DRTMultiAVRgen(M)
#Performance counter - clock with highest available resolution #Square root #Provides would-be index of element to insert #Dictionary that does not throw KeyErrors #http://www.numpy.org/ #RNG #Command-line arguments #Directory check and creation #If no directory exists... #Create output directory #Units #a_max, a_min : revolutions / min^2 #speed, peak_speed, speed_new : revolutions / minute (RPM) #Acceleration equal in magnitude to deceleration per Bijinemula et al. Sec. III.A Para. 4 # Multiple AVR Tasks generation # start1 = perf_counter() #in micro seconds # print(executionTimes) # print() # print(boundarySpeeds) # print() # end1 = perf_counter() # print('Time Taken to compute Random values is ', (end1-start1)*1000) #if M<0, run the algorithm on a predefined taskset #Start timer #Boundary Speeds #Squares of right boundary speeds used to avoid repeated sqrt computation later #Speeds in the first step are not counted per Lemma 2 - start from first right boundary speed #Setup nodeSpeeds #Iterate through all boundary speeds, creating nodes for all reachable speeds #Select boundary speed #While speed has not exceed maxspeed #Create node for speed #Add subsequent, increasing speed via a_max #Select boundary speed #While speed has not fallen below the lowest boundary speed #Create node for speed #Add subsequent, decreasing speed via a_min #Nodes sorted by speed increasing order #Setup nodes #For each tuple in a list of created tuples #Add each tuple to nodes #Empty adjacency matrix for creating DRT graph #For every node tuple #Select right boundary #Iterate through all possible next right boundaries #Select next right boundary #Calculate maximum next speed #If next right boundary is not reachable, break #If next right boundary is reachable via constant a_max... #Calculate minimum interarrival time - Mohaqeqi et al. Sec. 3.2 Case 2 #...otherwise, the next right boundary is reachable via variable acceleration #Calculate the prospective peak speed #If the prospective peak speed does not exceed maximum speed... #Calculate3 minimum interarrival time in seconds - Mohaqeqi et al. Sec 3.2 Eqn 20 #...otherwise, the prospective peak speed exceeds maximum speed #Calculate the time to reach maximum speed - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_1^* #Calculate the time over which maximum speed is maintained - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_2^* #Calculate the time to descend from maximum speed to final speed - Mohaqeqi et al. Sec. 3.2 Eqn 21 t_3^* #Sum individual time segments and conver to seconds #Mirror adjacency matrix #Dictionary for logging speeds, completion times - Supports dynamic programming. #Function for calculating maximum demand given a starting node over a set duration of time #Initialization #Stored demand check #If node is not the first right bounday speed... #Assign execution time based on index #...otherwise, assign first execution time #Iterate through all nodes #If the node pair selected is reachable from each other... #Extract the time required to reach the next node #If insufficient time remains to reach the next node... #Do not count towards demand, try next node #If time remaining is exactly the required time... #Update demand if necessary and continue #Calculate new demand #Update demand if necessary #Update hash table with max demand and return #f = open('DRTAlgoOutput.txt','w') #print('Max demand for time: {} = {}'.format(tot_time,max_demand)) #f.write('Max demand for time: {} = {}\n'.format(tot_time,max_demand)) #f.close() #End performance counting, calculate and log total time #print('Total time is ',total_time)