import numpy as np
import random
import math
import matplotlib.pyplot as plt
import time

n = 1
N = 200
N_GENERATIONS = 500
cross_rate = 0.5
mutation_rate = 0.1

distance = [
    [0.13, 0.18, 0.13, 0.24, 0.28, 0.28, 0.31, 0.34, 0.54, 0.51, 0.57, 0.54, 0.77, 0.77],
    [0.09, 0.18, 0.13, 0.21, 0.25, 0.25, 0.28, 0.31, 0.51, 0.47, 0.54, 0.51, 0.77, 0.73],
    [0.09, 0.18, 0.13, 0.18, 0.21, 0.21, 0.25, 0.28, 0.47, 0.43, 0.51, 0.47, 0.77, 0.73],
    [0.62, 0.59, 0.60, 0.80, 0.62, 0.19, 0.19, 0.19, 0.44, 0.41, 0.47, 0.41, 0.44, 0.41],
    [0.09, 0.18, 0.13, 0.18, 0.21, 0.21, 0.25, 0.28, 0.47, 0.43, 0.51, 0.47, 0.77, 0.73],
    [0.58, 0.55, 0.58, 0.77, 0.13, 0.15, 0.13, 0.19, 0.41, 0.38, 0.44, 0.38, 0.41, 0.38],
    [0.55, 0.52, 0.55, 0.74, 0.09, 0.12, 0.13, 0.19, 0.38, 0.34, 0.41, 0.34, 0.37, 0.34],
    [0.55, 0.52, 0.55, 0.74, 0.09, 0.12, 0.13, 0.19, 0.38, 0.34, 0.41, 0.34, 0.37, 0.34],
    [0.44, 0.41, 0.41, 0.44, 0.06, 0.06, 0.09, 0.03, 0.19, 0.19, 0.19, 0.22, 0.23, 0.20],
    [0.41, 0.38, 0.38, 0.41, 0.03, 0.03, 0.03, 0.01, 0.13, 0.15, 0.19, 0.19, 0.20, 0.17],
    [0.31, 0.28, 0.28, 0.31, 0.01, 0.01, 0.03, 0.01, 0.09, 0.12, 0.19, 0.15, 0.17, 0.13],
    [0.81, 0.78, 0.81, 1.00, 0.35, 0.38, 0.35, 0.41, 0.55, 0.58, 0.55, 0.15, 0.13, 0.19],
    [0.75, 0.72, 0.75, 0.93, 0.29, 0.32, 0.29, 0.34, 0.49, 0.52, 0.49, 0.09, 0.06, 0.12],
    [0.78, 0.75, 0.78, 0.96, 0.32, 0.35, 0.32, 0.38, 0.52, 0.55, 0.52, 0.12, 0.09, 0.15]
]
# distance = [
#     [0.61,0.66,0.14,0.61,0.87,0.20,0.03,0.83,0.75,0.63],
#     [0.64,0.20,0.76,0.46,0.12,0.59,0.48,0.88,0.63,0.52],
#     [0.05,0.41,0.72,0.65,0.46,0.45,0.20,0.47,0.14,0.98],
#     [0.95,0.30,0.04,0.86,0.93,0.25,0.14,0.54,0.30,0.43],
#     [0.35,0.10,0.02,0.91,0.35,0.05,0.42,0.90,0.33,0.05],
#     [0.11,0.36,0.37,0.58,0.12,0.07,0.46,0.64,0.99,0.92],
#     [0.67,0.46,0.33,0.40,0.59,0.88,0.15,0.12,0.03,0.40],
#     [0.93,0.35,0.38,0.58,0.49,0.20,0.57,0.70,0.81,0.16],
#     [0.27,0.02,0.03,0.92,0.87,0.70,0.54,0.04,0.17,0.38],
#     [0.90,0.28,0.66,0.19,0.13,0.21,0.03,0.97,0.67,0.35]
# ]
# distance = [
#     [0.36, 0.16, 0.53, 0.98, 0.19, 0.22, 0.49, 0.76, 0.32, 0.02, 0.11, 0.38, 0.02, 0.47, 0.37, 0.15, 0.64, 0.17, 0.82,
#      0.14],
#     [0.37, 0.04, 0.80, 0.56, 0.09, 0.88, 0.08, 0.39, 0.19, 0.97, 0.60, 0.52, 0.48, 0.05, 0.79, 0.78, 0.76, 0.49, 0.44,
#      0.50],
#     [0.91, 0.01, 0.51, 0.17, 0.87, 0.78, 0.57, 0.98, 0.03, 0.34, 0.50, 0.93, 0.92, 0.26, 0.41, 0.22, 0.72, 0.01, 0.57,
#      0.88],
#     [0.49, 0.09, 0.46, 0.34, 0.65, 0.53, 0.07, 0.40, 0.91, 0.05, 0.72, 0.97, 0.51, 0.33, 0.20, 0.54, 0.56, 0.07, 0.03,
#      0.07],
#     [0.40, 0.29, 0.88, 0.12, 0.61, 0.99, 0.51, 0.39, 0.77, 0.97, 0.72, 1.00, 0.96, 0.56, 0.07, 0.47, 0.99, 0.40, 0.11,
#      0.18],
#     [0.63, 0.97, 0.94, 0.11, 0.96, 0.80, 0.03, 0.01, 0.55, 0.64, 0.35, 0.69, 0.95, 0.05, 0.86, 0.83, 0.85, 0.82, 0.57,
#      0.05],
#     [0.80, 0.62, 0.10, 0.94, 0.39, 0.63, 0.84, 0.44, 0.52, 0.07, 0.14, 0.03, 0.84, 0.06, 0.63, 0.36, 0.95, 0.67, 0.45,
#      0.79],
#     [0.44, 0.14, 0.04, 0.84, 0.20, 0.51, 0.14, 0.51, 0.42, 0.73, 0.12, 0.50, 0.88, 0.68, 0.53, 0.35, 0.22, 0.02, 0.11,
#      0.28],
#     [0.63, 0.38, 0.45, 0.65, 0.92, 0.19, 0.05, 0.34, 0.25, 0.59, 0.07, 0.43, 0.36, 0.20, 0.58, 0.85, 0.99, 0.53, 0.58,
#      0.24],
#     [0.30, 0.42, 0.86, 0.10, 0.36, 0.03, 0.28, 0.97, 0.04, 0.69, 0.27, 0.72, 0.67, 0.89, 0.95, 0.68, 0.90, 0.91, 0.02,
#      0.53],
#     [0.54, 0.20, 0.68, 0.09, 0.46, 0.14, 0.74, 0.11, 0.17, 0.32, 0.07, 0.41, 0.44, 0.05, 0.91, 0.45, 0.39, 0.59, 0.98,
#      0.54],
#     [0.55, 0.04, 0.26, 0.57, 0.08, 0.75, 0.94, 0.04, 0.97, 0.50, 0.78, 0.23, 0.60, 0.32, 0.31, 0.58, 0.45, 0.18, 0.55,
#      0.59],
#     [0.97, 0.85, 0.29, 0.65, 0.25, 0.69, 0.67, 0.00, 0.31, 0.69, 0.59, 0.54, 0.79, 0.15, 0.28, 0.52, 0.47, 0.47, 0.40,
#      0.92],
#     [0.10, 0.49, 0.01, 0.83, 0.27, 0.29, 0.57, 0.91, 0.83, 0.88, 0.40, 0.52, 0.39, 0.96, 0.95, 0.73, 0.76, 0.81, 0.63,
#      0.78],
#     [0.44, 0.44, 0.77, 0.90, 0.24, 0.85, 0.08, 0.60, 0.09, 0.57, 0.51, 0.73, 0.04, 0.34, 0.04, 0.64, 0.60, 1.00, 0.94,
#      0.62],
#     [0.31, 0.68, 0.48, 0.54, 0.33, 0.73, 0.64, 0.85, 0.57, 0.40, 0.05, 0.52, 0.73, 0.52, 0.50, 0.54, 0.10, 0.58, 0.34,
#      0.50],
#     [0.32, 0.86, 0.28, 0.69, 0.75, 0.50, 0.90, 0.86, 0.67, 0.36, 0.50, 0.83, 0.54, 0.36, 0.87, 0.46, 0.41, 0.10, 0.53,
#      0.16],
#     [0.28, 0.37, 1.00, 0.63, 0.87, 0.46, 0.57, 0.91, 0.18, 0.69, 0.42, 0.45, 0.84, 0.23, 0.68, 0.44, 0.15, 0.93, 0.40,
#      0.94],
#     [0.98, 0.61, 0.75, 0.72, 0.28, 0.83, 0.16, 0.67, 0.32, 0.94, 0.63, 0.03, 0.91, 0.32, 0.73, 0.48, 0.23, 0.72, 0.98,
#      0.56],
#     [0.56, 0.39, 0.95, 0.49, 0.94, 0.80, 0.48, 0.84, 0.46, 0.73, 0.29, 0.84, 0.31, 0.04, 0.30, 0.94, 0.19, 0.95, 0.83,
#      0.81]
# ]


def esc(X, num):
    all_e = []
    all_s = []
    all_c = []
    for _ in range(N):
        x = X[_]
        e2 = c2 = R = 0
        for s in range(num):
            for t in range(num):
                dis = distance[s][t]
                e1 = x[s][t] * dis * 0.28
                e2 += e1
                c1 = x[s][t] * dis * 2.4
                c2 += c1
        l = [(p, q) for p in range(num) for q in range(num) if x[p][q] == 1]
        l2 = []
        for i in range(len(l)):
            l1 = l[i][1]
            l2.append(l1)
        for s1 in range(14):
            for s2 in range(s1 + 1, num - 1):
                for s3 in range(s2 + 1, num - 1):
                    dis1 = distance[s1][l2[s1]]
                    dis2 = distance[s2][l2[s2]]
                    dis3 = distance[s3][l2[s3]]

                    r = 1 - (abs((dis1 - dis2) * (dis1 - dis3) * (dis2 - dis3)) / (max(dis1, dis2, dis3) ** 3))
                    R += r

        e = e2 - (0.6 * R / 364)
        s = 1 - (R / 364)
        c = c2 + (0.3 * R / 364)
        
        all_e.append(e)
        all_s.append(s)
        all_c.append(c)

    return all_e, all_s, all_c


def max_esc(X, num):
    e, s, c = esc(X, num)
    max_e = max(e)
    max_s = max(s)
    max_c = max(c)
    return max_e, max_s, max_c


def getfitness(X):
    pred = []
    all_e, all_s, all_c = esc(X, num)

    max_e = max(all_e)
    max_s = max(all_s)
    max_c = max(all_c)

    for idx in range(N):
        all_e[idx] = all_e[idx] / max_e
        all_s[idx] = all_s[idx] / max_s
        all_c[idx] = all_c[idx] / max_c

        cos = all_e[idx] * (all_e[idx] ** 2) / (math.sqrt(all_e[idx] ** 4 + all_s[idx] ** 4 + all_c[idx] ** 4))
        pred.extend([cos])

    return pred


def gene(num):
    data = [i for i in range(1, num)]
    random.shuffle(data)
    return data


def getPop(num):
    pop = []
    for i in range(N):
        pop.append(gene(num))
    return pop


def get_matrix(pop, num):
    X = []
    for _ in range(len(pop)):
        k = pop[_]
        x = np.zeros((num, num), dtype=int)
        for i in range(num - 1):
            j = k[i]
            x[i][j-1] += 1
        X.append(x)
    X = np.array(X)
    return X


def cross_and_mutation(X, num):
    new_x = []

    for idx, father in enumerate(X):
        if random.random() >= cross_rate:
            idx_point = random.randint(1, num - 2)
            x_one = father[0:idx_point + 1]
            x_two = father[idx_point + 1: num]
            father = np.r_[x_two, x_one]
        if random.random() <= mutation_rate:
            father = mutation(father, num)
        new_x.append(father)
    new_x = np.array(new_x)
    return new_x


def mutation(X, num):
    a1 = random.randint(0, num - 1)
    a2 = random.randint(0, num - 1)
    while True:
        if a1 != a2:
            X[[a1, a2], :] = X[[a2, a1], :]
            break
        else:
            a2 = random.randint(0, num - 1)
    return X


def select_new(old_X, new_X, fitness1, fitness2):
    X = []
    all_fit = np.concatenate((fitness1, fitness2), axis=0)
    fitness_idx = np.argsort(-all_fit)
    fitness_idx = fitness_idx[:N]
    all_pop = np.concatenate((old_X, new_X), axis=0)
    for i in range(200):
        x = all_pop[fitness_idx[i]]
        X.append(x)
    X = np.array(X)

    return X


def print_gene(X, num):
    fitness = getfitness(X)
    max_fitness_index = np.argmax(fitness)
    best_x = X[max_fitness_index]

    best_idx = []
    print('max_fitness:', fitness[max_fitness_index])
    a = [[i, j] for i in range(num) for j in range(num) if best_x[i][j] == 1]
    for index in range(len(a)):
        best = a[index][1] + 1
        best_idx.append(best)
    print("the optimal feasible solution is obtained as：", best_idx)
    return fitness[max_fitness_index]


if __name__ == '__main__':
    time_start = time.time()  # Record start time
    num = len(distance[0])
    max_fitness = []
    pop = getPop(num)
    X = get_matrix(pop, num)
    E = []
    S = []
    s_one = []
    C = []

    for i in range(1, N_GENERATIONS + 1):
        print('generation', i)
        old_X = X
        new_X = cross_and_mutation(old_X, num)
        fitness1 = getfitness(old_X)
        fitness2 = getfitness(new_X)
        X = select_new(old_X, new_X, fitness1, fitness2)

        m_e, m_s, m_c = max_esc(X, num)
        E.append(m_e)
        s_one.append(m_s)
        S.append(m_s)
        C.append(m_c)
        m = print_gene(X, num)
        max_fitness.extend([m])

    score = [max_fitness[0]]
    S_one = [s_one[0]]
    max_E = [E[0]]
    max_S = [S[0]]
    max_C = [C[0]]
    for i in range(1, 500):
        e1 = score[i-1]
        e2 = max_fitness[i]
        if e1 >= e2:
            score.append(e1)
        else:
            score.append(e2)
    print('score:', score[-1])

    for i in range(1, 500):
        e1 = S_one[i-1]
        e2 = s_one[i]
        if e1 >= e2:
            S_one.append(e2)
        else:
            S_one.append(e1)
    print('S:', S_one[-1])

    for i in range(1, 500):
        e1 = max_E[i-1]
        e2 = E[i]
        if e1 >= e2:
            max_E.append(e1)
        else:
            max_E.append(e2)
    print('E:', max_E[-1])

    for i in range(1, 500):
        s1 = max_S[i-1]
        s2 = S[i]
        if s1 >= s2:
            max_S.append(s1)
        else:
            max_S.append(s2)

    for i in range(1, 500):
        c1 = max_C[i-1]
        c2 = C[i]
        if c1 >= c2:
            max_C.append(c1)
        else:
            max_C.append(c2)
    print('C:', max_C[-1])

    x = np.arange(1, 501)
    plt.figure()
    plt.plot(x, score)
    plt.xlabel('times')
    plt.ylabel('score')

    plt.figure()
    plt.plot(x, E)
    plt.xlabel('times')
    plt.ylabel('E')

    plt.figure()
    plt.plot(x, S_one)
    plt.xlabel('times')
    plt.ylabel('1-S')

    plt.figure()
    plt.plot(x, C)
    plt.xlabel('times')
    plt.ylabel('C')

    plt.figure()
    plt.plot(x, max_E)
    plt.xlabel("times")
    plt.ylabel("E")

    plt.figure()
    plt.plot(x, max_S)
    plt.xlabel("times")
    plt.ylabel("1-S")

    plt.figure()
    plt.plot(x, max_C)
    plt.xlabel("times")
    plt.ylabel("C")
    plt.show()
    myset = set(max_E)
    for item in myset:
        print("the %d has found %d" % (item, max_E.count(item)))
    # Record end time
    time_end = time.time()
    # The calculated time difference is the execution time of the program in seconds/s
    time_sum = time_end - time_start
    print(time_sum)