代码存档_古地图电子化_地物精度概率矩阵DAG

1_基本要素

共有 n 个节点(地物)，每个节点都含有 n 个参数，第 i 个参数用于表示“本节点(地物)绘制时的相对位置关系参照于第 i 个节点”的概率，此时我们获得了一个 \(n \times n\) 的概率矩阵

2_概率的获取

通过角度精度、距离精度、距离三个参数的加权计算取得，其中角度精度、距离精度与概率呈正相关，距离与概率呈负相关

\[F_d[i,j] = exp(-α \cdot \frac{D_{real}[i,j]}{max(D_{real})}) \]

\(D_{real}[i,j]\)：真实两点间距离

\(max⁡(D_{real})\)：研究区域的最大距离，用于归一化

\(\alpha > 0\)：控制衰减速度

\[S[i,j] = (w_d \cdot E_d[i,j] + w_a \cdot E_a[i,j]) \cdot F_d[i,j] \]

3_DAG生成

通过这种概率矩阵试图生成一个多入边的有向无环图 (OBS、Greedy)

4_做健壮性检验

（方法未定）

点击查看代码

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from typing import List, Tuple, Dict
import random
import matplotlib.colors as mcolors

class EnhancedGraphAnalysisFixed:
    def __init__(self, n_nodes: int = 15, alpha: float = 2.0, w_d: float = 0.3, w_a: float = 0.7, prob_threshold: float = 0.05):
        self.n_nodes = n_nodes
        self.alpha = alpha  # 距离衰减系数
        self.w_d = w_d      # 距离精度权重
        self.w_a = w_a      # 角度精度权重
        self.prob_threshold = prob_threshold  # 概率阈值，控制图的稀疏程度
        self.nodes_data = []
        self.probability_matrix = None
        self.coordinates = []
        self.distance_matrix = None
        
    def generate_data(self) -> List[Dict]:
        """生成模拟数据（重构版）"""
        print("生成模拟数据...")
        
        # 生成节点坐标
        np.random.seed(32)
        self.coordinates = np.random.rand(self.n_nodes, 2) * 100
        
        # 计算距离矩阵
        self.distance_matrix = np.zeros((self.n_nodes, self.n_nodes))
        for i in range(self.n_nodes):
            for j in range(self.n_nodes):
                self.distance_matrix[i, j] = np.linalg.norm(self.coordinates[i] - self.coordinates[j])
        
        self.nodes_data = []
        for i in range(self.n_nodes):
            node_data = {
                'node_id': i,
                'coordinates': self.coordinates[i],
                'angle_precisions': [],
                'distance_precisions': []
            }
            
            # 生成角度精度：0.5-1.0之间的高精度值（第i个为0）
            angle_prec = np.random.uniform(0.5, 1.0, self.n_nodes)
            angle_prec[i] = 0.0
            node_data['angle_precisions'] = angle_prec
            
            # 生成距离精度：0.5-1.0之间的高精度值（第i个为0）
            dist_prec = np.random.uniform(0.5, 1.0, self.n_nodes)
            dist_prec[i] = 0.0
            node_data['distance_precisions'] = dist_prec
            
            self.nodes_data.append(node_data)
            
        print(f"生成了 {self.n_nodes} 个节点")
        return self.nodes_data
    
    def calculate_probabilities(self) -> np.ndarray:
        """计算概率矩阵（基于重构的综合权重公式）"""
        print("计算概率矩阵（重构版）...")
        
        self.probability_matrix = np.zeros((self.n_nodes, self.n_nodes))
        
        # 计算最大距离用于归一化
        max_distance = np.max(self.distance_matrix)
        
        for i in range(self.n_nodes):
            # 计算综合权重S[i,j]
            S = np.zeros(self.n_nodes)
            
            for j in range(self.n_nodes):
                if i == j:
                    S[j] = 0.0
                    continue
                
                # 获取精度数据
                angle_prec = self.nodes_data[i]['angle_precisions'][j]
                dist_prec = self.nodes_data[i]['distance_precisions'][j]
                
                # 距离衰减函数 F_d[i,j] = exp(-α * (D_real[i,j] / max(D_real)))
                distance_factor = np.exp(-self.alpha * (self.distance_matrix[i, j] / max_distance))
                
                # 综合权重 S[i,j] = (w_d * E_d[i,j] + w_a * E_a[i,j]) * F_d[i,j]
                # 使用精度值作为权重（值越小精度越高，所以用1/精度）
                E_d = 1.0 / (dist_prec + 0.01)  # 避免除零
                E_a = 1.0 / (angle_prec + 0.01)
                
                S[j] = (self.w_d * E_d + self.w_a * E_a) * distance_factor
            
            # 归一化概率矩阵 P[i,j] = S[i,j] / Σ_{k ≠ i} S[i,k]
            row_sum = np.sum(S)
            if row_sum > 0:
                self.probability_matrix[i] = S / row_sum
            else:
                self.probability_matrix[i] = np.zeros(self.n_nodes)
                
        return self.probability_matrix
    
    def obs_method_single_root(self) -> List[Tuple[int, int]]:
        """改进的OBS方法：基于概率选择最优父节点，允许多个父节点，确保无环"""
        edges = []
        
        # 步骤1：为每个子节点选择概率最高的前k个父节点
        for child in range(self.n_nodes):
            # 收集所有可能的父节点及其概率
            parent_probs = []
            for parent in range(self.n_nodes):
                if parent != child and self.probability_matrix[parent, child] > self.prob_threshold:
                    parent_probs.append((self.probability_matrix[parent, child], parent, child))
            
            # 按概率排序，选择前k个（k=2或3，避免过多边）
            parent_probs.sort(reverse=True, key=lambda x: x[0])
            max_parents = min(2, len(parent_probs))  # 限制每个节点最多2个父节点
            
            # 步骤2：按概率从高到低添加边，确保不形成环
            for prob, parent, child_node in parent_probs[:max_parents]:
                temp_edges = edges + [(parent, child_node)]
                if not self._has_cycle(temp_edges):
                    edges.append((parent, child_node))
        
        return edges
    
    def greedy_method_with_root(self) -> List[Tuple[int, int]]:
        """改进的贪心方法：允许多个父节点，确保无环，使用概率阈值控制稀疏性"""
        edges = []
        
        prob_edges = []
        for i in range(self.n_nodes):
            for j in range(self.n_nodes):
                if i != j and self.probability_matrix[i, j] > self.prob_threshold:
                    prob_edges.append((self.probability_matrix[i, j], i, j))
        
        prob_edges.sort(reverse=True, key=lambda x: x[0])
        
        # 贪心添加边，允许多个父节点，只检查是否形成环，限制边数以避免过密
        max_edges = min(len(prob_edges), self.n_nodes * 2)  # 限制边数约为节点数的2倍
        for idx, (prob, parent, child) in enumerate(prob_edges[:max_edges]):
            temp_edges = edges + [(parent, child)]
            if not self._has_cycle(temp_edges):
                edges.append((parent, child))
        
        return edges
    
    def _has_cycle(self, edges: List[Tuple[int, int]]) -> bool:
        """检查是否有环"""
        if not edges:
            return False
            
        G = nx.DiGraph()
        G.add_edges_from(edges)
        
        try:
            # 检查是否有环
            nx.find_cycle(G)
            return True
        except nx.NetworkXNoCycle:
            return False
    
    def find_root_nodes(self, edges: List[Tuple[int, int]]) -> List[int]:
        """找出最高级节点（无入度节点）"""
        all_nodes = set(range(self.n_nodes))
        child_nodes = set(child for _, child in edges)
        root_nodes = list(all_nodes - child_nodes)
        return root_nodes
    
    def analyze_graph(self, edges: List[Tuple[int, int]], method_name: str) -> Dict:
        """分析图结构"""
        all_nodes = set(range(self.n_nodes))
        
        # 计算入度和出度
        in_degree = {node: 0 for node in all_nodes}
        out_degree = {node: 0 for node in all_nodes}
        
        for parent, child in edges:
            out_degree[parent] += 1
            in_degree[child] += 1
        
        # 无入度节点（根节点）
        root_nodes = [node for node in all_nodes if in_degree[node] == 0]
        
        # 无出度节点（叶子节点）
        leaf_nodes = [node for node in all_nodes if out_degree[node] == 0]
        
        # 孤立节点
        isolated_nodes = [node for node in all_nodes 
                         if in_degree[node] == 0 and out_degree[node] == 0]
        
        return {
            'root_nodes': root_nodes,
            'leaf_nodes': leaf_nodes,
            'isolated_nodes': isolated_nodes,
            'in_degree': in_degree,
            'out_degree': out_degree,
            'total_edges': len(edges)
        }
    
    def create_colored_edges(self, edges: List[Tuple[int, int]]) -> Dict:
        """为边创建不同颜色"""
        colors = list(mcolors.TABLEAU_COLORS.values()) + list(mcolors.CSS4_COLORS.values())
        edge_colors = {}
        
        for i, (parent, child) in enumerate(edges):
            color = colors[i % len(colors)]
            edge_colors[(parent, child)] = color
            
        return edge_colors
    
    def plot_enhanced_graphs(self, obs_edges: List[Tuple[int, int]], 
                           greedy_edges: List[Tuple[int, int]]):
        """增强版可视化（彩色边 + 最高级节点特殊标注）"""
        fig, axes = plt.subplots(1, 2, figsize=(20, 8))
        
        methods = ['OBS Method (natural DAG)', 'Greedy Method (natural DAG)']
        edges_list = [obs_edges, greedy_edges]
        
        for idx, (method, edges) in enumerate(zip(methods, edges_list)):
            G = nx.DiGraph()
            G.add_nodes_from(range(self.n_nodes))
            G.add_edges_from(edges)
            
            pos = {i: self.coordinates[i] for i in range(self.n_nodes)}
            
            # 创建彩色边
            edge_colors = self.create_colored_edges(edges)
            edge_color_list = [edge_colors[edge] for edge in edges]
            
            # 找出最高级节点（无入度节点）
            root_nodes = self.find_root_nodes(edges)
            regular_nodes = [n for n in range(self.n_nodes) if n not in root_nodes]
            
            # 绘制普通节点
            if regular_nodes:
                nx.draw_networkx_nodes(G, pos, ax=axes[idx], 
                                     nodelist=regular_nodes,
                                     node_color='lightblue', node_size=500)
            
            # 绘制最高级节点（特殊标注：红色大节点）
            if root_nodes:
                nx.draw_networkx_nodes(G, pos, ax=axes[idx], 
                                     nodelist=root_nodes,
                                     node_color='red', node_size=800, 
                                     node_shape='o', linewidths=2,
                                     edgecolors='darkred')
            
            # 绘制边，每条边使用不同颜色
            nx.draw_networkx_edges(G, pos, ax=axes[idx], 
                                 edge_color=edge_color_list, 
                                 arrows=True, arrowsize=20, 
                                 connectionstyle="arc3,rad=0.1")
            
            # 绘制标签
            nx.draw_networkx_labels(G, pos, ax=axes[idx], font_size=10,
                                  font_color='black')
            
            # 添加图例
            legend_elements = [
                plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='red', 
                         markersize=12, label=f'Root Nodes ({len(root_nodes)})'),
                plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='lightblue', 
                         markersize=10, label='Regular Nodes')
            ]
            axes[idx].legend(handles=legend_elements, loc='upper right', 
                           bbox_to_anchor=(1.0, 1.0))
            
            axes[idx].set_title(f'{method}\nEdges: {len(edges)}, Root Nodes: {len(root_nodes)}', 
                              fontsize=14)
            axes[idx].set_xlim(-10, 110)
            axes[idx].set_ylim(-10, 110)
            axes[idx].grid(True, alpha=0.3)
            axes[idx].set_aspect('equal')
        
        plt.tight_layout()
        plt.savefig('enhanced_graphs_with_roots.png', dpi=300, bbox_inches='tight')
        plt.show()
    
    def run_enhanced_analysis(self):
        """运行增强版分析"""
        print("运行增强版图关系分析...")
        
        # 1. 生成数据
        self.generate_data()
        
        # 2. 计算概率
        self.calculate_probabilities()
        
        # 3. 两种方法
        obs_result = self.obs_method_single_root()
        greedy_result = self.greedy_method_with_root()

        print(f"\n增强版分析结果:")
        print(f"OBS: {len(obs_result)} 条边")
        print(f"Greedy: {len(greedy_result)} 条边")
        
        # 4. 分析图结构
        obs_analysis = self.analyze_graph(obs_result, "OBS")
        greedy_analysis = self.analyze_graph(greedy_result, "Greedy")
        
        print(f"\nOBS方法 - 根节点: {obs_analysis['root_nodes']} (共{len(obs_analysis['root_nodes'])}个)")
        print(f"Greedy方法 - 根节点: {greedy_analysis['root_nodes']} (共{len(greedy_analysis['root_nodes'])}个)")
        
        # 5. 可视化
        self.plot_enhanced_graphs(obs_result, greedy_result)
        
        return {
            'obs': obs_result,
            'greedy': greedy_result,
            'obs_analysis': obs_analysis,
            'greedy_analysis': greedy_analysis
        }

if __name__ == "__main__":
    analyzer = EnhancedGraphAnalysisFixed(n_nodes=15)
    results = analyzer.run_enhanced_analysis()
    
    print("\n概率矩阵示例（前5x5）:")
    print(np.round(analyzer.probability_matrix[:5, :5], 3))