import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 
from sklearn.metrics import silhouette_score, davies_bouldin_score, completeness_score
import warnings

warnings.filterwarnings("ignore")

# 加载数据
data = pd.read_csv(r'C:\Users\29244\Desktop\data.csv') 

# 处理航迹数据
flight_ids = data['flight_id'].unique()
trajectories = {}

for flight_id in flight_ids:
    flight_data = data[data['flight_id'] == flight_id]
    flight_data_sorted = flight_data.sort_values(by='t')
    trajectories[flight_id] = flight_data_sorted

trajectory_features = []

for flight_id, trajectory in trajectories.items():
    trajectory_points = trajectory[['x', 'y', 'z']].values
    trajectory_features.append(trajectory_points) 

trajectory_features = np.array(trajectory_features)

# 计算距离矩阵（欧几里得距离）
N = 196
distances = np.zeros((len(trajectories), len(trajectories)))

for i in range(len(trajectory_features)):
    for j in range(i + 1, len(trajectory_features)):
        trajectory_i = trajectory_features[i]
        trajectory_j = trajectory_features[j]
        
        # 对齐轨迹长度
        if len(trajectory_i) > N:
            trajectory_i = trajectory_i[:N]
        if len(trajectory_j) > N:
            trajectory_j = trajectory_j[:N]
        
        dist = np.sum(np.sqrt(np.sum((trajectory_i - trajectory_j) ** 2, axis=1)))  
        distances[i, j] = dist
        distances[j, i] = dist

# 定义单向距离函数
def unidirectional_distance(trajectory_i, trajectory_j):
    min_len = min(len(trajectory_i), len(trajectory_j))
    # 计算从 trajectory_i 开始到 trajectory_j 结束的距离
    dist = np.sum(np.sqrt(np.sum((trajectory_i[:min_len] - trajectory_j[-min_len:]) ** 2, axis=1)))
    return dist

# 计算单向距离的距离矩阵
distances_unidirectional = np.zeros((len(trajectory_features), len(trajectory_features)))

for i in range(len(trajectory_features)):
    for j in range(i + 1, len(trajectory_features)):
        trajectory_i = trajectory_features[i]
        trajectory_j = trajectory_features[j]
        
        # 对齐轨迹长度
        if len(trajectory_i) > N:
            trajectory_i = trajectory_i[:N]
        if len(trajectory_j) > N:
            trajectory_j = trajectory_j[:N]
        
        # 计算单向距离
        dist = unidirectional_distance(trajectory_i, trajectory_j)
        distances_unidirectional[i, j] = dist
        distances_unidirectional[j, i] = dist

# 手肘法寻找最佳K（欧几里得距离）
inertia = []
K_range = range(2, 11)
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(distances)
    inertia.append(kmeans.inertia_)

# 绘制手肘法图（欧几里得距离）
plt.figure(figsize=(8, 6))
plt.plot(K_range, inertia, marker='o', label='Euclidean Distance')
plt.title('Elbow Method for Optimal K (Euclidean Distance)')
plt.xlabel('Number of Clusters (K)')
plt.ylabel('Inertia')
plt.legend()
plt.show()

# KMeans聚类（欧几里得距离）
best_k = 6  
kmeans = KMeans(n_clusters=best_k, random_state=42)
kmeans_labels = kmeans.fit_predict(distances)

# 计算轮廓系数、DB 指数和聚类精度（欧几里得距离）
sil_score = silhouette_score(distances, kmeans_labels)
db_index = davies_bouldin_score(distances, kmeans_labels)
completeness = completeness_score(kmeans_labels, kmeans.labels_)

print(f"Silhouette Score for K={best_k} (Euclidean): {sil_score:.4f}")
print(f"Davies-Bouldin Index for K={best_k} (Euclidean): {db_index:.4f}")
print(f"Completeness Score for K={best_k} (Euclidean): {completeness:.4f}")

# 手肘法寻找最佳K（单向距离）
inertia_unidirectional = []
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(distances_unidirectional)
    inertia_unidirectional.append(kmeans.inertia_)

# 绘制手肘法图（单向距离）
plt.figure(figsize=(8, 6))
plt.plot(K_range, inertia_unidirectional, marker='o', label='Unidirectional Distance')
plt.plot(K_range, inertia, marker='x', label='Euclidean Distance')  # 比较
plt.title('Elbow Method for Optimal K')
plt.xlabel('Number of Clusters (K)')
plt.ylabel('Inertia')
plt.legend()
plt.show()

# KMeans聚类（单向距离）
best_k_unidirectional = 6  
kmeans_unidirectional = KMeans(n_clusters=best_k_unidirectional, random_state=42)
kmeans_labels_unidirectional = kmeans_unidirectional.fit_predict(distances_unidirectional)

# 计算轮廓系数、DB 指数和聚类精度（单向距离）
sil_score_unidirectional = silhouette_score(distances_unidirectional, kmeans_labels_unidirectional)
db_index_unidirectional = davies_bouldin_score(distances_unidirectional, kmeans_labels_unidirectional)
completeness_unidirectional = completeness_score(kmeans_labels_unidirectional, kmeans_unidirectional.labels_)

print(f"Silhouette Score for K={best_k_unidirectional} (Unidirectional): {sil_score_unidirectional:.4f}")
print(f"Davies-Bouldin Index for K={best_k_unidirectional} (Unidirectional): {db_index_unidirectional:.4f}")
print(f"Completeness Score for K={best_k_unidirectional} (Unidirectional): {completeness_unidirectional:.4f}")

# 定义颜色列表
colors = plt.cm.get_cmap("tab10", best_k)  # 使用tab10色图以获取最多10种颜色

# 绘制三维图（欧几里得距离）
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制所有航迹，使用聚类颜色
for i, flight_id in enumerate(trajectories):
    trajectory = trajectories[flight_id]
    ax.plot(trajectory['x'], trajectory['y'], trajectory['z'], 
           color=colors(kmeans_labels[i]), alpha=0.5, linewidth=0.8)  # 根据聚类选择颜色

ax.set_title('3D KMeans Clustering of Flight Trajectories (Euclidean)')
ax.set_xlabel('X Coordinate')
ax.set_ylabel('Y Coordinate')
ax.set_zlabel('Z Coordinate')
plt.show()

# 绘制整体二维图（欧几里得距离）
plt.figure(figsize=(8, 6))
for i, flight_id in enumerate(trajectories):
    trajectory = trajectories[flight_id]
    plt.plot(trajectory['x'], trajectory['y'], color=colors(kmeans_labels[i]), alpha=0.6, linewidth=0.6)  # 根据聚类选择颜色

plt.title('2D Projection of Flight Trajectories - All Clusters (Euclidean)')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.grid()
plt.show()


# 绘制三维图（单向距离）
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制所有航迹，使用聚类颜色（单向距离）
for i, flight_id in enumerate(trajectories):
    trajectory = trajectories[flight_id]
    ax.plot(trajectory['x'], trajectory['y'], trajectory['z'], 
           color=colors(kmeans_labels_unidirectional[i]), alpha=0.5, linewidth=0.8)  # 根据聚类选择颜色

ax.set_title('3D KMeans Clustering of Flight Trajectories (Unidirectional)')
ax.set_xlabel('X Coordinate')
ax.set_ylabel('Y Coordinate')
ax.set_zlabel('Z Coordinate')
plt.show()

# 绘制整体二维图（单向距离）
plt.figure(figsize=(8, 6))
for i, flight_id in enumerate(trajectories):
    trajectory = trajectories[flight_id]
    plt.plot(trajectory['x'], trajectory['y'], color=colors(kmeans_labels_unidirectional[i]), alpha=0.6, linewidth=0.6)  # 根据聚类选择颜色

plt.title('2D Projection of Flight Trajectories - All Clusters (Unidirectional)')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.grid()
plt.show()



# 输出每个聚类包含的航迹代号和数量（欧几里得距离）
cluster_dict = {label: [] for label in range(best_k)}

for flight_id, label in zip(trajectories.keys(), kmeans_labels):
    cluster_dict[label].append(flight_id)

for cluster_num, flight_ids in cluster_dict.items():
    print(f"Cluster {cluster_num + 1} (Euclidean):")
    print(f"  Number of trajectories: {len(flight_ids)}")
    print(f"  Trajectories: {', '.join(map(str, flight_ids))}\n")
# 输出每个聚类包含的航迹代号和数量（单向距离）
cluster_dict_unidirectional = {label: [] for label in range(best_k_unidirectional)}

for flight_id, label in zip(trajectories.keys(), kmeans_labels_unidirectional):
    cluster_dict_unidirectional[label].append(flight_id)

for cluster_num, flight_ids in cluster_dict_unidirectional.items():
    print(f"Cluster {cluster_num + 1} (Unidirectional):")
    print(f"  Number of trajectories: {len(flight_ids)}")
    print(f"  Trajectories: {', '.join(map(str, flight_ids))}\n")