BHI_assembly_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import os\n",
    "\n",
    "# Load the .mat files\n",
    "data1 = sio.loadmat('BHI_assems.mat')\n",
    "data2 = sio.loadmat('BHI_inferred_params.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "relabus = data1['relabus'].T\n",
    "relabus = relabus[:,105:]\n",
    "relabus[np.isnan(relabus)] = 0\n",
    "\n",
    "R_inf = data2['R_inf']\n",
    "G = (R_inf > 0).astype(float)\n",
    "G = G[:, 13:]\n",
    "\n",
    "# g2f = np.array([1, 2, 3, 4, 3, 3, 3, 1, 3, 1, 3, 4, 3, 4, 5])\n",
    "\n",
    "# G= np.zeros((15, 5))\n",
    "# for i in range(15):\n",
    "#     G[i, g2f[i] - 1] = 1  # 调整索引，MATLAB 索引从 1 开始，Python 索引从 0 开始\n",
    "\n",
    "inds2del = (relabus[1, :] == 0)\n",
    "relabus = relabus[:, ~inds2del]\n",
    "\n",
    "\n",
    "relabus = relabus / np.sum(relabus, axis=0, keepdims=True)\n",
    "N, mm = relabus.shape\n",
    "N, M = G.shape\n",
    "\n",
    "mx = np.mean(relabus, axis=1)\n",
    "n1_orig = np.sum(G, axis=0)\n",
    "dm_origs = np.zeros(M)\n",
    "\n",
    "for j in range(M):\n",
    "    dm_origs[j] = np.log10(np.mean(mx[G[:, j] == 1])) - np.log10(np.mean(mx[G[:, j] == 0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 15)"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 设置要绘制的柱状图数目\n",
    "mm2plot = 30\n",
    "\n",
    "# 对 relabus 的行求和并按降序排序，获取排序后的索引\n",
    "sorted_inds_x = np.argsort(np.sum(relabus, axis=1))[::-1]\n",
    "\n",
    "# 随机选择 mm2plot 个索引进行绘制\n",
    "inds2plot = np.random.permutation(mm)[:mm2plot]\n",
    "\n",
    "\n",
    "# 转置 relabus 并选择要绘制的部分\n",
    "xthis = relabus.T\n",
    "xthis = xthis[inds2plot, :][:, sorted_inds_x]\n",
    "\n",
    "\n",
    "xthis.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbn0lEQVR4nO3dbYyV5ZkH8GsAlZcykYiroiCuH1ZwSWMXIi+2BQoRB2uddDVaQ6VaJQY3imhKU4sIKrUY2q1CujVsLKKt3aZYgyJrq6CUviyWadS6UA1YWEJTa1kpSl+G2Q8WVurMyMyZ6zzn5ff7BOfM3Nd9zpy3/7nu534a2tra2gIAAADocb2KngAAAADUKqEbAAAAkgjdAAAAkEToBgAAgCRCNwAAACQRugEAACCJ0A0AAABJhG4AAABI0qeowgcPHozdu3fHwIEDo6GhoahpAAAAQJe1tbXFvn37YsiQIdGrV8f97MJC9+7du2Po0KFFlQcAAICS7dy5M0477bQOry8sdA8cODAi3plgY2NjSWMdeOqlnphSh/pOPbviav/PbZtS60ZEnHr7+HYv/8aNG1LrXvvVj7Z7eZG3ees/jU6v/Q/Pb2738qIeY9l1O6ud/bfu6O8cEXHvzEtSa//LA//R7uX//tq/ptaNiLjq9Bvav+KZBbmFJ7U/fpGPsSLv73p83yjqeRVRn6+hP338v1Prnjv9rHYvz/58ENHxZ4Tsx1hEJ4+zgl5Di7zNRX72rsfPoUV9Lsp+LYno+PWkyM/eXfHmm2/G0KFDD2fbjhQWug8tKW9sbCw5dL/5X3/piSl16O8+2fH8iqr95nEDUutGRId/l37H5tbuqG6Rt3nMFb9Jrx0d1D62/wdSy/bt6P5OfmxHdPz4vn/3fal15zZO6/C6vscck1q7w+fVwL6pdTurHQOOyy7c7sVFPsaKvL+Lek5nP68iOn5uFfW8iijufTr77xzR8d96QHLtoj4fdFY7+zHWWe2iXkOLvM1FvY5FFPc59IV/fDW1bkTEiA5eQ4v6XJT9WhLR8f39gd69C6vdHe93uHRhoZvSPLL97vQac+PD6TW6oh5vc0TE6xv2po5/2gXtX16v9/e4T3256ClQ44p6ThepHp9Xyx74fHqNuResSa8B7anH17F6NH7LuPwin/jfdi8ecd/F+bXLSOgGgDLK/kKrEr/MKlI93t+DPvfJ3AKfeDl3fCpePT6voBQ1EbqLfOJ70SGbxxi1ql5XUwBAV9TjaqFaUxOhGwDonA9tAFAMoRuoOMIBAAC1QugG2iX4Aj3FMcbl9eyTV6WOPyJ1dIDaI3RXqXoMRPV4myPq93YXRTggm+d0ebm/gWq3Ze2e1PHHf2Jk6vgI3QBQF7I/tEX44PZuwj7Z/n7OisJqF/n4nn1yc3KF9k9hBaWoidDtjQ0AgEJsXJ47/pTFueMD6WoidANQfXxhCgDUA6G7RH0H3VT0FKDmOHYJaouNvQCoZ72KngAAAADUKp1ugHfRkQMAoCcJ3QAAlGzIzLOLngJARRK6qRpOd1Ne7m/IYc8CAKgvQjdQcWafszW5wuTk8QEA4B1CNwDUgcnrZ5ehystlqAFQXzQjqp/QXaUs/S2v1dO2pddo7uByS1EBqo/3aQAOEbpL5JsnAAAAOiJ002Wz//mp5Artf9GQ/wVHx7Upr5e/8rXU8UfMui51fI6OTiBAdbIKD7pG6Aba5UsOyJF/bHX7x1Vnn4M+ojLPQz9+z4bkCsIBkEszovr1KnoCAAAAUKtqotNdj0tcdCGBnrL6Azemjt/RJoFQy7xPA3CITjcAAAAkqYlON1Bbso89rcTjTqGW3bdhZer4s4bfkjo+QJF8Lqp+QjdUOKelAwAqic8m0DU1Ebo98QGoFjoWAFBfaiJ0Az0v+/QUEU5RUe9sNAUA1AOhGwAAgMqxcXl+jSmL82v8ld3LAQAAIIlON1XDcmcAAKDaCN0A7+LURlBbsr+w9WUtAO9H6AYAoKrNvWZU0VMA6JDQXaUstQaqndcxAHh/VuFVP6G7RJatUauyzyUc4XzCAADUPqEbAKCGzHrljdwCw3OHB6g1ThkGAAAASYRuAAAASGJ5OQDUgeyNeCJsxvNuNgoE4BChm6phYy8AgOLZSBi6RuiGCueNDQAAqldNhO5/O/8PqePPSh0dAMiSvUrKCikor9XTtqWO35w6OvWqJkJ3PbLUGgAAoPIJ3SXyDToAAAAdEboBKIQVOwBAPRC6gYqTfWojpzUCAKBcehU9AQAAAKhVQjcAAAAksbwcoALMeuWN/CLD80sAAHAkoRsAoIfZKBCAQywvBwAAgCQ63VDhnAu+Tmxcnl9jyuL8GgDUvEs/nxshXkgdHcpPpxsAAACS6HSXyPmEyyf7vo5wfwMAAD1LpxsAAACS6HQDADXLijQAiiZ0VylLrQEAACqf5eUAAACQRKcbAMrIcmcAqC863QAAAJCkJjrdugYAAABUopoI3QAA1K/hPzgvvcaOKeklgBpleTkAAAAk0ekGoBBOfQgA1AOhG4C6M+uVN/KLDM8vAQBUPsvLAQAAIInQDQAAAEksLwfaZfktAACUTqcbAAAAkgjdAAAAkMTycgAg1Qvbf130FACgMDrdAAAAkESnGwCglmxcnjv+lMW54wPUGJ1uAAAASCJ0AwAAQBKhGwAAAJII3QAAAJDERmpQ4e7bsDJ1/FnDb0kdv9o4tVH5vHDqhUVPAQAgnU43AAAAJNHpBiqODigAALVC6AYA6GHZhwZFODwIyukn/7k2vcaEq2en16AYlpcDAABAEp1uAAAAjjD8wMPpNXakV6gMOt0AAACQRKcbAACgE+M+9eWip0AV0+kGAACAJEI3AAAAJLG8nC4b/oPzUsffMSV1eAAAgLLR6QYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKnDAOg7mSf+jDC6Q8BgHcI3QAAAJ3YsnZPeo3xnxiZXoNiCN0AADVk9bRtqeM3p44OUHsc0w0AAABJdLoBgFQ6rwDUM6EbAACgE7PP2VqGKpPLUIMiCN0A76IjBwBAT3JMNwAAACQRugEAACCJ5eUAANBNww88nDr+jtTRgXIQuqHCvXDqhUVPAQAA6CahG2jfxuX5NaYszq8BAAAFEroBAACoGNmHbUSU99ANoRsA6oBDVQC67+WvfC29xohZ16XXoBh2LwcAAIAkOt1UDV0aAACg2uh0AwAAQBKd7iql6wsAAFD5dLoBAAAgiU43AABAhbLCtfoJ3QAAAJ149smr0muMSK9AUYRuAKBm6RABUDTHdAMAAEASnW4AAOCoWUFSZhuX544/ZXHu+Oh0AwAAQBahGwAAAJII3QAAAJBE6AYAAIAkNlIDqADDDzycXmNHegXgEBtNAXCITjcAAAAkEboBAAAgidANAAAASRzTDQAA0In7NqxMrzFr+C3pNSiGTjcAAAAk0ekukd1JAQAA6IhONwAAACQRugEAACCJ0A0AAABJhG4AAABIInQDAABAEqEbAAAAkjhlGACUkVNNAkB9EboBKMTqD9yYXqM5vQIAQOcsLwcAAIAkQjcAAAAkEboBAAAgidANAAAASWykBrRr9bRt6TVscgUAQK3T6QYAAIAkNdHpds5TAAAAKlFNhG4AAIo159HN6TWax05PrwHQ0ywvBwAAgCQ63QAAABUqe3NbG9vm0+kGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSOE83AABw9DYuzx1/yuLc8aHMhG6odN7YAACgalleDgAAAEmEbgAAAEgidAMAAEASoRsAAACS2EitWmVvrhVhgy0AAIASCd0AAEBVmPPo5tTxm8dOTx2f+mR5OQAAACTR6QYAgCqT3fGN0PWFniJ0Uz0cxw4AAFQZy8sBAAAgidANAAAASYRuAAAASCJ0AwAAQBIbqZUqe3MvG3sBAABULaEbAMrJl7UAUFeEbgCgdvmSA4CCOaYbAAAAkgjdAAAAkMTycgCAnpa9rD3C0naAKqHTDQAAAEmEbgAAAEhieTlAnZvz6ObU8ZvHTk8dHwCgkul0AwAAQBKhGwAAAJII3QAAAJCkNo7pzj4th1NyAAAA0A21EboBgIplsz4A6pnl5QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgiY3UAN7Fhk8AAPQkoRsAADhqq6dtSx2/OXV0KD/LywEAACCJTneVyv6GMcK3jJXCt8kAAFC9dLoBAAAgiU43AABAJ17Y/uuip0AV0+kGAACAJEI3AAAAJLG8HADqwcbl+TWmLM6vAQBVRugukZ2ly8eO7QAAQLWxvBwAAACSCN0AAACQROgGAACAJEI3AAAAJLGRGgCUkQ04AaC+6HQDAABAEqEbAAAAklheDgDULMv5ASia0A1HYc6jm9NrNI+dnl4DAAAoL8vLAQAAIIlONwBAD8te1h5haTtAtdDpBgAAgCRCNwAAACQRugEAACCJ0A0AAABJhG4AAABIYvdyAACACjXn0c2p4zePnZ46PjrdAAAAkEboBgAAgCRCNwAAACRxTDcAAEAnVk/bll6jOb0CRenxTvdf/vKXnh4SAAAAqlKXQ/c111wTb731VrvXtbS0xOjRo0ueFAAAANSCLi8vf+ihh+LZZ5+NVatWxZgxYyIiorW1Ne644464884749RTT+3xSQIApbE0EgCK0eVO989//vNobGyMCRMmxKJFi2LLli0xZsyYuOOOO+L666+Pl156KWOeAAAAUHW6HLrPOuus+MlPfhILFiyIhQsXxujRo+OPf/xj/OxnP4ulS5dG//79M+YJAAAAVadbG6m9+uqr8fTTT0dra2uccsopsX379njyySejtbW1p+cHAAAAVavLx3QvXLgwFi9eHCeccEKsWbMmpk6dGrfddlvMnz8/Hnnkkbj//vsPH+tdLtnHqTlGDQAAgO7ocqd7wYIFcckll8SLL74YTU1Nccwxx8Rdd90VmzZtij/96U8xfvz4jHkCAABA1elyp/v73/9+fPzjH3/P5WPGjImWlpb44he/2CMTAwAAgGrX5U53e4H7kGOPPTbuvvvukiYEAAAAtaJbG6kBAAAA70/oBgAAgCRCNwAAACTp8kZqUJQ5j25Or9E8dnp6jWpR6P29cXlu4SmLc8cHAIC/0ukGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASpwyrUk6fBQAAUPl0ugEAACCJ0A0AAABJhG4AAABI4phuAADgqGXvLWRfIWqN0A0VzhsbAABUL8vLAQAAIInQDQAAAEmEbgAAAEgidAMAAEASG6kBAAB0Intj2wib29YynW4AAABIInQDAABAEqEbAAAAkgjdAAAAkMRGagAANSR7wyebPQF0TU2Ebm8uAAAAVKKaCN1AbVk9bVvq+M2powMAwP8TugGgjKzOKi/3NwBFs5EaAAAAJNHpLpFv0AEAAOiITjcAAAAk0ekGAOhh2SvhIqyGA6gWOt0AAACQRKcbgELoBAIA9UCnGwAAAJII3QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKhGwAAAJII3QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKhGwAAAJII3QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKhGwAAAJII3QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKhGwAAAJII3QAAAJCkT9ETAADyzXl0c3qN5rHT02sAQLXR6QYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKhGwAAAJII3QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgidANAAAASYRuAAAASCJ0AwAAQBKhGwAAAJII3QAAAJBE6AYAAIAkQjcAAAAkEboBAAAgidANAAAASfoUPQGAvzXn0c2p4zePnZ46PgAAHKLTDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJH2KKtzW1hYREW+++WbJYx3841slj9GZzuZYVO3sukXWrsfbXGRtt7kyatfj/V3Lt7nI2m5zZdSux/u7lm9zkbXd5sqoXY/3dy3f5s5qd2eMQ9m2Iw1t7/cTSXbt2hVDhw4tojQAAAD0iJ07d8Zpp53W4fWFhe6DBw/G7t27Y+DAgdHQ0FDEFAAAAKBb2traYt++fTFkyJDo1avjI7cLC90AAABQ62ykBgAAAEmEbgAAAEgidAMAKRzBBgBCNwCUbPPmzTFjxowYNmxY9OvXL84888y49tprY/v27UVPrSx27NgRDQ0N8cADD0RExN69e+PTn/50PPfcc8VODAAqgNANACVYtmxZjBs3Ln7zm9/El770pVi7dm3Mmzcv1q9fH6NHj45f/OIXRU8x3SmnnBI//vGPY/r06RER0dLSEg8++GAcPHiw4JkBQPHsXg4A3fSjH/0oPvrRj8b1118fX/3qV4+47re//W2cc845cdJJJ8Xzzz9fzAQLsn79+pg0aVI888wzMXHixKKnAwCF0ukGgG5asmRJHH/88XHXXXe957oTTzwxli5dGhdffHHs378/WltbY/ny5TFq1Kjo169fDBs2LObNmxcHDhw4/DszZ86MadOmxTe+8Y0488wzo1+/fjFhwoTYtm1brFmzJkaNGhX9+/ePc889N1paWkr+vYkTJ74nFK9fvz4aGhpi/fr1ERHxwAMPRJ8+feKnP/1pjBs3Lvr27Runn3563HPPPYd/593Lyw8F7oiISZMmxcSJE2PZsmXR0NAQ27ZtO6LWqlWronfv3rFz585u/gUAoPL1KXoCAFCN2traYt26dXHRRRdF//792/2ZSy+99PC/P/vZz8bKlStj3rx58eEPfzi2bNkSt99+e2zZsiWefPLJaGhoiIiITZs2xe7du2Pp0qXx9ttvx3XXXRdNTU3R0NAQCxcujAEDBsSsWbPiiiuuiJdeeunw+N39vaNx8ODBuPTSS+Omm26KO++8M1asWBG33HJLjBo1Ks4///wjfvZDH/pQLFu2LGbPnh3Lli2LiRMnxpAhQ+Lmm2+OBx98MBYtWnT4Z7/5zW/Gxz72sRg6dGiX5gMA1UToBoBueP311+PAgQNxxhlnvO/P/vKXv4wVK1bE4sWLY968eRERMXXq1BgyZEjMmDEj1q5dG01NTRERsW/fvvjOd74TZ511VkREbNiwIb7+9a/HD3/4w5g8eXJERPzqV7+Km2++Ofbu3RvHH398Sb93NNra2mL+/Plx9dVXR0TEhAkT4nvf+16sWbPmPaG7sbExRo4cGRERI0eOPPzv5ubmWLVqVSxcuDAaGhpi165d8fTTT8eqVauOeh4AUI0sLweAbujT553vrVtbW9/3Zzds2BAREZdffvkRl1922WXRu3fvw0u5IyIGDRp0ODhHRJx00kkREXHuuecevuyEE06IiHd2CS/1947WuHHjDv/7uOOOixNPPDH2799/1L9/9dVXx44dOw7vaL5y5coYOHBgNDc3d3kuAFBNhG4A6IZBgwbFwIED47XXXuvwZ/bv3x+///3v44033oiIiJNPPvmI6/v06RODBw8+IgQ3Nja2O9aAAQM6nU93f+9o/e0S+l69enVpd/LJkyfHGWecEStXroyId0L3ZZddFn379u2R+QFApRK6AaCbzj///HjmmWeO2Azt3e6///4YPHjw4f/v2bPniOv//Oc/x+uvv37Ez5RTQ0PDezr1f/jDH9JqzZw5M1avXh3PP/98bN26NT7zmc+k1AKASiJ0A0A3zZ07N373u9/Frbfe+p7r9uzZE/fcc0+MHDny8BLqb33rW0f8zLe//e1obW2N8847ryzz/VuNjY2xa9euIy7buHFjyeP27t273ctnzpwZe/fujblz58aIESOOWPoOALXKRmoA0E1jx46NRYsWxa233hovv/xyXHnllTF48OB48cUXY8mSJfH222/HU089FSNGjIgrr7wy5s+fH2+99VZ85CMfiZaWlliwYEFMmjQppk2bVsj8L7zwwnjsscfipptuiosuuiiee+65w8u/S3Fok7bHH388Bg0aFB/84AcjImLYsGExderUWLduXdx9990l1wGAaqDTDQAl+MIXvhBPPPFERETceOON0dTUFPfee29ceOGF0dLSEiNGjIiIiBUrVsRtt90WDz30UDQ1NcWyZcvihhtuiCeeeCJ69Srm7fiqq66Kz33uc/Hwww9HU1NTbNq0Kb773e+WPO7ZZ58dl19+edx3331xxRVXHHHd9OnTo3fv3jFjxoyS6wBANWhoa2trK3oSAEB9uOCCC6Jv376xevXqoqcCAGVheTkAkG7RokWxdevWWLduXY8cNw4A1ULoBgDSPfbYY/HKK6/EkiVLYvz48UVPBwDKxvJyAAAASGIjNQAAAEgidAMAAEASoRsAAACSCN0AAACQROgGAACAJEI3AAAAJBG6AQAAIInQDQAAAEmEbgAAAEjyf37avhEC25udAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 假设 data 是一个形状为 (30, 101) 的二维数组\n",
    "data = xthis\n",
    "\n",
    "# 归一化数据，使得每列的和为1\n",
    "\n",
    "# 创建颜色映射对象，包含30种颜色\n",
    "colors = []\n",
    "colors.extend(plt.cm.tab20(np.linspace(0, 1, 20)))\n",
    "colors.extend(plt.cm.tab20b(np.linspace(0, 1, 20)))\n",
    "colors.extend(plt.cm.tab20c(np.linspace(0, 1, 20)))\n",
    "colors.extend(plt.cm.Set1(np.linspace(0, 1, 9)))\n",
    "colors.extend(plt.cm.Paired(np.linspace(0, 1, 12)))\n",
    "colors.extend(plt.cm.Set2(np.linspace(0, 1, 8)))\n",
    "colors.extend(plt.cm.Accent(np.linspace(0, 1, 8)))\n",
    "colors.extend(plt.cm.Set3(np.linspace(0, 1, 12)))\n",
    "colors = colors[:100]\n",
    "colors.append(([0,0,0,1]))  # 确保总共101种颜色\n",
    "\n",
    "# 创建堆叠柱状图\n",
    "fig, ax = plt.subplots()\n",
    "x = np.arange(30)  # x轴上的点\n",
    "bottom = np.zeros(30)  # 底部值的初始化\n",
    "\n",
    "# 绘制堆叠柱状图\n",
    "for i in range(N):\n",
    "    ax.bar(x, data[:, i], bottom=bottom, color=colors[i], edgecolor='none')\n",
    "    bottom += data[:, i]  # 更新底部值，使得下一个柱子在上一个柱子的顶部\n",
    "\n",
    "# 设置坐标轴和图形属性\n",
    "ax.set_xlim([-1, 30])\n",
    "ax.set_ylim([0, 1])\n",
    "ax.set_yticks([])\n",
    "ax.set_xticks([])\n",
    "ax.set_ylabel('x', fontsize=12, fontname='Arial')\n",
    "ax.set_xlabel('Community', fontsize=12, fontname='Arial')\n",
    "ax.tick_params(width=0.5, length=0.0)\n",
    "ax.minorticks_off()\n",
    "ax.tick_params(direction='out')\n",
    "fig.set_size_inches(10, 5) \n",
    "fig.tight_layout() \n",
    "\n",
    "# 显示图形\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8QAAAHbCAYAAAANs1OQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAY/ElEQVR4nO3dbYxV9Z3A8d8dxgjiTECwtBRQ9MU6NKSx66489GGga9WBuKXpGpCgFF0p0Y1d1JSklFqp2i7GuurYrl2sIq1N15YNwQdilKFauyayQPrgYtuohRJ8aCUC6qaB2RfWWafOKDNzz5y59/f5vIK5d87/f+6de+/53nPuuZXOzs7OAAAAgGQayp4AAAAAlEEQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKTUWNbAR44cib1790ZTU1NUKpWypgEAAECd6ezsjAMHDsT48eOjoaH3/cClBfHevXtj4sSJZQ0PAABAndu9e3dMmDCh18tLC+KmpqaIeHOCzc3N/V7Ok/f/T7Wm1KMz55w2pMb9r3V3FDrutAsv7fHn2W7nIscdauu666/PKHTcv9r2lHFLGres54u4898LHTeWXNLzsM//a7HDnnTFkBq3rL/lItd3qK1rPb4Gvdu4mZ4fix53qP1NlXUbZ9tuvXXxPxQ25j/d9R+DPuZQHLcvXn311Zg4cWJXd/amtCB+6zDp5ubmAQXxWa9uq9aUetb8tz3+eORxxxc7bC+3ycgRI8oZt6z1rcNxh9q6Hj9smHHrdNxn53YWOuanenvuLvh5KnoZd0TT8IKHHVrj/vfD/1jouH9TwvoOpcdPRHnPy2Vt22R6fix63KH2N1XWbfzzxzoKHfdTl1/V48/Lup2HH3NMXY05FMftj/f6eG5pQQxAfXr6m7cUuvyWpcsKXT4AkIcg5qjM2Le14BGmFLx8AKgN3lQCGDw1H8ReNACAIngzGKD+1XwQl6WsF8nK2FkFj0tZbHgBAMDgEsQAVNVPHlpS6PJbCl06MNg8ZwBlEsQAAKRTZIiLcKgdDWVPAAAAAMpgDzEk51A1AACyEsQ1JtuJl8pa32LH7XnMfzv7YIFjRiwtdOkA9cc3WVBt2bbjyuJ2pi8EMUfFRgEAZXEkC1ALho9enmLMeiOI+8nePKhNNqwBAHiLIAYAgDrmzWDonSCuMWXtmfZESrVl+5vKtr4AALVAEAMA4I07ICXfQwwAAEBKghgAAICUHDINAADAUZt+wb+UPYWqsYcYAACAlGp+D7ETQAAAAG95+pu3FLr8lqXLCl0+g6vmgxiKUOTXW/X21VYAEOHNfoDBJIgZ0sr63mUAAOiLy07fVeDSZ5cwZu/j1hNBDAAAg8ChvPWtyPu3t/vW39TAOakWAAAAKdlDDADQA5/lBah/ghgAAKi64aOXlz0FeE8OmQYAACAlQQwAAEBKghgAAICUfIYYAACoG5lOiJdpXYsiiAEAYBCIl/pW5P3rvi2OQ6YBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASKmx7AkAAABQO7Y/uK/Q5c/4+ymFLv/t7CEGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAIKXGsicAAADUn8tO31XwCLMLXj4Z2EMMAABASvYQAwAAVff0N28pdPktS5cVunxysIcYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASMlZpgEAADhq9fQd0/YQAwAAkJI9xAAAQNX95KElhS6/pdClk4U9xAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAAp+R5iAOrCbVvXFbr8pSdfXejyAYDBJ4gBoAZ5AwAABk4QA8AAZAvTItd3qK0rAPVPEAPUqWyhBgDQV06qBQAAQEqCGAAAgJQcMg1AVTlUG6gFPg8PRAhiAABK5E00oEwOmQYAACAle4hhiCjrHXLvzANAffNaD70TxP0kXgAAAGqbIOaoCHHqhb9lqD0etwAURRADpci2gZttfQEAaoEgZkgTEQAwOLzmAhkJYgAAAI7a09+8pdDltyxdVujy304QQw+KfJfcO+QAADA01HwQO7wHAKgntm0ABk9D2RMAAACAMghiAAAAUqr5Q6YBAIrg0GWA+mcPMQAAACnZQwwAANQNR3fQF4IYAAAGgVCjXvzkoSWFLr+l0KV355BpAAAAUrKHGAAAoAY56mDgBDEAAMAAFRmnGcK0LA6ZBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUmosewIAAED9uW3rukKXv/TkqwtdPjnYQwwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUmosewIAAADUjtu2rit0+UtPvrrQ5b+dPcQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUmosewIAQN/9/INzy54CANQ8QQxAXRCIAEBfOWQaAACAlAQxAAAAKQliAAAAUvIZ4hqT7TNy2dYXqsnjB6B3mZ4jM60r9JUgBoAByLahmW19AahvghiAqhJMAECtEMT9lG2DL9v6lsFtDEBGXv+AMgliSM6GCAAAWQli6IFIrF/uW6g9HrcAFEUQA6XItoGbbX2B2uN5Csio5oPYkzcAUATbGAD1r+aDGACgngjx+uW+HRxuZ/qioewJAAAAQBkEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSr10CAAAYoDK+7slXTA2cPcQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASKmx7AkAAABQO37+wbllT6Fq7CEGAAAgJXuIAQCAqqunvYjUL3uIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACClxrInAO/q8duLXf7f3VDs8gEAgCHLHmIAAABSEsQAAACkJIgBAABIyWeIoSdFfnbZ55aBanCOBag9HreDw+1cvDq6jQUxR6eO/ugBAFKxHQe96lcQ/+53v+v1soaGhjj++ONj1KhR/Z1TbfDEQrX5m4KB8Rii2vxNAdS9fgXxySefHJVK5V2vc8IJJ8QVV1wRK1eu7NfE6IUXZ2Co8zwF1AIfjwKin0F89913x6WXXhqtra2xYMGCGDduXLz44ovxox/9KDZt2hRf/vKX4+DBg3HdddfFmDFjYtmyZdWe9/+z4QUA1BPbNgCDpl9BfO+998b8+fPju9/9brefL1q0KJYtWxbbtm2LjRs3xgknnBDf+ta3ig1iYGDK2vAybvHj2qgeHNlu5zL2qmW7jcuS6fkR4M/69bVLHR0dccEFF/R42Wc+85l45JFHIiJixowZ8dvf/rb/swMAAICC9CuIx4wZEzt37uzxsp07d0Zzc3NERBw8eDBGjhzZ/9kBAABAQfp1yPTChQtj1apVccwxx8RnP/vZeN/73hcvvPBCbNiwIa655pr4/Oc/H6+88krcfPPNMW3atGrPGQAAAAasX0H8ta99LV588cVYvnx5LF++vOvnDQ0NsWTJkrj++uvjvvvui+3bt8ejjz5atckCAABAtfQriBsbG+POO++MlStXxpYtW+Kll16KCRMmxMyZM2Py5MkREXHuuefG73//+zj22GOrOmEAAACohn4F8VtOOeWUOOWUU3q8bPTo0QNZNAAAQO1wFv6a1K+TagEAAECtE8QAAACkJIgBAABISRADAACQkiAGAAAgpQGdZRoAAKBHzoBMDbCHGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqNZU8AAACA2rHhnGcKXf68QpfenT3EAAAApCSIAQAASEkQAwAAkJLPEAMAQB2rp897QrXVfBB7gNc39y8AZSnrNchrX/1y3w4OtzN9UfNBnI0HOABQBNsYQEaCGChFtj0vNjQBhpYin5c9J0PtEMTQg0wvkkINalOmx26mdYV64rFLLRDE/ZTtAZ5tfYHa43kKgDKVsUPFa9/ACWIYIjyhUW0ODwdqgecMoEy+hxgAAICU7CEGgAHItncr0zkWAKh/9hADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKTWWPQEAAKD+/PN/PlXo8udNm1Po8uldPd239hADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKTnLNADAEFJPZ28FGOpqPoi9aAAARbCNAVD/aj6IAWqBDWtgqPM8Vb/ct4OjyNvZbVwcQQyUoqwXZxsFAEC9sF0zcIIYABjSbPABfeE5g74QxJCcFw0AGBxec2HoEcQAAKTj855AhO8hBgAAICl7iAEAKI3DiIEyCeIak+1FI9v6AgAAg0cQQw98rggY6rxhCAADJ4j7KduGSLb1LYPbmGrL9l3PHkMAQF8JYgAYgGwh7ggaAOqJs0wDAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSEsQAAACkJIgBAABISRADAACQkiAGAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAICVBDAAAQEqCGAAAgJQEMQAAACkJYgAAAFISxAAAAKQkiAEAAEhJEAMAAJCSIAYAACAlQQwAAEBKghgAAICUBDEAAAApCWIAAABSaixr4M7OzoiIePXVVwe0nCP/+1o1ptOr3uZnXOPW0pjGzTlupnU1bn2Mm2ldjVvf42ZaV+MOzriZ1vXdxu3PMt7qzt5UOt/rGgXZs2dPTJw4sYyhAQAASGD37t0xYcKEXi8vLYiPHDkSe/fujaampqhUKmVMAQAAgDrU2dkZBw4ciPHjx0dDQ++fFC4tiAEAAKBMTqoFAABASoIYAACAlAQxAPCefMIKgHokiAHgbZ566qlYtGhRTJo0KUaMGBGnnnpqXHrppfHss8+WPbVB8dxzz0WlUom77rorIiL2798fF154YTz22GPlTgwACiCIAeDP2tvbY/r06fHCCy/E17/+9XjwwQdjxYoV0dHREWeccUbs3Lmz7CkW7gMf+ED87Gc/izlz5kRExI4dO+Kee+6JI0eOlDwzAKg+Z5kGgIj46U9/Gp/4xCfi8ssvj5tvvrnbZS+99FKcfvrpMW7cuNi2bVs5EyxJR0dHzJo1K7Zs2RKtra1lTwcAqsoeYgCIiDVr1sSoUaPi+uuvf8dlJ554Ytx0003x6U9/Og4dOhSHDx+O22+/PaZOnRojRoyISZMmxYoVK+KNN97o+p3FixfHOeecE3fccUeceuqpMWLEiJg5c2Y888wzsWnTppg6dWocd9xxceaZZ8aOHTsG/Hutra3vCNaOjo6oVCrR0dERERF33XVXNDY2xpNPPhnTp0+P4cOHx0knnRQ33nhj1++8/ZDpt2I4ImLWrFnR2toa7e3tUalU4plnnuk21vr162PYsGGxe/fuft4DADD4GsueAACUrbOzMzZv3hznnXdeHHfccT1e5/zzz+/69yWXXBLr1q2LFStWxMc+9rHYvn17fPWrX43t27fHQw89FJVKJSIinnjiidi7d2/cdNNN8frrr8eyZcuira0tKpVKXHvttTFy5MhYunRpLFy4MH75y192Lb+/v3c0jhw5Eueff34sX748rrvuuli7dm1cffXVMXXq1Dj77LO7XfcjH/lItLe3x2WXXRbt7e3R2toa48ePj6uuuiruueeeWL16ddd177777vjkJz8ZEydO7NN8AKBMghiA9F5++eV44403YvLkye953V/96lexdu3auOGGG2LFihUREXHWWWfF+PHjY9GiRfHggw9GW1tbREQcOHAgfvjDH8Zpp50WERFbt26Nb3/72/HII4/E7NmzIyLi17/+dVx11VWxf//+GDVq1IB+72h0dnbGqlWr4uKLL46IiJkzZ8aPf/zj2LRp0zuCuLm5OaZMmRIREVOmTOn697x582L9+vVx7bXXRqVSiT179sSjjz4a69evP+p5AMBQ4JBpANJrbHzz/eHDhw+/53W3bt0aERELFizo9vP58+fHsGHDug5PjogYPXp0V9RGRIwbNy4iIs4888yun40ZMyYi3jyb80B/72hNnz6969/HHntsnHjiiXHo0KGj/v2LL744nnvuua4zT69bty6amppi3rx5fZ4LAJRJEAOQ3ujRo6OpqSmef/75Xq9z6NCheOWVV+KPf/xjRES8//3v73Z5Y2NjjB07tlugNjc397iskSNHvut8+vt7R+svDwtvaGjo01mkZ8+eHZMnT45169ZFxJtBPH/+/Bg+fHhV5gcAg0UQA0BEnH322bFly5ZuJ8Z6u+985zsxduzYrv/v27ev2+V/+tOf4uWXX+52ncFUqVTesYf74MGDhY21ePHi2LBhQ2zbti127doVn/vc5woZCwCKJIgBICKuvPLK+MMf/hArV658x2X79u2LG2+8MaZMmdJ1WPC9997b7To/+MEP4vDhw/HRj350UOb7l5qbm2PPnj3dfvb4448PeLnDhg3r8eeLFy+O/fv3x5VXXhktLS3dDucGgFrhpFoAEBHTpk2L1atXx8qVK+Ppp5+Oiy66KMaOHRu/+MUvYs2aNfH666/Hww8/HC0tLXHRRRfFqlWr4rXXXouPf/zjsWPHjrjmmmti1qxZcc4555Qy/7lz58bGjRtj+fLlcd5558Vjjz3WdUjzQLx1wq77778/Ro8eHR/+8IcjImLSpElx1llnxebNm+Mb3/jGgMcBgDLYQwwAf/alL30pHnjggYiI+MIXvhBtbW1x6623xty5c2PHjh3R0tISERFr166Nr3zlK/G9730v2traor29Pa644op44IEHoqGhnJfWJUuWxBe/+MX4/ve/H21tbfHEE0/EfffdN+DlfuhDH4oFCxbEbbfdFgsXLux22Zw5c2LYsGGxaNGiAY8DAGWodHZ2dpY9CQCg9px77rkxfPjw2LBhQ9lTAYB+ccg0ANAnq1evjl27dsXmzZur8jllACiLIAYA+mTjxo3xm9/8JtasWRMzZswoezoA0G8OmQYAACAlJ9UCAAAgJUEMAABASoIYAACAlAQxAAAAKQliAAAAUhLEAAAApCSIAQAASEkQAwAAkJIgBgAAIKX/A8mzryroF5jOAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gs = np.dot(relabus.T, G).T\n",
    "# 排序数据并选择绘制部分\n",
    "gs = gs / np.sum(gs, axis=0, keepdims=True)\n",
    "sorted_inds_g = np.argsort(np.mean(gs, axis=1))[::-1]\n",
    "gs = gs[sorted_inds_g][:, :30]\n",
    "\n",
    "colors = []\n",
    "colors.extend(plt.cm.tab20(np.linspace(0, 1, 20)))\n",
    "colors.extend(plt.cm.tab20b(np.linspace(0, 1, 20)))\n",
    "colors.extend(plt.cm.tab20c(np.linspace(0, 1, 20)))\n",
    "colors.extend(plt.cm.Set1(np.linspace(0, 1, 9)))\n",
    "colors.extend(plt.cm.Paired(np.linspace(0, 1, 12)))\n",
    "colors.extend(plt.cm.Set2(np.linspace(0, 1, 8)))\n",
    "colors.extend(plt.cm.Accent(np.linspace(0, 1, 8)))\n",
    "colors.extend(plt.cm.Set3(np.linspace(0, 1, 12)))\n",
    "colors = colors[:100]\n",
    "\n",
    "# 创建图形和坐标轴\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_position([0.1, 0.1, 0.8, 0.8])\n",
    "\n",
    "# 绘制堆叠柱状图\n",
    "x = np.arange(30)  # x轴上的点\n",
    "bottom = np.zeros(30)  # 底部值的初始化\n",
    "bh = []\n",
    "for i in range(gs.shape[0]):\n",
    "    b = ax.bar(x, gs[i], bottom=bottom, color=colors[i], edgecolor='none')\n",
    "    bh.append(b)\n",
    "    bottom += gs[i]  # 更新底部值，使得下一个柱子在上一个柱子的顶部\n",
    "\n",
    "ax.set_xlim([-1, 30])\n",
    "ax.set_ylim([0, 1])\n",
    "ax.set_yticks([])\n",
    "ax.set_xticks([])\n",
    "ax.set_ylabel('g', fontsize=12, fontname='Arial')\n",
    "ax.set_xlabel('Community', fontsize=12, fontname='Arial')\n",
    "ax.tick_params(width=0.5, length=0.0)\n",
    "ax.minorticks_off()\n",
    "ax.tick_params(direction='out')\n",
    "fig.set_size_inches(10, 5) \n",
    "fig.tight_layout() \n",
    "\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\lenovo\\AppData\\Local\\Temp\\ipykernel_13460\\3426175745.py:18: RuntimeWarning: invalid value encountered in divide\n",
      "  return temp / np.sum(temp, axis=0, keepdims=True)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "plt.rcParams['font.family'] = 'Arial'\n",
    "\n",
    "def shuffle_matrix_rows(matrix):\n",
    "    # 对矩阵进行复制，防止原矩阵被修改\n",
    "    shuffled_matrix = np.copy(matrix)\n",
    "    # 对于矩阵中的每一行进行随机排列\n",
    "    for row in shuffled_matrix:\n",
    "        np.random.shuffle(row)\n",
    "    return shuffled_matrix\n",
    "\n",
    "\n",
    "def ConvertToG(G, relabus):\n",
    "    temp = np.dot(relabus.T, G).T\n",
    "    return temp / np.sum(temp, axis=0, keepdims=True)\n",
    "\n",
    "#一种相似度计算方法\n",
    "def CalculateBC(xthis, ythis):\n",
    "    return 1 - np.sum(np.minimum(xthis, ythis))\n",
    "\n",
    "# 计算每列与均值的相似度度量\n",
    "def CalculateDs(xs):\n",
    "    mm = xs.shape[1]\n",
    "    xm = np.mean(xs, axis=1, keepdims=True)\n",
    "    ds = np.zeros(mm)\n",
    "    for iii in range(mm):\n",
    "        ds[iii] = CalculateBC(xs[:, iii], xm.flatten())\n",
    "    return ds\n",
    "\n",
    "dg_orig = CalculateDs(ConvertToG(G, relabus))\n",
    "dg_m_orig = np.mean(dg_orig)\n",
    "\n",
    "# 模拟\n",
    "ns = int(1e4)\n",
    "\n",
    "dg_m_shfls = np.zeros(ns)\n",
    "for rep in range(ns):\n",
    "    # 随机按行打乱矩阵G\n",
    "    G_shuffled = G[np.random.permutation(N), :]\n",
    "    dg_m_shfls[rep] = np.mean(CalculateDs(ConvertToG(G_shuffled, relabus)))\n",
    "\n",
    "\n",
    "dg_m_shfls_row = np.zeros(ns)\n",
    "for rep in range(ns):\n",
    "    # 打乱每一行G的分布\n",
    "    #G_shuffled = G[np.random.permutation(N), :]\n",
    "    G_shuffled = shuffle_matrix_rows(G)\n",
    "    dg_m_shfls_row[rep] = np.mean(CalculateDs(ConvertToG(G_shuffled, relabus)))\n",
    "\n",
    "#fig, ax = plt.subplots()\n",
    "\n",
    "dg_m_shfls_col = np.zeros(ns)\n",
    "for rep in range(ns):\n",
    "    # 打乱每一列G的分布\n",
    "    #G_shuffled = G[np.random.permutation(N), :]\n",
    "    G_shuffled = shuffle_matrix_rows(G.T).T\n",
    "    dg_m_shfls_col[rep] = np.mean(CalculateDs(ConvertToG(G_shuffled, relabus)))\n",
    "\n",
    "#fig, ax = plt.subplots()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhmklEQVR4nO3deZBU9b3w4W8j+6pcBGQRjUkU4ygxLgS9QK4SEY2IYBSjgmBckpSWVjSVqMW4xaSMWyVSrldxQwkkiiBRUceI4oavOgY1LheCKCpqxCA7/f5hmDgO6HRP/7pnep6nyoJezunvpM7l8uF3+pxMNpvNBgAAAFBwLUo9AAAAAJQr0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACTSslQfvHHjxnj77bejU6dOkclkSjUGAAAA5CybzcYnn3wSvXr1ihYttryeXbLofvvtt6Nv376l+ngajW9HxHMRsWdE/L8Sz0I56lbRLZaPXh7dZnSL5dXLSz0OAABlZsmSJdGnT58tvl6y6O7UqVNEfDZg586dSzUG/3bJJaX53GXLWsTNN0eMH//X6NlzY0RE9Fj2fEy4eUj87/hH492eAwr2Wb/85X9+/8Wf9/OvUV4efvnhGPWXUXH99dfH//T/n1KPAwBAmVixYkX07du3pm23pGTRvemU8s6dO4vuRqBNm9J8buvWm37tWDNDu9Ydo/O/f23TpnDHxucPsy/+vA7B8tWufbuaX/1ZAwBAoX3V16VdSA0AAAASEd0AAACQSMlOL4fGpLLyyx8DAADkw0o3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAIm0LPUAUCyVlaWeAAAAaG6sdAMAAEAiVrqbMSu/AAAAaVnpBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACCRlqUeACIiXnvttVi+fHVERPzro9dqnvu/5a3qvY+KiookswEAAOTLSjcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASKRlqQeAQqmuri7IfioqKgqyHwAAACvdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIpGWpB4DGqLJy878HAADIhZVuAAAASMRKdzNixRYAAKC4rHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACTSstQDAOWjsqoy/22H5r8tAAA0Vla6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiLqQG1NGQC6IBAAD/YaUbAAAAErHSDY1cvqvObsEFAAClZ6UbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJuGUYlKl8bzUGAAAUjpVuAAAASMRKN9Ao5LsyXzk0v+0AAKAYrHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACARtwyjQaqrqxu0/UcftY2IbxRmGAAAgEbGSjcAAAAkIroBAAAgEaeXQ5FUVlWWegQAAKDIrHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIi1LPQA0NtXV1bUeT536Ul77GTt2bCHGAQAAmjAr3QAAAJCI6AYAAIBERDcAAAAk4jvdQJNWWVX5pa+/8+47ERHxp4V/ivkfzP/PdkO/fDsAACgEK90AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIJGWpR4AmpLpy6fX+72vVr2acBIAAKApsNINAAAAiYhuAAAASMTp5WWusrLUEwAAADRfVroBAAAgEdENAAAAiYhuAAAASMR3uoFmqbKqMv9th+a/LQAAzYuVbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAi9Y7uyy+/PN59992UswAAAEBZqXd0n3XWWbF48eKax9lsNn7yk5/EW2+9lWQwAAAAaOpa1veN2Wy21uONGzfGtddeGyeeeGL06dOn4IMBNFaVVZX5bTc0v+0AAGi6GvSd7i+GOAAAAPAfLqQGAAAAiYhuAAAASCSn6M5kMvV6DgAAAMjhQmoREQMHDqzz3F577VXnuUwmE+vXr89/KgAAACgD9Y7uSZMmpZwDAAAAyo7oBiiSfG81FuF2YwAATVVOp5d/3pIlS2LZsmWRyWSiT58+0bNnz0LOBQAAAE1eThdSW7t2bVxyySXRr1+/2GGHHWLgwIGx7777Ru/evWPnnXeOq666KjZs2JBqVgAAAGhS6r3SvXr16hg2bFg8/vjjMXDgwDjhhBOiZ8+ekc1m46233oq5c+fGGWecEbNmzYr77rsvWrVqlXJuAAAAaPTqHd2XXnppPPvss3H33XfHYYcdVuf1iy++OGbPnh0//OEP4w9/+EOcccYZBR0UAAAAmpp6n17+xz/+Mc4888zNBvcmhxxySJx22mkxderUggwHAAAATVm9o/vNN9+MwYMHf+X7hgwZEq+88kqDhgIAAIByUO/oXrVqVXTp0uUr37f11lvHypUrGzQUAAAAlIN6R3c2m40WLb767ZlMpkEDAQAAQLnI6ZZhghoAAADqr95XL4+IOPzww6NNmzZf+p41a9Y0aCAAAAAoF/WO7uOPP95KNwAAAOSg3tF98803R0TE008/HYsXL46ddtop9txzz1RzAQAAQJNX7+j+5z//GYceemjMnz+/5rlBgwbFHXfcEX379k0yHAAAADRl9b6Q2rnnnhvPPfdcnH/++TF79uy47LLL4pVXXomTTz455XwAAADQZNV7pfvee++NSy65JE4//fSIiBg+fHj07t07jjnmmFi5cmV06NAh2ZAAAADQFNV7pXvZsmXxne98p9ZzQ4cOjQ0bNsQ//vGPgg8GAAAATV29o3vdunXRunXrWs917do1IiJWr15d2KkAAACgDNQ7ur9MNpstxG4AAACgrOQU3Vu6T7f7dwMAAEBd9b6QWkTEqaeeGp07d655vGmF+6STTopOnTrVPJ/JZOKhhx4q0IgAVFZV5rfd0Py2AwCgMOod3YMHD45MJlPnVPIhQ4ZERO1TzJ1uDgAAADlEd1VVVcIxAAAAoPwU5EJqAAAAQF2iGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCItSz0AlKvq6uoG76OioqIAkwAAAKVipRsAAAASsdJNs7Swx+R6v3f68vcSTgIAAJQz0Q1QxiqrKvPbbmh+2wEAUJvTywEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERalnoASqe6urrUIwAAAJQ1K90AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQSMtSDwD5WthjcqlHAAAA+FJWugEAACAR0Q0AAACJiG4AAABIRHQDAABAIi6kVmYqK0s9AQAAAJtY6QYAAIBERDcAAAAkIroBAAAgEd/pBqCOyqrK/Lcdmv+2AADlxko3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEjEhdSAJqm6urpe71u1ZlVERCxZsiSWv7e81msVFRUFnwsAAD5PdANQUPle+dxVzwGAcuT0cgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAIm0LPUAQPNSXV1d6hFqFGqWioqKguwHAIDyI7rhKyx8uXutx7v2f69EkwAAAE2N08sBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBE3KcbgEahsqoy/22H5r8tAEBKohsaserq6oLsp6KioiD7AQAAciO6ARrIP44AALAlvtMNAAAAiYhuAAAASER0AwAAQCKiGwAAABJxITWgXgp1sTAAAGhOrHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIpGWpBwCApmLq1KkF2c/YsWMLsh8AoPET3QBQZIWId+EOAE2D08sBAAAgESvd0AxUV1eXegQAAGiWrHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACTi6uWQo4Uvd6/5/a793yvhJECpTF8+Pe9tx3QbU5AZCnGv7wj3+waA1EQ3AE1eZVVlftsNzW87AID6Et0ANAubu1/91HcKs1oMALAloptG4f/+64/xzrZLPnvQ4v2a5xZuu20JpwIof05TB4C0RDc0wOe/3x3hO97Q1DTku9kAAPXh6uUAAACQiOgGAACARJxeDgA0Cr5fDkA5Et1NUGXl5n8PAABA4yK6m6DP3/Zm6tSXvvDabsUeBwAAgC3wnW4AAABIRHQDAABAIqIbAAAAEvGdbgCgwQp15XEAKDdWugEAACAR0Q0AAACJOL0coJH4/O0A81VRUVGASQAAKBQr3QAAAJCIlW4AGrVCnAEAAFAqVroBAAAgESvdTdz06buVegQAaFQKdfuysWPHFmQ/ADRvVroBAAAgESvdUEALX+5e8/td+79Xwkmg9HwXGwBAdANAUU1fPj2v7cZ0G1PgSQCAYhDdANRhlRoAoDB8pxsAAAASEd0AAACQiOgGAACAREQ3AAAAJOJCagAAmzF16tQG72Ps2LEFmASApsxKNwAAACRipRugjLjVFwBA42KlGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABJpWeoBAICvNn359Ly3HdNtTAEnAQByIboBoMzlG+xiHQAazunlAAAAkIjoBgAAgERENwAAACQiugEAACARF1KjYBb2mJzzNqtb9I2IXxV+GABoBKZOnVqQ/YwdO7Yg+wGg+Kx0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEbcMAwBo5Apx6zG3HQMoDdHdBFRWlnoCAAAA8uH0cgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACARFy9vJFyxfKmb+HL3Ws93rX/eyWaBAAKc9uxCLceA8iVlW4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERcvRwA2Kzpy6fntd2YbmMKPAkANF2iu4hyuVVHdfVuCScBAACgGEQ3tSzsMbnUIwAAAJQN0Q0AFFS+p6VHODUdgPLjQmoAAACQiJVuAADqLZdr1HyZsWPHFmQ/AI2d6AYAoOgKEe/CHWgKnF4OAAAAiVjpbiSmT3eLMAAAgHIjugGARiPfK5+76jkAjZXTywEAACARK90AADRrLuoGpCS6AQBokgp1+zKAlJxeDgAAAIlY6YYiWfhy95rf79r/vRJOAlB+8r0AW4SLsAGQlugGAIAGKtSp7r4bDuVHdAMA5MHtzQCoD9ENADRrDTk1HQC+iuiGRsj3vwGgeXKaOpQf0V2mFvaYXOoRAIDNaGor606HB2gY0Q2NwOdXtgGgMfHddYCGEd31VKhTfQAAILVC/N3VKepQGKK7hKZP363UIwAAJGGFnE18T53mTnQDANBolOI770J/85zpCYUhuqEEcvkO95e915XNAQDqz6o7pVCy6M5msxERsWLFilKNkJNPP/204Ptct+5fX/r6K91vyH/na/PftJg2rFtX8+uGtZ8NvW7duljx7183PcfmrVu1rtQjNHrr16//7Ne162PdRv97AVDX1CXFX9Ed+V8ji/6ZTd2NN95Y6hFqFGKWI488sgCTUEqbWnZT225JJvtV70jkrbfeir59+5biowEAAKAglixZEn369Nni6yWL7o0bN8bbb78dnTp1ikwmU4oRAAAAIC/ZbDY++eST6NWrV7Ro0WKL7ytZdAMAAEC523KOAwAAAA0iumn2Vq9eHRMnToytt946tttuu7jsssu+cpt58+bF1772tSJMBw2XyzE+e/bsGDBgQHTs2DF23333mDlzZhEnhfzkcozffvvt8c1vfjPatWsXgwYNiqeffrqIk0J+8vm7yqJFi6Jjx45RVVWVfkBooFyO8ZEjR0Ymk6n136xZs4o4be7cMoxm76yzzopnn302Hn744Vi8eHGMGzcu+vXrF2PGbP6endXV1TFmzJho27ZtkSeF/NT3GH/xxRfjiCOOiEsvvTRGjBgR999/f4wZMyaeeeaZ2GOPPUo0PXy1+h7jjz32WEycODFuuOGGGDRoUEyePDkOPvjgWLx4cXTs2LFE08NXy/XvKhERp556aqxcubKIU0L+cjnGFy5cGLfddlsccMABNc9ts802xRw3d1loxv71r39l27Ztm33kkUdqnrvwwguzQ4YM2ez7r7nmmmzHjh2zu+++e7Zfv35FmREaIpdj/Be/+EV2+PDhtZ77/ve/n/3Vr36VeErIXy7H+LRp07IXXXRRzeOPP/44GxHZp556qgiTQn5y/btKNpvN3nbbbdn99tsvGxG1toPGKJdjfPXq1dmtttoq++qrrxZvwAJwejnN2gsvvBDr1q2LQYMG1Ty3//77x1NPPRUbN26s8/45c+bElClT4owzzijmmJC3XI7xcePGxW9+85s6+/j444+Tzwn5yuUYP/LII+Occ86JiIhVq1bFFVdcEd27d49dd921qDNDLnL9u8oHH3wQZ599dlx77bXFHBPylssx/uqrr0Ymk2lyX/MU3TRr77zzTnTr1i1at25d81yPHj1i9erV8cEHH9R5/9133x1HHHFEMUeEBsnlGO/fv3+t08j/9re/xUMPPVTr9C1obHL9czwi4qGHHoqOHTvG+eefH1deeaVTy2nUcj3GzzzzzBg3blx861vfKuaYkLdcjvGXX345unTpEscdd1xst912sc8++8ScOXOKPXLORDfN2qeffhpt2rSp9dymx2vWrCnFSFBQ+R7jy5cvj9GjR8d+++0XI0eOTDojNEQ+x/huu+0WCxYsiAsuuCDGjx8fTz75ZPI5IV+5HONz586NefPmxXnnnVe0+aChcjnGX3nllfj000/joIMOir/85S8xYsSI+MEPfhDPPvts0ebNhwup0ay1bdu2zv8xb3rcvn37UowEBZXPMf7uu+/GsGHDYuPGjTF9+vRo0cK/z9J45XOM9+jRI3r06BEDBgyIJ598Mq655poYOHBg8lkhH/U9xletWhUnn3xyTJ48Odq1a1fUGaEhcvlz/LzzzovTTjut5sJpe+yxRyxYsCCuu+662GuvvYozcB78TYpmrXfv3rF8+fJYv359zXPLli2Ldu3axdZbb126waBAcj3Gly5dGoMHD441a9ZEVVVVbLvttkWcFnKXyzH+zDPPxHPPPVfruV133TWWL19ejFEhL/U9xp9++ul48803Y/To0dGxY8ear00cfPDBccoppxR7bKi3XP4cb9GiRZ0rlffv3z+WLl1ajFHzJrpp1gYMGBCtWrWqdWrhvHnzYu+997a6R1nI5RhfuXJlDB8+PFq0aBGPPvpo9OrVq9jjQs5yOcZvvPHG+OUvf1nruQULFkT//v2LMivko77H+D777BOvvfZaPP/88zX/RUTccMMNccEFFxR7bKi3XP4cHz9+fEyYMKHWc88//3zssssuRZk1X04vp1lr3759jBs3Lk455ZS46aabYunSpfG73/0ubrrppoj47F/ZunTp4jQtmqxcjvFf//rX8cYbb0RVVVXNaxER7dq1iy5dupTqR4AvlcsxftJJJ8W+++4bV111VYwYMSJuu+22ePrpp+OWW24p8U8BW5bLMf71r3+9zva9e/eO7t27F3tsqLdcjvHDDjssjj766Bg6dGgMGjQo7rjjjpg3b15cd911Jf4pvkKp71kGpbZy5crs8ccfn+3QoUO2V69e2SuuuKLmtYjI3nTTTXW2uemmm9ynmyajvsf4zjvvnI2IOv+NGzeuJHNDfeXy5/i9996braioyLZt2za71157ZR9//PHiDww5yufvKptec59umoJcjvHrr78++41vfCPbpk2b7J577pl99NFHiz9wjjLZbDZbot4HAACAsuZLqwAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwDISzabLfUIANDoiW4AKKBhw4ZF165dY+3atVt8T0VFRQwePDiv/S9atCgymUzcfPPNBd9mhx12iPHjx0dERFVVVWQymaiqqoqIiMrKyshkMjXvffzxx+OQQw7JcXoAaH5alnoAACgnEyZMiLlz58Z9990Xhx9+eJ3Xn3vuuXjppZdiypQpee1/u+22i/nz58dOO+3UwEnr+vOf/xydO3fe7GsnnnhiDB8+vObx9ddfHwsXLiz4DABQbkQ3ABTQqFGjYptttonbb799s9E9ZcqU6Ny5c4wZMyav/bdp0yYGDhzYwCk379vf/vYWX+vTp0/06dMnyecCQDlzejkAFFDbtm3jmGOOiVmzZsWKFStqvbZu3bqYOnVqjB07Nj799NP46U9/Gv369YvWrVtH165dY9SoUbFo0aKa9w8dOjSOPfbYGDNmTHTo0CGGDRu22VPF//rXv8ZBBx0U22yzTbRu3Tp23HHHqKysjI0bN9b6/KVLl8ahhx4a7dq1i759+8akSZNiw4YNNa9//vTyL/r86eXjx4+PKVOmxOLFi2tm2XvvvWPQoEF1tjvwwANj2LBhOf6vCADlQ3QDQIFNmDAhVq9eHTNmzKj1/Jw5c+L999+PiRMnxiGHHBIPPPBA/Pa3v40HHnggKisr46GHHopTTjml1jZ33XVXdOrUKWbOnBlnn312nc964YUX4oADDohu3brFXXfdFffee2/893//d5x//vkxbdq0Wu+dNGlS9OjRI+6+++6YMGFCXHzxxfHzn/8855/vvPPOixEjRkTPnj1j/vz5ccghh8TEiRNj/vz58frrr9e8b8mSJfHII4/ECSeckPNnAEC5cHo5ABTYnnvuGQMGDIjbb7+9VnBOmTIldt999+jdu3d06NAhLrvssth///0j4rNV7ddffz2uu+66Wvtq3bp1XHPNNdGmTZuIiFor4RERL774YgwbNixuvfXWaNHis39LHzZsWMycOTOqqqri6KOPrnnv8OHD48Ybb4yIiIMOOihWrFgRkydPjvPOOy+6du1a759vp512im233bbWqe5jx46NM888M2699dY4//zzIyLi1ltvjU6dOsWoUaPqvW8AKDdWugEggQkTJsQjjzwSb7/9dkREfPjhhzFr1qyYOHFi9OrVKx5++OHYb7/9YtGiRfHggw/G73//+3j88cdjzZo1tfbTv3//muDenOOOOy7uu+++WLt2bbz44osxY8aMmDRpUqxfv77Ovo466qhaj4844ohYu3ZtPPnkkw3+ebt06RKjR4+O2267rea5KVOmxFFHHRXt2rVr8P4BoKkS3QCQwI9+9KNo1apV3HnnnRERceedd0Ymk4ljjz02IiJuv/326NevX+y4445x9NFHxz333BPt27evs5+OHTt+6eesWrUqTjzxxOjSpUsMGDAgzj777Fi0aFG0atWqzn20e/bsWetx9+7dIyLio48+yvvn/LyJEyfGm2++GfPmzYv58+fH3//+d6eWA9DsiW4ASKBr165x+OGHxx133BEREbfcckuMGjUqunbtGvPmzYvjjz8+Ro8eHW+99VZ88MEHMXfu3Pjud7+b8+ecfvrpMX369Jg2bVp88skn8cYbb8Stt94arVq1qvPeDz/8sNbjZcuWRcR/4ruhhgwZEjvttFNMmzYtpk2bFrvsskuyK60DQFMhugEgkQkTJsSCBQvi0UcfjaeeeiomTpwYERFPPPFEbNy4MSorK6N3794REbFhw4Z48MEHIyLqXHX8y8ybNy++973vxciRI6NDhw4REbFgwYJ4//336+xn9uzZtR7feeed0b59+9h3331z/tm22mqrOs9lMpk44YQT4u67746ZM2du8UroANCcuJAaACRy4IEHxvbbbx8nnXRS7LjjjnHAAQdERMQ+++wTERE/+9nPYsKECfHhhx/G1VdfHS+88EJERKxcuTI6depUr8/YZ599Ytq0aXHNNddE//7944UXXoiLLrooMplMrFy5stZ7Z8yYEb17945hw4bF/fffH9dee21ceOGF0blz55x/tq233jrefffdmDNnTgwYMCC22267iPjsdmKTJk2KiM++bw4AzZ2VbgBIpEWLFjF+/Pia7zZvus/10KFD4+qrr44nnngiDj744DjzzDNj++23jz/96U8REfHYY4/V+zMuv/zyGDVqVJx77rlx6KGHxg033BDnnntu/PjHP4758+fXug/3VVddFc8880yMGDEipk2bFldeeWWcc845ef1sJ5xwQuywww4xcuTImDJlSs3zvXv3jj322CMOOuig6NWrV177BoByksl+8SorAAB5Wrp0afTr1y9mzJgRI0eOLPU4AFByohsAaLDnn38+7rnnnpgxY0asX78+XnrppZr7hgNAc+b/GwIADbZ69eq4/PLLY/369XHnnXcKbgD4NyvdAAAAkIh/hgYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBE/j+fE3KLDpmZLAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 可视化\n",
    "fig, ax = plt.subplots()\n",
    "dg_m_shfls = dg_m_shfls[~np.isnan(dg_m_shfls)]\n",
    "xthis = dg_m_shfls\n",
    "xmin = np.min(xthis)\n",
    "xmax = np.max(xthis)\n",
    "bins = np.linspace(xmin * 0.9, xmax * 1.1, 50)\n",
    "h = ax.hist(xthis, bins=bins, density=True, edgecolor='none', facecolor='grey',alpha=0.7, zorder=1)\n",
    "\n",
    "# 添加原始 dg 的直线\n",
    "ymax = np.max(h[0])\n",
    "ax.plot([dg_m_orig, dg_m_orig], [0, ymax * 1.1], '-', linewidth=1, color='red')\n",
    "ax.plot([np.mean(dg_m_shfls), np.mean(dg_m_shfls)], [0, ymax * 1.1], '-', linewidth=1, color='grey', zorder=5)\n",
    "\n",
    "xthis = dg_m_shfls_row\n",
    "xmin = np.min(xthis)\n",
    "xmax = np.max(xthis)\n",
    "bins = np.linspace(xmin * 0.9, xmax * 1.1, 50)\n",
    "h = ax.hist(xthis, bins=bins, density=True, facecolor='blue',alpha=0.5, zorder=1)\n",
    "ax.plot([np.mean(dg_m_shfls_row), np.mean(dg_m_shfls_row)], [0, ymax * 1.1], '-', linewidth=1, color='blue', zorder=5)\n",
    "\n",
    "dg_m_shfls_col = dg_m_shfls_col[~np.isnan(dg_m_shfls_col)]\n",
    "xthis = dg_m_shfls_col\n",
    "xmin = np.min(xthis)\n",
    "xmax = np.max(xthis)\n",
    "bins = np.linspace(xmin * 0.9, xmax * 1.1, 50)\n",
    "h = ax.hist(xthis, bins=bins, density=True, facecolor='green',alpha=0.5, zorder=1)\n",
    "ax.plot([np.mean(dg_m_shfls_col), np.mean(dg_m_shfls_col)], [0, ymax * 1.1], '-', linewidth=1, color='green', zorder=5)\n",
    "\n",
    "# 设置图形属性\n",
    "ax.set_ylim([0, ymax * 1.1])\n",
    "ax.set_yticks([])\n",
    "#ax.set_xlim([0.025, 0.175])\n",
    "#ax.set_xticks([0.05, 0.1, 0.15])\n",
    "ax.set_ylabel('PDF', fontsize=12, fontname='Arial')\n",
    "ax.set_xlabel('Variability', fontsize=12, fontname='Arial')\n",
    "ax.tick_params(width=0.5, length=0.0)\n",
    "ax.minorticks_off()\n",
    "ax.tick_params(direction='out')\n",
    "fig.set_size_inches(10, 5) \n",
    "fig.tight_layout() \n",
    "\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07839236474495727"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(dg_m_shfls_row)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "function_xys",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}