import networkx as nx
from collections import defaultdict
import pyodbc
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def compute_horton_strahler(graph):
    # 定义一个字典来存储每个节点的Horton-Strahler顺序
    hs_order = {node: 1 for node in graph.nodes()}  # 初始所有节点的HS顺序为1
    degrees = graph.in_degree()  # 计算入度

    # 从源到出口遍历通过拓扑排序
    nodes = list(nx.topological_sort(graph))  # 获取拓扑排序的节点

    for node in nodes:
        # 从源节点到当前节点的顺序计算
        if degrees[node] == 0:  # 源节点
            continue

        # 获取所有进入当前节点的前驱节点
        incoming = list(graph.predecessors(node))
        orders = [hs_order[pre] for pre in incoming]

        if len(set(orders)) == 1:  # 相同顺序相遇
            hs_order[node] = orders[0] + 1
        else:  # 不同顺序相遇
            hs_order[node] = max(orders)

    return hs_order


def compute_tokunaga_indices(graph, hs_order):
    # 统计N_ij
    n_ij_count = defaultdict(int)

    for node in graph.nodes():
        order = hs_order[node]
        for pre in graph.predecessors(node):
            pre_order = hs_order[pre]
            if pre_order < order:  # 只有i < j的情况
                n_ij_count[(pre_order, order)] += 1

    # 计算Tokunaga指数T_ij
    tokunaga_indices = {}
    for (i, j), count in n_ij_count.items():
        n_j = sum(count for (k, j) in n_ij_count.keys() if j == j)  # N_j
        T_ij = count / n_j if n_j != 0 else 0
        tokunaga_indices[(i, j)] = T_ij

    return tokunaga_indices

def check_self_similarity(tokunaga_indices):
    """检查Tokunaga指数以判断是否为自相似树"""
    ratios = []
    keys = sorted(tokunaga_indices.keys())

    for i in range(1, len(keys)):
        # 计算相邻Tokunaga指数的比率，并加入列表
        prev_key = keys[i - 1]
        current_key = keys[i]
        if tokunaga_indices[prev_key] > 0:
            ratio = tokunaga_indices[current_key] / tokunaga_indices[prev_key]
            ratios.append(ratio)

    # 检查比率是否相等
    return len(set(ratios)) == 1, ratios


def check_tokunaga_tree(tokunaga_indices,threshold=1e-5):
    """检查Tokunaga指数以判断是否为Tokunaga树"""
    # 将Tokunaga指数整理为T_k格式
    tokunaga_k = defaultdict(list)
    for (i, j), T_ij in tokunaga_indices.items():
        k = j - i
        tokunaga_k[k].append(T_ij)

    # 计算每个k对应的平均Tokunaga指数
    for k, values in tokunaga_k.items():
        tokunaga_k[k] = np.mean(values) if values else 0

    k_keys = sorted(tokunaga_k.keys())

    if len(k_keys) < 2:
        return True  # 只有一个k值时，视为稳定

    # 计算相邻T_k的差异
    for i in range(1, len(k_keys)):
        if abs(tokunaga_k[k_keys[i]] - tokunaga_k[k_keys[i - 1]]) > threshold:
            return False  # 超过阈值则不稳定

    return True  # 所有k值均在阈值内稳定

def main(graph):
    # 计算Horton-Strahler顺序
    hs_order = compute_horton_strahler(graph)
    # print("Horton-Strahler Order:", hs_order)

    # 计算Tokunaga指标
    tokunaga_indices = compute_tokunaga_indices(graph, hs_order)
    # print("Tokunaga Indices:", tokunaga_indices)

    # 检查自相似树
    is_self_similar, ratios = check_self_similarity(tokunaga_indices)
    if is_self_similar:
        print("该管网是自相似树。")
    else:
        print("该管网不是自相似树，Tokunaga指数比率:", ratios)

    # 检查Tokunaga树
    is_tokunaga_tree = check_tokunaga_tree(tokunaga_indices)
    if is_tokunaga_tree:
        print("该管网是Tokunaga树")
    else:
        print("该管网不是Tokunaga树。")

#管网数据
.....

G=nx.Graph()
query='SELECT gid,stnod,ednod FROM hzps_lin'
for row in cursor.execute(query):
    G.add_edge(row[1],row[2])
subgraphs = [G.subgraph(c).copy() for c in nx.connected_components(G)]
largest_subgraph = max(subgraphs, key=len)
tree = nx.minimum_spanning_tree(largest_subgraph)
directed_tree = nx.DiGraph()
for u, v in tree.edges():
    directed_tree.add_edge(u, v)

# 关闭数据库连接
conn.close()

# 运行主函数
main(directed_tree)

自相似树是一种特殊的树结构，其Tokunaga指数T_ij仅取决于相应分支顺序的差异，即对于所有i，定义T_k≔T_i(i+k)  (k=1,2,…)。满足特定条件——两个连续分支指数比值为常数T_(k+1)/T_k=c或T_k=ac^(k-1) (a,c>0)的自相似树称为Tokunaga树。请根据这段话，完善判断自相似树和Tokunaga树的代码