# Program Name: NSGA-II.py
# Description: This is a python implementation of Prof. Kalyanmoy Deb's popular NSGA-II algorithm
# Author: Haris Ali Khan
# Supervisor: Prof. Manoj Kumar Tiwari

#Importing required modules
import math
import random
import matplotlib.pyplot as plt
import numpy as np

import Users
import tool
import GDOP
import time



#计算用户平均snr
def function_snr(users,uavs_pos,c_association,power):
    sum=0
    for user_index in range(len(c_association)):
        snr=tool.cal_snr(uavs_pos[c_association[user_index]],users[user_index],power[c_association[user_index]][user_index])
        sum+=snr
    return (sum/len(users))

#计算用户平均gdop
def function_gdop(users,uavs_pos,c_association,p_association):
    # p_association=tool.get_p_association(users,uavs_pos,c_association,len(users)*2)
    sum=0
    for user_index in p_association:
        uavs=[]
        uavs.append(uavs_pos[c_association[user_index]])
        for uav_index in p_association[user_index]:
            uavs.append(uavs_pos[uav_index])
        gdop=GDOP.calculate_gdop(uavs,users[user_index])
        # print(gdop)
        sum+=gdop
    return -sum/len(users)

def punish(power,p_rest):
    f = 0
    i=0
    for p in power:
        dis = sum(p) - p_rest[i]
        if dis > 0.001:
            f += 10 ** (dis * 100)
        for p_singe in p:  # 功率要大于0
            dis1 = 0 - p_singe
            if dis1>0:
                f += 10 ** max(0, dis1 * 100)
        i+=1
    return f

#Function to find index of list
def index_of(a,list):
    for i in range(0,len(list)):
        if list[i] == a:
            return i
    return -1

#Function to sort by values
def sort_by_values(list1, values):
    sorted_list = []
    while(len(sorted_list)!=len(list1)):
        if index_of(min(values),values) in list1:
            sorted_list.append(index_of(min(values),values))
        values[index_of(min(values),values)] = math.inf
    return sorted_list

#Function to carry out NSGA-II's fast non dominated sort
def fast_non_dominated_sort(values1, values2):
    #支配的节点
    S=[[] for i in range(0,len(values1))]
    front = [[]]
    n=[0 for i in range(0,len(values1))]#被支配的数量
    rank = [0 for i in range(0, len(values1))]#等级

    for p in range(0,len(values1)):
        S[p]=[]
        n[p]=0
        for q in range(0, len(values1)):
            if (values1[p] > values1[q] and values2[p] > values2[q]) or (values1[p] >= values1[q] and values2[p] > values2[q]) or (values1[p] > values1[q] and values2[p] >= values2[q]):
                if q not in S[p]:
                    S[p].append(q)
            elif (values1[q] > values1[p] and values2[q] > values2[p]) or (values1[q] >= values1[p] and values2[q] > values2[p]) or (values1[q] > values1[p] and values2[q] >= values2[p]):
                n[p] = n[p] + 1
        if n[p]==0:
            rank[p] = 0
            if p not in front[0]:
                front[0].append(p)

    i = 0
    while(front[i] != []):
        Q=[]
        for p in front[i]:
            for q in S[p]:
                n[q] =n[q] - 1
                if( n[q]==0):
                    rank[q]=i+1
                    if q not in Q:
                        Q.append(q)
        i = i+1
        front.append(Q)

    del front[len(front)-1]
    return front,rank

#Function to calculate crowding distance
def crowding_distance(values1, values2, front):
    distance = [0 for i in range(0,len(front))]
    #从小到大排序
    sorted1 = sort_by_values(front, values1[:])
    sorted2 = sort_by_values(front, values2[:])
    distance[0] = 4444444444444444
    distance[len(front) - 1] = 4444444444444444
    for k in range(1,len(front)-1):
        distance[k] = distance[k]+ (values1[sorted1[k+1]] - values1[sorted1[k-1]])/(max(values1)-min(values1))
    for k in range(1,len(front)-1):
        distance[k] = distance[k]+ (values2[sorted2[k+1]] - values2[sorted2[k-1]])/(max(values2)-min(values2))
    return distance


#初始化无人机位置
def init(pop_size,k,users):
    min_pos=0
    max_pos=1000
    pos_solution=[]
    p_solution=[]
    for i in range(pop_size):
        uavs_pos=[[min_pos+(max_pos-min_pos)*random.random(),min_pos+(max_pos-min_pos)*random.random(),100] for i in range(k)]
        p_alloc=init_p(uavs_pos,users)
        pos_solution.append(uavs_pos)
        p_solution.append(p_alloc)
    # print(pos_solution)
    return pos_solution,p_solution

def get_dis(i,uavs_pos,users,uav_to_user):
    dis=[]
    for k in uav_to_user[i]:
        dis.append(tool.calculate_3d_distance(uavs_pos[i],users[k]))
    return dis

def init_p(uavs_pos,users):
    c_association=tool.get_c_association(users,uavs_pos,math.ceil(len(users)/len(uavs_pos)))
    uav_to_user=tool.get_uav_to_user(c_association,uavs_pos)
    # print(uav_to_user)
    p_association = tool.get_p_association(users, uavs_pos, c_association, len(users) * 2)
    p_rest=tool.get_rest_p(users,uavs_pos,p_association)
    # print(p_rest)
    p_alloc=[]
    for i in range(len(uavs_pos)):
        #归一化生成 功率
        dis=get_dis(i,uavs_pos,users,uav_to_user)
        # p=[random.randint(1,len(uav_to_user[i])) for i in range(len(uav_to_user[i]))]
        p_new=[0 for i in range(len(users))]
        for j in range(len(uav_to_user[i])):
            p_new[uav_to_user[i][j]]=dis[j]/sum(dis)*p_rest[i]
            # p_new.append(p[j]/sum(p)*p_rest[i])
        p_alloc.append(p_new)
        # p_alloc.append(np.random.uniform(0, p_rest[i] / len(uav_to_user[i]), (len(uav_to_user[i]))))
    return p_alloc

def cross_over_old(pos_solution1,pos_solution2):
    child1 = [pos_solution2[j] if random.random() < 0.5 else pos_solution1[j]
              for j in range(len(pos_solution2))]
    return child1

def cross_over_stage(pos_solution1,pos_solution2,iter,iter_max):
    p_cross=get_beta_stage(iter,iter_max)
    child1 = [pos_solution2[j] if random.random() < p_cross else pos_solution1[j]
              for j in range(len(pos_solution2))]
    return child1

#交叉 pos_solution1=[[],[],[]]
def cross_over(pos_solution1,pos_solution2,iter,iter_max):
    p_cross=get_beta(iter,iter_max)
    child1 = [pos_solution2[j] if random.random() < p_cross else pos_solution1[j]
              for j in range(len(pos_solution2))]
    # for i in range(len(pos_solution1)):
    #     pos1 = []
    #     pos2=[]
    #     beta = get_beta()
    #
        # for j in range(len(pos_solution1[i])):
        #     if j<len(pos_solution1[i])-1:
        #         x1 = (pos_solution1[i][j] + pos_solution2[i][j]) / 2
        #         x2 = abs((pos_solution1[i][j] - pos_solution2[i][j]) / 2)
        #         x1=x1 + beta * x2
        #         x2=abs(x1 - beta * x2)
        #         pos1.append(x1)
        #         pos2.append(x2)
        #     else:
        #         pos1.append(100)
        #         pos2.append(100)
        # child1.append(pos1)
        # child2.append(pos2)
    return child1

def cross_over2(pos_solution1,pos_solution2,iter,iter_max):
    p_cross=get_beta2(iter,iter_max)
    child1 = [pos_solution2[j] if random.random() < p_cross else pos_solution1[j]
              for j in range(len(pos_solution2))]
    return child1

def get_beta_old(crossover_param=2):
    u = random.random()
    if u <= 0.5:
        return (2 * u) ** (1 / (crossover_param + 1))
    return (2 * (1 - u)) ** (-1 / (crossover_param + 1))
def get_beta_stage(iter,iter_max):
    if iter/iter_max<=0.34:
        return 0.9
    elif iter/iter_max>0.34 and iter/iter_max<=0.67:
        return 0.4
    else:
        return 0.1
def get_beta(iter,iter_max,p_max=0.9,p_min=0.1,crossover_param=2):
    return p_max-(p_max-p_min)*iter/iter_max
    # u = random.random()
    # if u <= 0.5:
    #     return (2 * u) ** (1 / (crossover_param + 1))
    # return (2 * (1 - u)) ** (-1 / (crossover_param + 1))
def get_beta2(iter,iter_max,p_max=0.9,p_min=0.1,crossover_param=2):
    return p_min+(p_max-p_min)/2*(1+math.cos(iter/iter_max*math.pi))

def get_delta(iter,iter_max,p_max=0.1,p_min=0.01,mutation_param=5):
    u = random.random()
    a=random.random()
    delta=p_max-(p_max-p_min)*iter/iter_max
    if a<0.5:
        return u,delta
    return u,-delta
    # if u < 0.5:
    #     return u, (2 * u) ** (1 / (mutation_param + 1)) - 1
    # return u, 1 - (2 * (1 - u)) ** (1 / (mutation_param + 1))

def get_delta2(iter,iter_max,p_max=0.1,p_min=0.01,mutation_param=5):
    u = random.random()
    a=random.random()
    delta=p_min+(p_max-p_min)/2*(1+math.cos(iter/iter_max*math.pi))
    if a<0.5:
        return u,delta
    return u,-delta

def get_delta_stage(iter,iter_max):
    u = random.random()
    a=random.random()
    if iter / iter_max <= 0.34:
        delta = 0.1
    elif iter / iter_max > 0.34 and iter / iter_max <= 0.67:
        delta = 0.05
    else:
        delta = 0.01
    if a<0.5:
        return u,delta
    return u,-delta
    # if u < 0.5:
    #     return u, (2 * u) ** (1 / (mutation_param + 1)) - 1
    # return u, 1 - (2 * (1 - u)) ** (1 / (mutation_param + 1))
def get_delta_old(mutation_param=5):
    u = random.random()
    return u,0.05
    # if u < 0.5:
    #     return u, (2 * u) ** (1 / (mutation_param + 1)) - 1
    # return u, 1 - (2 * (1 - u)) ** (1 / (mutation_param + 1))

def mutate_old(child,variables_range):
    u, delta = get_delta_old()
    child_new = []
    if u < 0.5:
        for i in range(len(child)):
            pos = []
            for j in range(len(child[i])):
                if j < len(child[i]) - 1:
                    x = child[i][j] + delta * (child[i][j] - variables_range[0])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    else:
        for i in range(len(child)):
            pos = []
            for j in range(len(child[i])):
                if j < len(child[i]) - 1:
                    x = child[i][j] + delta * (variables_range[1] - child[i][j])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    for i in range(len(child_new)):
        for j in range(len(child_new[i])):
            if j < len(child_new[i]) - 1:
                if child_new[i][j] < variables_range[0]:
                    child_new[i][j] = variables_range[0]
                elif child_new[i][j] > variables_range[1]:
                    child_new[i][j] = variables_range[1]
    return child_new
def mutate_stage(child,variables_range,iter,iter_max):
    u, delta = get_delta_stage(iter,iter_max)
    child_new=[]
    if u < 0.5:
        for i in range(len(child)):
            pos=[]
            for j in range(len(child[i])):
                if j < len(child[i]) - 1:
                    x=child[i][j]+ delta * (child[i][j] - variables_range[0])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    else:
        for i in range(len(child)):
            pos = []
            for j in range(len(child[i])):
                if j<len(child[i])-1:
                    x= child[i][j] +delta * (variables_range[1] - child[i][j])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    for i in range(len(child_new)):
        for j in range(len(child_new[i])):
            if j<len(child_new[i])-1:
                if child_new[i][j] < variables_range[0]:
                    child_new[i][j] = variables_range[0]
                elif child_new[i][j] > variables_range[1]:
                    child_new[i][j] = variables_range[1]
    return child_new

def mutate(child,variables_range,iter,iter_max):
    u, delta = get_delta(iter,iter_max)
    child_new=[]
    if u < 0.5:
        for i in range(len(child)):
            pos=[]
            for j in range(len(child[i])):
                if j < len(child[i]) - 1:
                    x=child[i][j]+ delta * (child[i][j] - variables_range[0])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    else:
        for i in range(len(child)):
            pos = []
            for j in range(len(child[i])):
                if j<len(child[i])-1:
                    x= child[i][j] +delta * (variables_range[1] - child[i][j])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    for i in range(len(child_new)):
        for j in range(len(child_new[i])):
            if j<len(child_new[i])-1:
                if child_new[i][j] < variables_range[0]:
                    child_new[i][j] = variables_range[0]
                elif child_new[i][j] > variables_range[1]:
                    child_new[i][j] = variables_range[1]
    return child_new

def mutate2(child,variables_range,iter,iter_max):
    u, delta = get_delta2(iter,iter_max)
    child_new=[]
    if u < 0.5:
        for i in range(len(child)):
            pos=[]
            for j in range(len(child[i])):
                if j < len(child[i]) - 1:
                    x=child[i][j]+ delta * (child[i][j] - variables_range[0])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    else:
        for i in range(len(child)):
            pos = []
            for j in range(len(child[i])):
                if j<len(child[i])-1:
                    x= child[i][j] +delta * (variables_range[1] - child[i][j])
                    pos.append(x)
                else:
                    pos.append(100)
            child_new.append(pos)
    for i in range(len(child_new)):
        for j in range(len(child_new[i])):
            if j<len(child_new[i])-1:
                if child_new[i][j] < variables_range[0]:
                    child_new[i][j] = variables_range[0]
                elif child_new[i][j] > variables_range[1]:
                    child_new[i][j] = variables_range[1]
    return child_new


def get_function_values(pop_size,pos_solution,p_solution,users):
    function1_values = []
    function2_values = []
    punish_values = []
    for i in range(pop_size):
        uavs_pos = pos_solution[i]
        c_association = tool.get_c_association(users, uavs_pos, math.ceil(len(users) / len(uavs_pos)))
        p_association = tool.get_p_association(users, uavs_pos, c_association, len(users) * 2 / len(uavs_pos))
        p_rest = tool.get_rest_p(users, uavs_pos, p_association)
        function1_values.append(function_snr(users, pos_solution[i], c_association, p_solution[i]))
        function2_values.append(function_gdop(users, uavs_pos, c_association, p_association, ))
        punish_values.append(punish(p_solution[i], p_rest))
    return function1_values,function2_values,punish_values


def tournament(pos_solution,non_dominated_sorted_solution,rank,distance):
    a1 = random.randint(0, len(pos_solution) - 1)
    b1=a1
    while b1==a1:
        b1 = random.randint(0, len(pos_solution) - 1)
    participants = [pos_solution[a1],pos_solution[b1]]
    best = None
    for participant in participants:
        if best is None or crowding_operator(a1,b1,non_dominated_sorted_solution,rank,distance) == 1:
            best = participant
    return best

def get_f_by_pos(pos,users):
    c_association = tool.get_c_association(users, pos, math.ceil(len(users) / len(pos)))
    p_association = tool.get_p_association(users, pos, c_association, len(users) * 2 / len(pos))
    p = init_p(pos, users)
    function1_value=function_snr(users, pos, c_association, p)
    function2_value=function_gdop(users, pos, c_association, p, )
    return function1_value,function2_value

def crowding_operator(index1,index2,non_dominated_sorted_solution,rank,distance):
    rank1=rank[index1]#index1的等级
    rank2=rank[index2]
    distance1=0
    distance2=0
    for i in range(len(non_dominated_sorted_solution[rank1])):
        if non_dominated_sorted_solution[rank1][i]==index1:
            distance1=distance[rank1][i]
    for i in range(len(non_dominated_sorted_solution[rank2])):
        if non_dominated_sorted_solution[rank2][i]==index2:
            distance2=distance[rank2][i]
    if rank1<rank2 or (rank1==rank2 and distance1>distance2):
        return 1
    else:
        return 0

def drawResultAndIter(resultAndIter):
    f, ax = plt.subplots(1, 1)
    plt.plot(resultAndIter)
    ax.set_xlabel('iteration')
    ax.set_ylabel('fitness')
    plt.show()

if __name__ == '__main__':
    num=[5, 7, 4, 4, 4, 4, 4, 8, 9, 7, 9, 8, 8, 9, 9, 9, 11, 13, 14, 12, 15, 14, 15, 14, 14, 14, 12, 14, 14, 15, 15, 15, 19, 19, 19, 19, 18, 18, 19, 15, 17, 17, 17, 17, 17, 18, 19, 19, 17, 17, 18, 20, 15, 16, 16, 16, 17, 17, 17, 18, 18, 19, 19, 19, 20, 22, 25, 25, 25, 26, 26, 26, 26, 27, 28, 30, 21, 22, 23, 24, 25, 26, 27, 27, 28, 29, 29, 29, 28, 28, 29, 29, 29, 29, 31, 30, 30, 26, 26, 27, 27, 27, 27, 27, 28, 29, 29, 29, 30, 33, 35, 35, 34, 36, 36, 39, 41, 40, 40, 41, 40, 42, 42, 42, 43, 43, 43, 46, 46, 48, 48, 52, 54, 54, 58, 59, 59, 60]
    drawResultAndIter(num)

