Hadza_Nepal_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import os\n",
    "\n",
    "# Load the .mat files\n",
    "data = sio.loadmat('Hadza_converted.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomes = np.array([gen_name[0] for gen_name in data['genomes'].flatten()])\n",
    "usamples = np.array([sample_name[0] for sample_name in data['usamples'].flatten()])\n",
    "relabus = data['relabus']\n",
    "kegg_per_genome = data['keggs_per_genome']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "inds2del = np.all(relabus <= 1e-6, axis=1)  # Identify rows with all values <= 1e-6\n",
    "relabus = relabus[~inds2del, :]  # Remove these rows from relabus\n",
    "G = kegg_per_genome > 0\n",
    "G = G[~inds2del, :]\n",
    "\n",
    "# samples filter\n",
    "coverage_genome = np.any(kegg_per_genome>0, axis=1)\n",
    "coverage = np.sum(relabus[coverage_genome,:], axis=0)\n",
    "inds2del = coverage < 0.8\n",
    "relabus = relabus[:,~inds2del]\n",
    "usamples = usamples[~inds2del]\n",
    "\n",
    "\n",
    "relabus = relabus[coverage_genome,:]\n",
    "G = G[coverage_genome,:]\n",
    "\n",
    "relabus = np.nan_to_num(relabus)  # Replace NaNs with 0\n",
    "relabus = relabus / np.sum(relabus, axis=0, keepdims=True)  # Normalize relabus\n",
    "\n",
    "keggs = np.array([i for i in range(0,int(data['max_ko_num']))]).T\n",
    "inds2del = np.sum(G, axis=0) == 0  # Identify columns in G where the sum is 0\n",
    "G = G[:, ~inds2del]  # Remove these columns from G\n",
    "keggs = keggs[~inds2del]  # Remove corresponding elements from keggs\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_m = np.zeros(len(usamples)).astype(int)\n",
    "for i, name in enumerate(usamples):\n",
    "    samp_name = name.split('_')[0]\n",
    "    if samp_name == 'Nepal':\n",
    "        group_m[i]=1\n",
    "relabus = relabus[:,group_m>0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3722, 57)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relabus.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KEGGs\\M00001.txt\n",
      "KEGGs\\M00004.txt\n",
      "KEGGs\\M00009.txt\n",
      "KEGGs\\M00022.txt\n",
      "KEGGs\\M00023.txt\n",
      "KEGGs\\M00024.txt\n",
      "KEGGs\\M00025.txt\n",
      "KEGGs\\M00049.txt\n",
      "KEGGs\\M00050.txt\n",
      "KEGGs\\M00051.txt\n",
      "KEGGs\\M00052.txt\n",
      "KEGGs\\M00115.txt\n",
      "KEGGs\\M00120.txt\n",
      "KEGGs\\M00127.txt\n",
      "KEGGs\\M00144.txt\n",
      "KEGGs\\M00958.txt\n",
      "KEGGs\\M00959.txt\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "\n",
    "G_orig = G.copy()  # Make a copy of G\n",
    "files = os.listdir('KEGGs')  # List all files in the KEGGs directory\n",
    "ko_all = []\n",
    "ko_all_inds = []\n",
    "\n",
    "# Loop through the files starting from the third file (index 2 in Python)\n",
    "for filenum in range(0, len(files)):\n",
    "    filename = os.path.join('KEGGs', files[filenum])\n",
    "    ko2look = []\n",
    "\n",
    "    # Open the file and read line by line\n",
    "    with open(filename, 'r') as fid:\n",
    "        for line in fid:\n",
    "            sthis = line.strip().split()[0] \n",
    "            ko2look.append(int(sthis[1:]))  \n",
    "\n",
    "    # ko2look: KO from each .txt\n",
    "    \n",
    "    \n",
    "    Mthis = len(ko2look)\n",
    "    inds2look = np.zeros(Mthis, dtype=int)\n",
    "\n",
    "    # Find indices in keggs that match the ko2look values\n",
    "    for j in range(Mthis):\n",
    "        ind = np.where(keggs == ko2look[j])[0]\n",
    "        if ind.size > 0:\n",
    "            inds2look[j] = ind[0]\n",
    "\n",
    "    inds2look = inds2look[inds2look != 0]  # Remove zero entries\n",
    "\n",
    "\n",
    "    # Calculate the proportion of samples with G == 1\n",
    "    G = G_orig[:, inds2look]  # Select columns from G_orig\n",
    "    N, M = G.shape\n",
    "    ps = np.sum(G, axis=0) / N  # Calculate the proportion of samples with G == 1\n",
    "\n",
    "    if np.any(ps > 0.3):  # Check if any proportion is greater than 0.3\n",
    "        print(filename)\n",
    "\n",
    "    # Append the corresponding keggs values to ko_all\n",
    "    ko_all.append(keggs[inds2look[np.argmax(ps > 0.3)]])\n",
    "     # record keggs index for each keggs value\n",
    "    ko_all_inds.append([inds2look[np.argmax(ps > 0.3)]])\n",
    "   \n",
    "\n",
    "\n",
    "G = G_orig[:,np.array(ko_all_inds).flatten()]\n",
    "N, M = G.shape\n",
    "mx = np.mean(relabus, axis=1)\n",
    "\n",
    "\n",
    "dmx = np.zeros(M)\n",
    "for j in range(M):\n",
    "    dmx[j] = np.log10(np.mean(mx[G[:, j] == 1])) - np.log10(np.mean(mx[G[:, j] == 0]))\n",
    "dmx[dmx < -6] = -6\n",
    "\n",
    "ps = np.sum(G, axis=0) / N\n",
    "n1 = np.sum(G, axis=0)\n",
    "n1_orig = n1\n",
    "dm_origs = dmx \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 101)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 设置要绘制的柱状图数目\n",
    "mm2plot = 30\n",
    "mm = relabus.shape[1]\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_vis 的值\n",
    "N_vis = 100\n",
    "\n",
    "# 转置 relabus 并选择要绘制的部分\n",
    "xthis = relabus.T\n",
    "xthis = xthis[inds2plot, :][:, sorted_inds_x]\n",
    "\n",
    "# 合并第 N_vis+1 列及之后的所有列，并将结果存储在第 N_vis+1 列\n",
    "xthis[:, N_vis] = np.sum(xthis[:, N_vis:], axis=1)\n",
    "\n",
    "# 只保留前 N_vis+1 列\n",
    "xthis = xthis[:, :N_vis+1]\n",
    "\n",
    "xthis.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCkklEQVR4nO3dfbhXZZkv8Ae0BmV8O5UhomMTWFBOkDpsJBMpL5VIxU5doZFtUTgMZJPGlZ4pNUtMa0ZTsenFmcLAYx3UtGTjOZIYETXYRkF2aIkoim+NqDExVynnD4+K097KfrnXs14+n79qAff33gjs373Ws56n37Zt27YlAAAAoM/1z90AAAAA1JWhGwAAAIIYugEAACCIoRsAAACCGLoBAAAgiKEbAAAAghi6AQAAIIihGwAAAILsnCv4+eefT4888kjabbfdUr9+/XK1AQAAAN22bdu29Oyzz6bBgwen/v27fp6dbeh+5JFH0n777ZcrHgAAAHrtoYceSkOGDOnyx7MN3bvttltK6YUGd999917V2mOPPfqipS49/fTTpcuOzs2Z3cSvOWe2r7kc2U38/a7z15wz29dcjuwm/n7X+WvOme1rLkd2E3+/6/w1v1p2dzzzzDNpv/32e2m27Uq/bdu2bet1Wg8888wzaY899khPP/10r4fu6OXpr/ZblCu7iCX5ubKb+DXnzPY1lyO7ib/fdf6ac2b7msuR3cTf7zp/zTmzfc3lyG7i73edv+ZXy+6OHZ1pbaQGAAAAQQzdAAAAEMTQDQAAAEEM3QAAABDE0A0AAABBDN0AAAAQxNANAAAAQQzdAAAAEMTQDQAAAEEM3QAAABDE0A0AAABBDN0AAAAQxNANAAAAQQzdAAAAEMTQDQAAAEEM3QAAABDE0A0AAABBDN0AAAAQxNANAAAAQQzdAAAAEMTQDQAAAEEM3QAAABDE0A0AAABBDN0AAAAQxNANAAAAQQzdAAAAEMTQDQAAAEF2zt0AAAAAvOhP/3hq7hb6lKEbAACA0tjpzNG5W+hThm6AhjvvF+flbgEA4CXP3742PKP/uPCIlxi6ARru/L8dnLuFRlly6725WwAACmToBmi46LvJRd5JBgCq79v3jgjPmDYuPOIlhm6AhvvbDbeF1l8ZWh0AqJt/efivwjOmhSe8zNAN0HBn3fKn2IBTYssDAPXyzod/UEDK0QVkvMDQDdBwm8edkbsFAICXjB19eu4W+pShGwAAgNLY6YjVBaQUdyyZoRsAAIDSWP3EffEhB8ZHvMjQDQAAQGnctv7++JCx8REvMnQDQIHa9/ldaP0j07DQ+gBA9xi6AaBAz/xHW3BCS3B96Ny2bd/I3QJQE5982z/mbqFPGbphB2zuf17uFoCaOPc/NuduAQBK7c6NW8IzWg8Nj3iJoRt2wKU/j884Pz4CKIFHTl4YWn/Iw5eF1oeuPL809tWJlFLqf0R4BECfM3TDDpj2xzflbgGoiRPm7Blaf2Vodeha/32uKyDlnAIygNzO2py7g75l6AaABlhy6725W6DmnvvRweEZO50ZHsEOOO8XXruD7jB0A0CBzrrlT7EBp3R+eeke82NzU0pHelGm0XY6c3TuFijI/1j6/tiAv40tD0UzdANAgSZfNzBL7ueXPRQf4oNyoz1/+9rwjP7jwiPYAYNmR/+3fk9wfSiWoRug4e68887cLVAATyGJdt0enwnPmByewI6Y+d5+ofXn3hFaHgpn6AZouLGjT8/dAgXwFJJwXz4qPuO6jvgMXtMVF9yTuwWoFEM3AECNbNv2jdwtAPTK5vZN8SGtw+Mz/j9DN0DDrR1xfXCCZc3bO+20kaH1v/3tzq//+IB9Q3NTSumD4QlAGTxy8sLQ+kMeviy0PhTN0A3QcLetvz82YGznl5t65Eyu5fwX/uvfhGd88AvhEeyAn9zQxV+6PnLkpNDyVMAJc/YMrb8ytDoUz9ANO6D/8nfFh8R+RoIurRz45twtNIqVBdRVrp35KV6uow+hqgzdANAA73z4BwWkHF1ABpDb195+WWh9u9SzcM0t4Rmj0vjwjBcZugHIoqnnRv/h8U9kybVLfXMs3Tf2BsuR6R2h9Sm/OUe8LncL1NysDxe3yVkRDN2wA/50WQH3XGcXMIBAiTg3ulhX7PLV8IzWVMTTdF7L7h0TYwO6uJl13I3xr6rcdEJ4BECfM3QDNNxj184Lrb/PpCs7vT7zvf1Cc1NKae4d4RFQOkMXT4kNOKXzs7K/PvPu2NyUUjohPgKgrxm6Acjiigvuyd1Co/z4psXxIR+Lj6C8zjs2/p3+Lk7Eo2ALFiwIrT9+fHHv2lJOy2+O3nQ0pRNnTA3PeJGhG3ZA9NEYKTkeg3wunPyL0PqdP+dO6fn7RoTmppRS/3HhEd121ubcHVB3E6fl2dDOvgEAnavF0N0y68bcLVBz5x3x8dwtNIq/08XKNQS2vP6K8IyVaVp4Rnddfc3VofW/1Br/7jbl9u17Y29oTRvX+fWdjlgdmvsCe0GUgRssROt49qncLfSpWgzdP//kecEJxwfXB7Z39JT24AR/p7e3uX1TbEBr5zuQTvvPT8bmllTddmSFMnCztlgddz2WuwWolFoM3enAv8vdATU38P6D40MOiI+AzuwzfH2W3H95+K/CM8r3nDufttlrwjNawxMos9VP3BcfcmDnl39+hSGwSNsm/yw44bjg+lCsWgzd0Tvg2v2W6A1DUrJpyPbO/9vBuVuA2rlz45bwjNZDwyMosdvW3x8fMjY+AqCv1WLotgMu0SafdE7uFhrFjTTqLHpH1q52Y12+Jvg1gpRSmmTpfJOtHBh/TjflEH6Dxc2Vxpt66cm5W+hTtRi6r/t67DEok8eFlgf+CzfSipVrCKRYC4fsk7sFCBN9szYlN2yBnqvF0A3Uy18vPiq0/gPjQsv3iE2AAHrOzVqgzAzdQOlcdP9nghM+EFy/+37239cWkGLX9jI47IIJuVsAAApUi6H7vct+n7sFoA9dO/mvQ+tPDq3eM4//8vDwjEFHdH4911mYc454XZbc3ObPiz2n+6xDZnZ6Pfp88JScEU4+0WeTp9T1+eRNtPLxR3O3AJVSi6EbqJfH28+IDTghtnxPDJpdxJPu9xSQseOcClAs54NTZ9OmfbOIlAIyquGxFXeF1rcDBXVj6IYdsPj6+HfF3jf+reEZVfHOh38QnHB0cP3ue/72+KG7/7jOr485NTy6U2NHn54nOLOOvxyQuwUAoECGbtgB/7flmgJSjisgA/5ce/vW0PrjW0PLV86YYfNztwC103LeheEZK74QHgH8f2tXx3+vHNTF61gRDN2wAz75tn/M3UKjNPEJ6CMnLwzPGPLwZZ1eH37qYeHZnem467EsuU0VfTRcSo6HI5+m7hGRywcOujG0/q9Cq1MFW0e9O3cLfaoWQ/e6a74bWn/f0OrAf2UYK9b134/9eDNhZOfX33nG46G5QHPYI6JYC4d465pY0Z9NUur680mEWgzdEG3cmi3xIYfGR0BnBmy5KncLFCDXLvVA/USfhuAkBOr22aQWQ3f03U13NoFo7/n4P4dnPBCe0D2rn7gvPuTA+Aggvya+lrTk+jXZsp2GAN1Ti6G7if/QUqxdti4vIMU3sBdtm/yz4ITybVo36PW1+Oe4W25bf398yNj4CCib3W67NDZgWvmOzhp+fOyGkECxDp61d+4W+lTzPuVBxSy59d7cLVCAhd8+LT7kCw91ennTI9FPS0YG12dHrD88dwdAX2p/5g/hGUeGJ0Az1GLoXjsiekfW0cH1Kbucu/6++5jgIxOeOz+2PqU3cY5XaJqgbu/H0bX3Lvt9ltwzUvzxO1eGJ1TH0MVT4kNO6ej0cvTnIichUDe1GLrDlyhanth4t4yM34DoxC6u/+VXOn86Sb1En8KQUvlOYlh5UPmW+QM9t3zNpviQSZ2/irXwvlvCo1v2PiI8ozsWbzwzPMO/0tA3ajF0Q7Qxw4KfNvMKbqQ1xP2/iM94VwFPgqBkzjv26ND63w6t3jNN3CPinN9eUECKfZOgL9Ri6Pa0hDrr9+7dcrdAARp5xuz6u3J30Ch125QGmu6W8z4fnlG+LfOgmmoxdEOdPXLywtD6Qx6+LLR+T6wc+ObcLRTOKQxAX/HvSTN8/NF8ZyJ2PBv/2h3UiaEbSm7Iw80bQB+7dl5o/X0m2YqnDB5bEf+ke5/whO5rX/V07haouSt2+Wpo/db0g9D67Jgb7pkRnjE5db6RGkQ7Z+Pj4RmtwwaHZ7yoHkN39HuB3gmEQl04OfbvdBlH7o67HsvdAlATTXztrolf8+TrBuZugYJs23Zz7hbopXoM3d4LJFj/9/xF7hYa5azNuTso3onn7pG7hcI9MOfB8IwyPumGcB5GQM08krsBeqkWQ/d3D4/dWbo1tDpVULclLlAGv37qmfCMMWnP8IzusqEZ0Zr4is5jF8afEb7P991sAHqmFkP3nRu3hNZvPTS0PBUwf9uHCkj5badXH943dtnxvg+HlqcCom8qdXVDKfrf7pT8+00zbT22gOP4yO64G+P3fLnphPAIdsBPbog9r+7ISaHlSTUZupev2RQbMGl4bH1K7223fic+pIvTnJp43urm9uC/063+TgPUSdvsNeEZZVv5+PWZd8eHnBAfAU1Qi6F74RBv7RHrh5v/GJ7R1cK4ySedE55Nfvduvi88o2XvI8IzeG25VhZAnV3R9pvwjNZD9w/P6I7om/IpdX1jfuqlJ4dn87JHR5wfnFC+EwmW7bekgJSRBWS8oBZD99XXXB1a/0utsUdvAAD0lSauFmriA5jT5njHvCk2rfhsbMDbYsv3xIYH54ZnDB16ZnjGi2oxdM/6cPn+8adeZo5aV0BKF+vLG2jhmltC64/yew1Axf3drfEb2v3qY1ZIlcHQxcE3WE5xHnu0WgzdAFW3+on45eXpwPiI7mji0XA5OYWhOTo2xn6AduOyHHKeTb52dezJQYMOid+NvkqOW/Dp3C3QS4Zu2AHfOPr34RnTu7i++Pp7QnPfN/6tofWhK+FLYFMq5TLY+NMQOj8JAeos+lXDlMr3umHOY9La27eG5o4v2651mc18b7/Q+nPvCC1PqsnQvfzm60Prnzhjamh9eDXvPCP+6RTkEP0aQUrlfCI3cNfYJ0SUw7Rp03K3AEBJ1GLohmhXLp0XnjH9gNnhGQBAPTww58HwjOZtTwcxDN0AJXDb+vvjQ8bGRwDNlXMH8SZuqvvrp54JzxiT9uz0esdfDgjP5mVfm7Itdwv0Ui2G7o5nn8rdAtCHmvjhKeeGOLk08b9zSikNufGPsQEtseXZMR895oxs2Ws3r8+S29R9GnK5c+OW8IzWQzu/PmaY12SK9INDJoTWnxxanZRqMnRDnV32y8Wh9T9+4Gmh9dlBt30zPuNdnW+Is2y/JcHBI4PrQ/ks3fcH4RlHpneEZ1BeOU+AaF/1dL5wqCBDN5Tcj2+KHbrTx2LLU34bHpwbWn/o0DND60MZvW2P2JMnyij6qLKUut4cMXpT3ZRsrAv0nKEboOFy7abtQzJ1tmnFZ+ND3tb55VyvbuRa1g619+WjYutfF3/DrOlqMXRPvfTk3C1QgPN+cV7uFrLYeuwvcrdQOMcAUmcbht0dWn9IOjy0PjvmU7/7egEp8WdTV4X9fYq16ZE1wQkjg+tXy+TrBuZugV6qxdBNM5y/541FpBSQ0T3hG9PYlIZMmvohud8bj8zdAjXnxiXRJs7pfJk/MVrOuzC0/oovhJYnGbqpkGu3DQvPsHtjOTRxGHtsxV3hGc5bBSI19UQCiDbniNflboFeMnQDQAPE71KfkiWhAH1vwYIFofXHj7dyIZqhm8o48qZPxYfMjo8AgEi5Vgs99W/t4RmDZoRHsAPO2fh4aP3WYYND61fN2NGn526BXjJ0UxkL9h4QnuFgI5qora0ttH5LS0tofSijH69cHp5xXHgCUAbDj9+auwV6ydBNZQxdPCU+5BRHJpTBmFNzdwBQXU51IVr86yojg+tXy9/demVo/V997IjQ+hi6gRJqb4+9ozu+NbR8j7TNjj5+JaWyfdkGAwB4bT++aXFswMdiy2PoBiiFK9p+E57Reuj+4RndsXb1/PCMQYfMDM+AsvnVzrGrtiaEVu+ZJq6QCj9SNCXHipZE9I35st2UryNDN5Vx3IJP524Bwiwc4kAvqJPFG+N3Cenqne6O1bHD2ISRoeV7JHqFVErlWyW1cM0t4Rmjkl2tyyD6xnzZbsrXkaGbyvi3mw8Pzzj0g+ER0ClPLKBerrjgntwtNMuhb8ndAYRxY776ajF0z593dWj9syxPJKPoO9llvIs9/gx/55qgiU+maI7n7xsRntF/XHhEt3Q8uj48w63Dl836sN8NqIpaDN00w8G7faOAlMsKyOieJn5T3fT72LNeR4VW7xnLBKFefnBI/JvPk8MTKDPnokN1GLqpjCY+NQCA7uroKOB1lc5yn30qS25TWVkA1WHopjL+tOnA8Iwy/oVYfvP1ofVPnDE1tD4AzZDzGMBcNxpycpMDqqOMM0a3DT/1sNwtAAA02qY9B4RnDApPAOh7tRi6m3g0RhMNOPczBaSsLCCje9zJpq7cMIV6id6PI6Vy7smRSxPPJoeqqsXQDXW2Pv6kNACgYvY63OthUBWGbirjuBvfHJ5x0wnhEQAAQIMYuqHkxgybn7sFCtDEo+EA6DnL+aE6DN0VNW3atNwtFO7x9jPiQ06IjwDye+72OblbgDADtlyVuwVqbsODc0PrDx16Zmj9qvn5zUtC67+l1Y3/aIZuKmPOEa/L3QJQE0v//bjwjCPDEyi1Lx8Vn3FdR3wG0FhLrl+Tu4XaMHRX1N2vn5C7hcItWLAgPGP8+PHhGdCZ6PPYU3Im+/aW7vuD8Iwj0zvCMyivayf/dXjG5PAEoMmOfMfJBaSs6vTqwF3r9XplLYbujo7YI8PK6OdXPJa7hcKNHX167hYgjKPhivU/lr4/PuRv4yMoL69EAVV37bZh4RlNuXlYi6G7iZ5f+rvwjP5HhEd0yxW7fDU8ozXFP/0C8luw94DwDG8kNtvCb58WH/KFh+IzAOi1WgzdNgyBvtcy68bcLUCYoYunxIecUq73baM3PkrJ5kfbW3fNd8Mz9g1PoMw6Vsev9JwwMjyCHbB28/rcLdBLtRi6aYaVA+PP6eZlP/vvawtIOb6AjGoYc2ruDgAAXrZpxWfjQ94WH1EGtRi6Nz0SvbPeyOD6UD7997mugJRzCsiAPzdx2tG5Wyhc3Tal2VHT7mjm1w3QW5/63dcLSLm6gIz8ajF0T5zTvB2nHXdTrNXrH8zdQuEe/eGV4RmDZodHQKeev29EeEb/ceERlJgTN4C+MuvDec7R/vHK5eEZXU00bW1t4dktLS3hGS+qxdANAFAm/21A/M1xXmZ/H6iXs0d/KHcLfcrQDTvghmPuDc+YFJ5QHdGbw9gYBoj2f1uuKSClXIO9jb0gxlP/1h5af9CMzq8ft+DTobmv5smlm8MzhhwbHvGSWgzd52x8PLR+67DBofUpv0/fuDI8Y1LLB8IzgBf84JAJ4RllO3u0bkv1YHvx+/ukZI8fmubfbj48POPQD4ZHlEIthu4mWrpv/HnSR6Z3hGdQXoNmF7F7+XsKyABSqt9SvbJbeVC5nkLXXRP396FYjlKlNwzdFXXI44fkbgEAoPGeu31O7hYowD8cvVPuFgp38G7fKCDlsgIy8jN0A52a+d5+4Rlz7+j8ekdH/HuB0DR1ez8OyuK6Qb8Krd/Vqyq+VxLNSR99pxZD97L9lgQnjAyu330TD/xcESkFZFBWV1xwT+4WACrrsQtnhmfs8/0p4Rm8tk0rPhsb8LbY8kC8WgzdTXTttmHhGWXbBIhiPXLywvCMIQ9fFp7Baztt6/tytwC188CcB8Mz9glPAJqs5fVXhGesTNPCM8rA0F1R4XdVU3JnFQDgNXzqd18PTri606vOJifasg3xQ3dXNgy7OzxjSIrfnf1Fhm4oufZVT+duAUJc//3Y9yBTcqYvEG+nM0dnyXVMGtFuuKeLA7z70OTUEZ5RBoZuKLmDZ+2dJbf9Z2eEZwzp4rq790DV/fqpZ8IzxqQ9wzN4bdFnGZfxHOOBu87P3QIFmHzdwNwt1Iahu5eW3Hpvltz4pUwpdbWcKZfHrp0XnrHPpCvDMwCovyvafhOe0Xro/uEZVXHOxsfDM1qHDe70+qhHfxqcPDy4PnTuuBvfHJ5x0wnhEaVg6K6oXEuZAOid9kWPhtY/7Pj4I154bQuH2OYMgBcYuisqeilTSuVbztQ2O/7dpdbwBHZEriX1AADQ1wzdvfQ/b1obWv/nR8UfDQbk196+NTxjvLtKAH3uX7bF1u/qQKWJc8bHBgN9xtANQPN8+aj4jOuasSMrAPX0eHv8prrphPiIMqjF0L3hwbmh9YcOPTO0PgDFunbyX4dnTA5PAIA473z4BwWkHF1ARn61GLop1k9Ojt/JkJdF78ja1W6sv3+qgJHhgPgI4AUbht0dnjEkxe83UhU/v3lJeMZbWu1qDcQZO/r03C3UhqEboASGn3pYtuzoYcwgBtD3ch6TRjOsHXF9ASnNOJHJ0A106pnzb48P+eHH4jMqomP1pvCMCSPDIyixFbcuCs8YO3VmeAbl1dER/+9YGe1226WxAdO62koNYt22/v74kLHxEWVg6K6og3f7RgEplxWQQVn968g3hWdMD0+AzuXcHGbmqHXBwXY0pnmW7Re/nD+lkQVkVENbW1to/ZaWltD67JiVBx2Xu4XaMHRX1PP3jQjP6D8uPAI6les99pya+oQIqI8BW67K3UIWk68bmLuFwp09+kO5W6AI9/8iPuNdUzq9XLfVWYbuioo+EzKlrs+FzOXOjVvCM1oPDY8AgFra9MiaAlJGFpDRPcfdGLvB7E0nhJYHCmDorqh/efivwjPKNnQD9JU5R7wudwsAQEMYuitq4bdPiw/5wkPxGbym+PfURgbXh86NGTY/W/a6ddHvVac0frx3q8kj/onzyOD67AjvsUN1GLorat013w3P2Dc8AWiy9lVPZ8sesPO7s2UDAM1i6AYogZxPfZto+PFbc7cA1ET4aQgnxJYH4hm6K2rBggXhGWVbGrl8TQG7O08aHp9REe98+AcFpBxdQAZldfCsvbNl37v5vvCMlr2PCM8A8ov/ful7JVSdobui/tsA5+YRy5+xYuVcat1Eq5+IH7rTgfERlNfazetzt0BBPvvZL+duASg5QzfQqRPP3SN3C43S1KN2mqjj0stD6w+fPiO0Pjtm1oetnGqKZbf/NrT+sGFvDK0PxKvF0D1wV+9CNsHCIfvkboGCNHHH9olzyvU6R93dtv7++JCx8RFAfmtHXB+cMDq4PhCtFkM3zbC5vYB3uls9mYAmWHmQ1yeItfzm6EEspRNnTA3PAJrrsWvnhWfsM+nK8IwyMHT30tFT2oMTjg+uD9BA9/8iPuNdU+IzumHMSZfkbgFq6eK/fGfuFoCSM3T30vl73hidEFy/+24/43+E1h/XxfV9htuUpkgL77slPKOr3Z03PDg3NHfo0DND61N+Oe/e39F2amhuV+t12hc9GpqbUkqHHT8iPKMqOp59KncLFOSxC2eG1t/n++W6gQd0n6GbbnvbHvfkbgEAAAjUNjt+k9fW8IRyMHRTGU19Py7XU18bTQH03NRLT87dAgAlYeiGksu1O//Kx+OXotJs52x8PDyjddjgTq8/MOfB8GznLQAAKdVk6G5rawut39LSElq/J3KeoXzkTZ+KDZjd+WXvxxXrsRV3hWcYSl6WcwBtol8/9Ux4xpi0Z6fXr1wa+z759AO6+EcUCBF9E8/3ynJ4dth3C0iZWEAGOdRi6KZYb9rpX4MT3hNcHyCf1fv6UFWkb9wWvMz78JNi61N60TfxurqBR7E2rfhsfMjb4iPIw9ANQONc0fab8IzWQ/cPzwAAys/QTWXYlAYAKJvom3hu4NFEY066JHcLfcrQDQANMHPUugJSxheQAQDVYuim2/q9e7fcLQAAQGGGLp4SH3JKR3wGWRi66bZHTl4YWn/Iw5eF1gcA6CtXP5End8ODc8Mzhg49MzyD8rpz45bwjNZDO7/evij+6NrDjh8RnvEiQzeVMX/e1eEZZx0yMzwDAKiPhWtuCa0/ymsbUHmGbrpt0Pmfz90C1M6y/ZYUkDKygIwdN3/bhwpI+W0BGQAAXTN0AzTcilsXhdYfO9UKEgCq7bgFn87dAhVm6Kbb/rTpwND6/lCWwwNzHgzP2Cc8gTIbuOv83C0AVJZ/Q4v1/O1rwzP6jwuP6JblazbFh0waHp9RArWYb84eXcQSxc49+sMrQ+sPmh1avkfGPfefofVXhFZnR/36qWfCM8akPcMzeG0D9rJRDgC8mufvi990q2xDN32nFkP3k0s3h9YfcmzXPzZodvRdr/cE1+++OUe8Lkvu+DMsUQXoqY5LLw/PGD59RngG5TVxTjM3/Jr14WY8qWu6HxwyITxjcngCudRi6M5p475fDK0/5OFpofUBcmlrawvPaGlpCc+AzqxdHb/0d5ATN0ph+c3Xh9Y/ccbU0PrsmPdOHBMf8vBD8RndsHCIFwH7iqG7l06Ys2do/ZWh1Xtm3bp1ofXHj2/mnfKyuaLtN+EZrYfu3+l176kBVbfX4Qalpuh49qncLVCAmz73P8MzyrZWaHN7Ae90tzZjpYihGwCgj236fXt4xqjwBHbEmFPz5FotRJ3NHBX7kO8FxT3oM3RX1L2b7wvPaNn7iE6v33nnneHZAFBl13//V+EZE0aGR7AD2tu3htYf3xpanh3UxM+/+wxfn7uF2jB099J5R3w8dwuF+28DjsvdAgXwHg8AwAt8/qU3DN1Ap7zHAwCvzekqzXDiuXvkbqFw0ZsEptScjQJrMXRvGHZ3aP0h6fDQ+j2x+on45eXpwPgIyqtjY0d4xqgC36WBprujLf7FU7fRXjZgy1W5WwBqouPR2GXe/u2OV4uhG+h7azd7jwegpzY9sqaAlJEFZOy4czY+Hp7ROmxweAZ0Jud+Srl2yLczf98xdPfSI58Kfu/1h51fvm39/bG5KaU0tvPLuZbXLPn5LeEZo4Y2Y4nLjpj1Yfc9myJ+h1ArGmieiXP8uQfgBYbuXvrXkW8KrT+9i+srB745NPfVRN/p6+ouHwAANM2AUbE75Hdl6qUnZ8mtI0N3L8054nV5gtfflSeXxrB5BgBAfjbrqz5Dd0U9tiJ+6C7bgVEdHQXspg0AAP/FZb9cHJ7x8QNP6/T6pt+3h+aOCq3eMx2XXh6eMXz6jPCMFxm6K6ptdvwGLa3hCZSZzTOKteHBueEZQ4eeGZ4BAHW08iDndNNzhm4qw/Er1NnAXefnbqFRrn4idwcAUG5rV8d/Nhl0SDOWzhu6Acji7NEfyt0CUBNWCxHu/l/EZ7xrSnxGN7S3x2/gNr4hS2sN3RV158Yt4Rmth3Z+ffUTwecUHtj55SaeeZpSSm1tbaH1W1paOr1ux0qiPbl0c3jGkGPDI4ASyLlaKPpIU8eZQvUZuum2i//ynblbAKiuZcGvyrz/otj67JBzNj4entE6bHB4Rncs229JASkjC8jgtUTfNHXDlLoxdEPJWYILfa9jY0d4xqg0PjwDAKI4qqzvGLor6qzNuTugKO4mEy36WI4ij+QAACgbQzcAAJUWvf9JSl3vgQLwWgzdFbW5fVN8SOvw+AwgpeQDI0BveBULKDNDdy8tWLAgtP748Z2/E7hwTexOmSl5HxEA2HGO7QLonKG7lyafdE7uFqB2ch2TBgAAfc3Q3UuLr78ntP77xr+10+uzPmzpN/VlmSAAwAs6Vse+VjphZGh5kqG7spbffH14xokzpnb+A7d9Mzb4XVNi61N6dmwn2trN63O3AAA0hKGbbvvu4fND67eGVgeAeMv2W1JAysgCMnbcwF1jPx9AU3V0FLCBMqEM3RXV8exTuVsAauKOtlND65fxZRiv6AAARTF099KJ5+6RuwVqbsOwu0PrD0mHh9aHrkT/2U6p6z/fWV/RAQAaxdDdS/duvi+0fsveR3R6fUzsgylKZMWti0Lrj506M7Q+ANSZpb/U1ZKfF3BE8dBm3KA2dAOdWrs6/t28QYcY+MnDKzoAdMv6u3J3QIUZugHIInoVR0pWcgBAT1nF0XcM3RW11+HNWIoBEGHAqK25WwCgQh5bEf+ke5/wBHKpxdDtnVcAumP8Gf5dB+gpm7wWa8CWq3K3QC/VYugG+l57e/yTwPEOZQeojba2tvCMlpaW8AzKrWXWjblbgG4zdNNtd27cElq/9dDQ8uyg4acelrsFAHpg4pzxuVtoFE8hi3XGrO8FJxwfXL86/NnuO7UYusecdEnuFgAAGu3s0R/K3UIWmx5ZE5wwMrg+dC7+z3ZKTfnzXYuhm2ItXxO8k+Gk4bH1oYSa+mEVAKDuDN29dNkvF4fW//iBp4XWh65c//1fhWdMGBkeURlPLt0cnjHk2PAI6NS0adNyt1C4DQ/ODc8YOvTM8AwAes/Q3UsrH380dwsQYsyw+dmy7YoK9XL36yfkbqFwA3fN928oAOVSi6G7fVHs4HvY8SNC6/fEpt+3h2eMCk8AAPraORsfD63fOmxwaH12jNeSoDr6fOj+05/+lHbeuRazPAA15cYlAN3xwJwHwzP2CU8gl25Px6effnr62te+lnbdddc/+7FVq1alT3ziE2nVqlV90RsldfUTeXIdgdIcK25dFFp/7NSZofUBKJZ9MYAy6/bQPX/+/HTHHXek733ve+nQQ184UPm5555LX/rSl9KFF16Y9t133z5vknLZZ/j63C0AQKmNfGpg7hagljat+GxswNtiy9NM3R66f/WrX6UpU6aksWPHps9//vNp4sSJaerUqWn16tXpk5/8ZPrSl74U0ScANTNgLzsvN8XP/vva4ITjg+tD1zZtiv7zDVRdt4fut7/97WnFihXp4osvTuedd146//zz09vf/vb0y1/+Mo0a5Q22Jlh+8/Wh9U+cMTW0PgDFumCX/wytf35odQDonR7tePbb3/42LVmyJD333HNp8ODBaf369amtrS39zd/8Tdppp536ukc60bF6U3hGV2codzz7VHh2Z6J3Y03JjqwAVN+y/ZYEJ4wMrs+O8B47VEe3h+4LLrggXXTRRekNb3hD+tGPfpSOOuqodN5556Vzzz03XXfddelb3/rWS+96U09TLz05dwsAVMjb9rgndwsAvfLrp54JzxiT9gzPII9uD93nn39++tjHPpYuv/zytOeee6aUUpozZ06aNGlSOuWUU9Jhhx2W/vjHP/Z1nwAAQMN96ndfD064Org+TdTtofuHP/xh+uAHP/hn1w899NC0atWq9PnPf75PGgMAoPs2PDg3tP7QoTZBJJ9fDftMaH3rdYnQ7aG7s4H7Ra9//evTxRdf3KuGqmbDXt8Mrf+WLq53dMS/092Vtavnh9YfdEjnZyjHv6OWkvfUoDgzR60rIGV8ARnVcOXSeeEZ0w+YHZ4BkMOdG7eEZ7Sa+GurRxup8bKrr4ldgvKl1q+G1u+JraPenbsFACrEubpE2zDs7vCMIenw8AygngzdvTTrw8Oz5A7YclWWXJqjfdXTuVuAWrrhmHtD608KrQ5ds5s2dXbW5twdUGWG7l7KdWb1pkfWhOa+YGQBGUBKntJQbzY+AqDJDN1Apw6etXfuFoCa2Lbb6NwtFM5TX6K5WQvVYeiuqE2b1uZuoXDRu7GmZEfWshiwl/8OAPBqVty6KDxj7NTON7dtos3tBWxi3JrntVXiGbp7qePZp3K3AAAA1FT8a6Ujg+tj6KYyBu4ae1QZr3TOxsfDM1qHDQ7PAGga3y8BysXQ3UtjTs3dAcSYv+1DBaT8toAMXkv0+c3Obgboe2NOuiR3C8AOMnT30oL7YuuPjy1fKW1tbeEZLS0t4RlV4UkJ0TouvTw8Y/j0GeEZ3fXpG1eG1p/U8oHQ+j3xp00Hhtb3YQaItnDNLeEZo7r45D9xTp6JIFduHfk+1UtjhhlMgGpbve/E3C0U7vrv/yo8Y8LI8IjKuOGe2Jsfk1NHaH0A6A1DNwAQ6mtvvyy0/uTQ6gDQO4ZugO385OQ3524BAKDR7miL3ziryAPaDN1Ap0Y+NTB3CwCVtWHY3eEZQ9Lh4RnAC2Z9ON8Z2tEnynR1moyTbPqOoRtKrom7k84ctS44wcYgUKSF3z4tNuALD8XWr5jojUdtOkoTLb/5+vCME2dMDc8gD0M33Ra9AZHNh6AZ6rZ0DJrO033qrOPZp3K3QIUZuoFOPbl0c3jGkGPDI6BTA7ZclbsFACpk6qUn526BCjN0022OSQMAgHpbtt+SAlJGdnr1yqXzwpOnHzA7PONFhm66rX3V07lboABNXSb43tM/lbuFxsj5DfXgWXuHZ/Oyddd8N7T+vqHVAVJauzr+odOgQ2Z2ej1++B0ZXB9DN93mwyp1dt2g2D0LnCcMQNWtuHVRaP2xUzsfPptqw4NzQ+sPHXpmaH0M3UAXor+hplTOb6of+fHusQFviy0P0ERN/Z4FVIOhG2A7P165PLT+caHV2VHOHgWAVxf9hD2l5jxlr8XQ7UxfACivBQsWhNYfP973aSjSgL3yDUqjHv1pcEL5DpscuKtNjKuuFkM3xYp+QuTpUDmMOemSbNkdl14eWn/49Bmh9XvC0kjqbOzo03O3QM3l/J4FdWXY7zuGbqBT7YseDc847PgR4RmU1+p9J2bLnr/tQwWk/LaADIBmaXn9FaH1V6ZpofVpJkN3RZ3+r7Nyt1C4s0cX8SEZAKrN90vq7Kxb/hQbcEpseZrJ0E1lPLl0c3jGkGPDI4ASsGQOoJqunfzXofUd7fmytra28IyWlpbwjDIwdFdUzp13l+23JDh5ZHB9AACq6PH2M2IDTogtTzMZuum26OMDmnJ0AABEiF4Z1tWqMBtCUoQ5R7wudwvQbYZuoFPxR/Gl5Dg+INq9J8X+W3ZgaHXgv9pvv7/J3ULhopd5N2WJd06G7oqKX+KdkmXeNNHijbErLY4LrQ6Uhae+UC/t7VvDM8a3hkeQiaGbytgw7O7wjCHp8PAMgKbpuOux3C0AQDaGbgCAPjbmpEtytwBASRi6AQAAXsX4M7zOQc8ZugEAAErq7NEfyt0CvWToptsG7jo/dwsUoOPSy8Mzhk+fEZ4BkEP7okfDMw47fkSn16P3QLH/CTn9bo97QusPS0eE1qeZDN0ANE708SspOYKFfKJ3TrdrOvHHijpSdHtPLt0cWn/IsaHlSYZuABrIUj0AoCiGbgAAqJicrzAA3WPoBoAGWL3vxNwtNEr88tuULMGliRbed0to/Za9vdP9IqvC+o6hG+jUHW2nhmcMD08AoAk89W2Oi//ynblbgG4zdNNt0RsQ2XwIAIBO3fbN2PrvmhJbv0KiN3BLqTmbuBm6gU5duXReeMb0A2aHZ3TXOb+9IDjh9OD6UD4nnrtH7haAmnhsxV2h9fcJrU5TGborasODc8Mzhg49s9Prp219X3g2L4teMme53CudMGfP0PorQ6tXzLKr4jPef1F8BgAE2jDs7tD6Q9LhofUxdANdsOlSsQbs1flNLgAAqs3QDbCdwXvsn7sFgF7JdRNvzEmXZMmlWdpmrwmt3xpavVqin7Cn1Jyn7Ibuihq46/zcLUAt/f3uZ+VuAWrn3s33hdYv4xE/HZdeHp4xfPqM8AzKy7F0UB2GboDtrFsX+yFm/HgfYAD6mgG0Oe7cuCW0fuuhoeVpKEM33WYzB6JFnxHufHAAqKapP9kUGzCpfJ8SVty6KLT+2KkzQ+tj6K6s6LOyU3JedlnE373v4s69naWbI/q/tf/OAECDGbor6uzRH8qW7W4b9D1LIwEA6snQDQAAFdPUzfo6NnaE1h/lBjUBDN0V9eTSzeEZQ44NjwCgAVY/Ebt7eTowtnxPRO9NkZL9KQCqon/uBgAAAKCuDN0AAAAQxPJyKLnod7bK+L4WAD2X7dQLADpl6AYogaZuiAPQF/wbSrRNv28PzxgVnkAuhu6K2jDs7vCMIenw8AwA6u+29ffHBoyNLQ8AvWHopjKizwdPyRnhpPT+9384dwsUwAkQxVp50HG5WwCAbAzdANv51teWh9a/5ErDBwBAkxi6K8pTX+rsyqXzQutPP2B2lz924rl7hGbzshuOuTc8Y1J4AgDAqzN0UxljTrokdwsAAFAon4Grz9ANANDHolfspPTqq3aa5o62U8MzhocnUGab2s+KDxkZH0Eehm6A7Sy875bQ+i17HxFaH0rptm/G1n/XlNj68CpO/9dZuVugAO2rns7dAhVm6AagcRy7CEB3HDxr79wtUGGGbgCAGum49PLQ+sOnzwitXzXnbHw8tH7rsMGdXrekvjnaFz0aWv+w40eE1sfQDUADOQECgO6IvrmSUtc3WKg+QzfAdm5bf39swNjY8gAUa/62DwUn/Da4PhCtFkO3ZVRAX1l50HG5W4DaeWzFXaH19wmtDpDSsv2WFJAysoAMcqjF0A0AAEDf8SpW3zF0AwCh2mavCa3fGlodIKUND84Nzxg69MzwDPIwdAOduuGYe8MzJoUnANAEVy6dF54x/YDZ4RmU18Bd5+dugQozdFfUmJMuaWR2E0UfCeI4EACAV9fW1hae0dLSEp7RHT7z9x1Dd0VFn9eXUtdn9jkrkFq77Zux9d81JbY+AFArM0etC04YH1wfQzfAdn5x3M2h9ceEVodyGrdmS2zAobHlKb/V6x/M3QJAlwzdvdS+6uncLVBz0e+pdfWO2qdvXBmam1JKk1o+EJ4BnRmwl81qirS5fVNsQKsXZQAoL0N3Lx08a+/cLVBzq/edmLsFoA6WXRWf8f6L4jMAoGIM3XSb90oAoLxybcCZc78ZitVx6eWh9YdPnxFaH4pm6KYyfDOnCL9+6pnQ+mPSnqH1ASiWo6Sa4ezRH8qW7SZH9Rm6e+mcjY+H1m8dNji0PgBEW7jmltD6o6yQIqPoo6TKdoxUUz25dHN4xpBjwyPIxNANlE7O99jv3Bi7y3KrXZYBABrF0N1L87dFLzX5badX49+rTqmrd6stcYF6ueGYe0PrTwqtDgBQbobuXvIeDwC8ulkfdqQXUG0bht0dnjEkHR6eQR6GboDtnLU5dwdQP1cvuDq0/j/MmNrlj+2+4ZzQ7HTgwtj6QCmsuHVReMbYqTPDM8jD0N1LNs+Aetll6/LgBE/8yiDnKzpNNPXSk3O3QAH8vSrWlUvnhWdMP2B2eAY0gaG7l07b+r7cLQB9aPnN14fWP7GLJ3LR5+qmZNwnn7WrY1/FGnSIp0NArDEnXZK7BSrM0N1LD455Q2j9IaHVAYC6iX4C6unnK+U8v5nitC96NDzjsONHhGeQh6G7l6L/Anb1ly96B/GU7CJOMx12wYTcLQBQIdHnNzu7mVzcaOg7hm6A7VgGCwBAXzJ0023R755675S07KrY+u+/KLY+AI1x3fqLQ+uf5RgpqDxDNwCN08RXdG445t7wjEldXG9v3xqaO741tDwA9Iqhu5fij8co39EYNmgBAADYMYZuAAB6rYkrSAB2RP/cDQAAAEBdedINAAA9NOakS3K3QAHiXylNqYyvldI3DN0AAAAl5eSg6jN0V1T0X76UyvcX0B1GAIBmW7t5fe4WoNsM3b0UvWmIDUMAAACqy9ANAABUwskHbM2Sa3d+esPQTbet3ndi7hYAABqtqZ/HOh6NXV5ettcrqYdaDN02FwAAgPrrePap3C1At9Vi6M4p18B/5dJ5obkppTT9gNnhGQA5NHEzSiBG+6JHQ+sfdvyI0PpVM/XSk3O3AN3WP3cDAAAAUFeedAMAALwKK6ToDUN3L03/zpdiA863xBsAAHJq4qudM0etKyBlfAEZ+VleDgAAAEE86QYAACph7er5ofUHHTIztD7NZOgGoHGauEwQAMjD0A2wna2j3p27BQAAasTQDQBAr9ndGaBzhm4qo+PSy8Mzhk+fEZ4B0DTjz/COJADNZegG2E7H6k2h9SeMDC0PQMHij1VqxpFKUGeGboDtdHTEDt2Uw+p9J+ZuAQBoCEM3leFdMYowYMtVuVsAAGi0ut0cN3QDbOfgWXvnbgEAgBoxdAMAQA9Fb/Rqk1eoPkN3RdVtyQWUxTkbHw+t3zpscGh9AKizvQ6fmruFxnByUN8xdAMA0GtXLp0XnjH9gNnhGQB9zdANsJ0bd9mauwUAKiR6o1ebvEL1GboBSsATIgCgM9GfEbr6fODkoL5j6AbYzpAb/xgb0BJbHspoyc9vCa0/aqh3PMkn10AEVIehG4AsPn3jyvCMSS0fCM+AsmniZqs3HHNveMak8ASgrgzdANvZMOzu0PpD0uGdXm/ih2QAqsl77NA9hm6gdKKfWLza04oVty4KzR47dWZofQAAysXQDbCdMSddkrsFAABqxNANsJ32RY+G1j/s+BGh9Xsi+t1q71UDAE1Wi6HbrpEAAFAMn72he2oxdAMAAMWw+Sd0j6EbAKBOll0VW//9F8XWB6gZQzcAAL3m6SdA5wzdvZTzaCMAqIKOjk25W6AI0U/YU/KUHaik/rkbAAAAgLrypBsAAHrIsnrgtRi6e8n5tlAvM0etC04Y3/llyzIBqAqb9UG3GLqrygd0AACA0jN0AwChBmwp4EYxAJSUoRsACHXwrL1ztwAA2Ri6AYBQ52x8PLR+67DBofUBoDcM3QBAqGX7LQlOGBlcHwB6rhZDt6MaAAAAaqJmm0bXYuimGa5cOi88Y/oBs8MzAACA5jB0AwChnnzysND6Q4eGlgfIyqre6jN0AwCh2traQuu3tLSE1geA3jB0A9A8NXtXDAAoL0M3AAC9dsMx94ZnTApPAF5kP6W+Y+gGAKiR6OHX4AvQPYZuAACgEjb9vj20/qjQ6jRVPYbu6HfzvJcHAD122tb35W4BALKpx9BNsTLd5GjscQluKgEVt2HY3aH1h6TDQ+sDQG8YuivKZiUAVMWKWxeF1h87dWZofQDoDUM3AEAfa+zqLAjWsXpTaP0JI0PL01CGboASsHqFOhtz0iW5Wyies+AhREdH7NDdFTfS6A1DNwA0wKdvXBmeManlA+EZABTDjYa+Y+gGAEK1L3o0tP5hx48IrQ8AvdE/dwMAAABQV550A9A43qEHAIpi6KY6bEoDUEkzR60LThgfXB8Aes7ycgAAAAjiSTfdFr0s05JMondZfrUdljsuvTw0e/j0GaH1oYz8vQL6ylP3/V3uFprDKtM+U4uh2xAIAAD1N3GO10monloM3VBnbioBVXdH26mh9YeHVgegaHXb8NTQXVHRy29TevUluAAANFPOBwLnbHw8NLt12ODQ+jSToRtgO57IQd+7cum80PrTD5gdWr9qcu6LAcCfM3QDAACVcOMuW3O3AN1m6AYAACir6F3EG7KDeE6GbgAAoBLa2tpC67e0tITWp5kM3QAAQCWcPfpDuVuAbqvF0G3DEACA5nKqS7F89obu6Z+7AQAAAKirWjzpBgAok+hzjFN69bOMASgPT7oBAAAgiCfdADSO9z8BgKIYugGAUKvXP5i7BQDIxvJyAAAACOJJNwAAUAlPLt0cWn/IsaHlaShPugEAACCIJ90AANBTy66Krf/+i2LrA+EM3XRb9K6/dvwFAHgN0cN+SgZ+6COGbirjhmPuDc+YFJ4AAEBPXbf+4tD6Z6XDQ+tXic/efcfQDQCEiv7g1pQPbUAz+Te0+gzdAABAJQzY68zcLUC3GboBgFD2AgGgyQzdUHI+rBbryqXzQutPP2B2aH0AIICN6+gFQzcAAFAJM0etC04YH1yfJuqfuwEAAACoK0M3AAAABLG8HAAAKsYZylAdhm4AAKASOi69PLT+8OkzQuvTTIZuAAB6Lfq0jZTKeeJG9BNnT5uh+rzTDQAAAEEM3QAAABDE0A0AAABBDN0AAAAQxEZqANtZvf7B3C0AAF24o+3U0PrDQ6vTVIZuAAAASqNupyEYugEAAEoqegAt41F8dWPoBtiO81YBoLym/+bfYwMO6Pxy9OeDlHxGqDMbqQEAAEAQT7oBAKCHLP0FXouhGwAAqIZlV8XWf/9FsfVpJEM3lVG3XQxhe/58AwDUk6EbAACAV/BAoO/YSA0AAACCGLoBAAAgiOXlAABAJUSfl+2sbCIYugEAoGK8bwvVYegG2I7zVgEA6Eve6QYAAIAgnnQDAAC8Csv56Q1PugEAACCIJ90AAH3MUzEAXuRJNwAAAATxpBsAAKgEp4xQRZ50AwAAQBBDNwAAAAQxdAMAAEAQQzcAAAAEMXQDAABAkGy7l2/bti2llNIzzzzT61rP/+d/9LrGq3m1HnNlR+fmzG7i15wz29dcjuwm/n7X+WvOme1rfqVnt/wxNHtXv9+F5TY129dcjuwm/n7X+Wt+teye1Hhxtu1Kv22v9TOCbNy4Me233345ogEAAKBPPPTQQ2nIkCFd/ni2ofv5559PjzzySNptt91Sv379crQAAAAAPbJt27b07LPPpsGDB6f+/bt+czvb0A0AAAB1ZyM1AAAACGLoBgAAgCCGbgAghDfYAMDQDQC9tnLlyjRlypS0//77p1122SW99a1vTdOmTUvr16/P3VohHnjggdSvX7/0ne98J6WU0ubNm9PHP/7x9NOf/jRvYwBQAoZuAOiFuXPnpjFjxqTHHnssffnLX06LFi1KZ599drr99tvTIYccku66667cLYbbZ5990s9//vP0gQ98IKWU0qpVq9I111yTnn/++cydAUB+di8HgB762c9+lo444og0a9asdNlll73ix5544ok0atSo9OY3vzndeeedeRrM5Pbbb09HHnlk+slPfpLGjRuXux0AyMqTbgDooa985Stpzz33THPmzPmzH3vTm96U/umf/imdcMIJacuWLem5555LV111VTrooIPSLrvskvbff/909tlnp61bt770az7xiU+kY445Jn3zm99Mb33rW9Muu+ySxo4dm+699970ox/9KB100EFp1113TaNHj06rVq3q9a8bN27cnw3Ft99+e+rXr1+6/fbbU0opfec730k777xz+sUvfpHGjBmTBgwYkP7qr/4qffWrX33p12y/vPzFgTullI488sg0bty4NHfu3NSvX7907733viLre9/7Xtppp53SQw891MP/AgBQfjvnbgAAqmjbtm1p8eLF6bjjjku77rprpz/nIx/5yEv/+7TTTkvz5s1LZ599djr88MNTe3t7+sIXvpDa29tTW1tb6tevX0oppeXLl6dHHnkk/dM//VP6wx/+kGbMmJEmTJiQ+vXrly644II0cODANH369HTyySene+6556X6Pf11O+L5559PH/nIR9KZZ56ZLrzwwnT11Ven2bNnp4MOOigdffTRr/i57373u9PcuXPTzJkz09y5c9O4cePS4MGD02c+85l0zTXXpC9+8Ysv/dzvfve76X3ve1/ab7/9utUPAFSJoRsAeuDJJ59MW7duTW95y1te8+euXbs2XX311emiiy5KZ599dkoppaOOOioNHjw4TZkyJS1atChNmDAhpZTSs88+m77//e+nt7/97SmllJYuXZr++Z//Od12221p/PjxKaWU7rvvvvSZz3wmbd68Oe255569+nU7Ytu2bencc89NU6dOTSmlNHbs2HT99denH/3oR382dO++++5pxIgRKaWURowY8dL/njRpUvre976XLrjggtSvX7+0cePGtGTJkvS9731vh/sAgCqyvBwAemDnnV+4b/3cc8+95s9dunRpSimlyZMnv+L6Rz/60bTTTju9tJQ7pZT22muvlwbnlFJ685vfnFJKafTo0S9de8Mb3pBSemGX8N7+uh01ZsyYl/73X/zFX6Q3velNacuWLTv866dOnZoeeOCBl3Y0nzdvXtptt93SpEmTut0LAFSJoRsAemCvvfZKu+22W9qwYUOXP2fLli3pqaeeSv/+7/+eUkpp0KBBr/jxnXfeOb3xjW98xRC8++67d1pr4MCBr9pPT3/djvqvS+j79+/frd3Jx48fn97ylrekefPmpZReGLo/+tGPpgEDBvRJfwBQVoZuAOiho48+Ov3kJz95xWZo2/vWt76V3vjGN770/x999NFX/Pgf//jH9OSTT77i5xSpX79+f/ak/ve//31Y1ic+8Yl0ww03pDvvvDOtW7cutba2hmQBQJkYugGgh84666z0u9/9Ln3uc5/7sx979NFH01e/+tU0YsSIl5ZQX3vtta/4Of/rf/2v9Nxzz6X3vOc9hfT7X+2+++5p48aNr7i2bNmyXtfdaaedOr3+iU98Im3evDmdddZZafjw4a9Y+g4AdWUjNQDooZaWlvTFL34xfe5zn0sdHR3plFNOSW984xvTmjVr0le+8pX0hz/8If2f//N/0vDhw9Mpp5ySzj333PQf//Ef6b3vfW9atWpVOv/889ORRx6ZjjnmmCz9T5w4Md10003pzDPPTMcdd1z66U9/+tLy7954cZO2H//4x2mvvfZK73rXu1JKKe2///7pqKOOSosXL04XX3xxr3MAoAo86QaAXviHf/iHdMstt6SUUvr7v//7NGHChHTFFVekiRMnplWrVqXhw4enlFK6+uqr03nnnZfmz5+fJkyYkObOnZs+9alPpVtuuSX175/n2/Gpp56aPvvZz6YFCxakCRMmpOXLl6f//b//d6/rvuMd70iTJ09OV155ZTr55JNf8WMf+MAH0k477ZSmTJnS6xwAqIJ+27Zt25a7CQCgGY499tg0YMCAdMMNN+RuBQAKYXk5ABDui1/8Ylq3bl1avHhxn7w3DgBVYegGAMLddNNN6Te/+U36yle+kg477LDc7QBAYSwvBwAAgCA2UgMAAIAghm4AAAAIYugGAACAIIZuAAAACGLoBgAAgCCGbgAAAAhi6AYAAIAghm4AAAAIYugGAACAIP8PIy0zXZ8ZIM8AAAAASUVORK5CYII=",
      "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(101):\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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8QAAAHbCAYAAAANs1OQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdhUlEQVR4nO3dfaxdZZ0v8N85raFQe0IpWO20lQImtqbXizGhpaItDgEKYcAxBCQVpAZsaoK3QGyutSIFfIEwZKRqNFUoRYxBe22QFo3QCqLe2LRlVAI6kQI2vCkNFGFi2jN/YDtUdsvZL2uvl9/n81e79znPetY6a629vs/vWWsPDA8PDwcAAAAkM1h2BwAAAKAMAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkNLosha8Z8+e2LFjR4wbNy4GBgbK6gYAAAANMzw8HC+++GJMmjQpBgcPXAcuLRDv2LEjpkyZUtbiAQAAaLgnnngiJk+efMD3SwvE48aNi4hXOzg0NNRxOxs3vbtXXWpp7ge2tXz9P//z5kKXe+yxn6zUcsvazk38+x5oGx/3s4cKW2ZExB/e/79avp5tnyprO5exvpnWNaK8fTnbsVvG+bGJ63qw5WY7hpq4vlVbV+ep+i/XPtW+F154IaZMmbIvdx5IaYF47zTpoaGhrgLxW466o1ddaulAfRs3bkwpy/3L86sKXe7xQ/+35etlbeexY4u9zb2Mv++BlvngO/9/Yct8dbnva/l6WftyWX/bwbFvLmW5ZaxvpnWNKG9fznbslnF+zPaZm+0YKmt9i/z7Vu16qqzPg7KOoSaeM6p2vihrn+rEG92eW1oghirb/vjKwto+7rglhbXdiSLXNeLA6zv2sNsLXS7l8bdttkznR/sydMcxRB0IxB0qK0Rk40QKnXtgyr0FL+F/F9x+PZT1eZBpMMtnbn/YzjSFfZl21D4Qb9iwodD2Z82a1fJ1QY2mKGtfLuvYFRKBqnONAdA/tQ/ENJsBDwAA6iDT9WOTCgwCMVREtuk9S0/417K7QEHKGsjKNIUYAOiN2gdiF9X9UdYFrr8vQGuZgnimdc2orL9vGcst63oKOLDaB+KymMoLVF22WQdAd1xj0GsGAKiDYr+4CgAAACpKhRiSM3rbH5kqL251AKBMmT5zy9KkWWi1D8TPbdpZaPuTTy+0eSqqyJAoIEIxsn2FGMBINemJwNBrtQ/EZclWAcm2vmXINpppMAvqyYBh8QyyACOV7fqxCAIxlSY00WtNmuIDFEswhe4oqFAHAjEAAKXJNPAgIEL1CMQdyla5zLa+mfhwbq6yplE5X9BrmQITUF9uKamn2gfi7e94qND2J8dJhbZPNQmJwEg5XwBAfdU+EAMA0D2VeHrNjCHqQCBmRFTi6TX7FAAAZROIqbSyQlORI5pGMwE4GNPw+8P9nkCEQNwx1S2gHaYiFi/bNi5rfYU16JwpxFA9AjEk58O5ucoKTAYMAYC6EIiBVMr6KiKay6ASdMesg+YyQEodCMRQEdmmewL1VMYzFgQm6M72x1cW2v5xxy0ptH0okkAMANCCIA7QfAIxAAAGAGgMs+5oh0AMLRR5z4v7XaBZ3CMHjJTzRbMZVKongbhmjHgBI+WDmSIYMASgSQRiIBWDSlA/nuQN0Jrrmu4JxJCc6VsAAGQlEAP0genLxbONAYB2CcQAAKSbmm4QDTo39rDby+5CzwyW3QEAAAAogwoxkIqKANRPWc86yFYxLYvtDJRJIAaAGjK401wCItRTkceu47Y4AjFUhItbes1FNQDAwbmHGAAAgJRUiAEA8L30DdakJwJDrwnEAH1g+nLxbGMAoF0CMQAA0Biey0I7BGJG5Jc/Xl9o+3MWLi60fQCA1/JEYHqtyNsO3HJQHIEYAKAF99QCNJ9ADJTCrIPiuZgH6sC5CiiTr10CAAAgJRViAIAKUTGl1zZs2FBo+7NmzSq0fSiSQAwANeRrpgDwRO3uCcQ1Y6cHaE1VDQBol0DcIQ8EotdUe/rDdgYAYC+BGAAAgBFr0n3ptQ/EKrVAHZjOC0A2rtOpA1+7BAAAQEq1rxADUC0qAgBAXQjEVJoLa4BqKfK87JwMxcj2LSUeoEk7BGJoIdMFn0EHAACyEogBAKAPDEJD9QjEANAFF7gAUF8CMQAA6RT5dXi+Cg/qw9cuAQAAkJIKcc14ah4AAEBvCMRQEUVO3YowfWsv25mmsC8DQPdMmQYAACAlgRgAAICUTJkGACrNV1sBUBSBGAAAoIY8cLd7AjEAQAsq01BPHjpIOwTimnGAAwAAZVp6wr+W3YWeEYgBAFARB1LylGkAAABSEogBAABIyZRpAABoME8i7o8ibztwy0FxBGKAhsp2P2C29QUAuicQAwBQGoNZQJkEYgCAChEQAfrHQ7UAAABISYUYoA9UfAAAqkcgZkRmf+TLZXcBAACgpwRiAACAGtr+jocKbX9ynFRo+1XgHmIAAABSUiHukCnEzebvCwAAzScQA6Uw6AAAQNlMmQYAACAlFWIAAKDnzAajDgRiqAjfU0tTuADqj7LOGf6+ADSJQAwAQDpFDioZUKLpntu0s9D2J59eaPP7qX0gdmIB2qESDwDN5rOedtQ+ENNsBjwAALqz/R0PFdr+5Dip0PahSAIxAHTBwB0A1JevXQIAACAlFWJITnULAICsBOKa8ZAAAACA3jBlGgAAgJQEYgAAAFIyZRoAqDTPOgCgKCrEAAAApCQQAwAAkJIp0wANZZop0A7nDCAjgRgAAKBLBpXqSSAGoBFciAAA7RKIAfpAWIP6cdwCNJ+HagEAAJCSCjEAAEAN/fLH6wttf87CxYW2XwUCMQAAACO2/R0PFdr+5Dip0PZfSyAGAKA07tUGyiQQAwBUiIDYH7YzEOGhWgAAACQlEAMAAJCSQAwAAEBK7iEGAAAaw/3htEMg7tCW9U8V2v6J/zKj0PYBAACyE4iBVIwaAwCwl0DMiGSriBe5vgdaV0ENaIdzBjBSv/zx+kLbn7NwcaHtQ5EEYgAAoOeyFVSoJ4G4ZlQEAAAAekMgBkph1BjqqYxbSgCgKAIxAD1lsAMAqAuBGAAA+sCAIVSPQAwAXXCBCwD1NVh2BwAAAKAMta8QG5mH7jiGgKpzngKgKLUPxDSbiyAAAKAopkwDAACQkgoxAI1gRgkA0C6BGAAAg0pASgIxAABAl4ocVDrQgNLsj3y5sGVm4R5iAAAAUlIhBuiDMqYimv4IAHBwKsQAAACkJBADAACQkinTAAAAjNgvf7y+0PbnLFxcaPuvJRADALTgPnyA5jNlGgAAgJRUiAEAKE1ZlfgyvjMWqB6BGAAAaAy3O9AOgRgAoEJczAP0j3uIAQAASEkgBgAAICVTpgEAoMFmf+TLZXcBKksgrhn3FQEQ4fMAAHpBIIaKcHHbH7YzAAB7uYcYAACAlFSIoYXFxz9SYOsnF9g2AAAwUirEAAAApCQQAwAAkJIp0x0qdkpthGm1AADAwXhYaPcEYgAA6INsBZVs60s9mTINAABASirEjIgRPgAiPIUfgIjZH/ly2V3oGYEYgJ4ygAYA1IVADABdMAAAAPUlEAP0gdAEnXP8AFAUgRgohQvc4tnGANA/nrFQTwIxleaCHgAAKErtA7HABADQPddUQEa1D8QAEOFiHgBon0AMAEA67vcEIgRiSE9VDQCArAbL7gAAAACUQYUYAABoDLPfaIcKMQAAACmpENeMES8AAIDeEIgBAFowCN0ftjNQJlOmAQAASEmFGAAAoIbMsOieQAwAUCEucOm1LeufKrT9E/9lRqHtQ5FMmQYAACAlFWIgFZUXmsK+DADdUyEGAAAgJYEYAACAlEyZhoow/REAAPpLhRgAAICUBGIAAABSMmUaAACAEWvSd1urEAMAAJCSCjEj8vC//Xuh7U+/dFGh7beryPWt2roCtMP5EYAmEYg7lC0gAgBAO1wvUwcCMQAA9IGACNXjHmIAAABSUiEGoKdUQACAuhCIgVIITQAAlE0gptKEJpqijH052/FT1vpm284A0CTuIQYAACAlFWIAAIAu+Z72ehKIAYBKMy0dgKIIxJCcC00AALISiAEAMEAKpCQQAwCQjvs9gQiBGAAAaBCzHWhH7QOxHR4AAIBO1D4QZ2MAAAAAoDcGy+4AAAAAlEGFGAAAoIbKmj26+PhHCl1uxMkFt/8/BGIAAGiwJoUX6DWBGACgBc/t6A/bGSiTe4gBAABISYUYAGpIVQ0AuicQAwBUiMEOgP4RiIFUXGgCALCXQAwVIagBAEB/eagWAAAAKQnEAAAApCQQAwAAkJJADAAAQEoCMQAAACl5yjQj8rMNFxfa/vRCWwcAAHg9FWIAAABSUiGGFoqsiKuGA3Xm/AhAk6gQAwAAkJIKMQAA0HOeQUMdqBADAACQkgpxh4x4QXccQwAAlE0gBuiDTAMAZa1rpm0MAPSGQAwAAH1g4A6qRyCG5Hw4Q3ccQwBE+Fq6uuooED/++OMHfG9wcDDe/OY3x+GHH95pn2AfF5rQOccPAMDBdRSIjz766BgYGDjozxxxxBFx2WWXxbJlyzrqGK25wAUgG599ABSlo0B86623xiWXXBJz586N888/PyZOnBjPPPNMfP/734+77rorPvvZz8auXbvi2muvjQkTJsSiRYt63W8AAADoSkeB+I477ojzzjsvvv3tb+/3+oIFC2LRokWxefPmWLduXRxxxBHxta99TSAGAKg4lXggo44C8caNG+OHP/xhy/c+9KEPxdlnnx0RESeeeGJcc801HXcOAADozsP/9u+Ftj/90moVvwzu0I7BTn5pwoQJsW3btpbvbdu2LYaGhiIiYteuXTF27NjOewcAAAAF6ahCfMEFF8Ty5cvjTW96U3z4wx+Ot7zlLfH000/H2rVr46qrropPfOIT8fzzz8dNN90Us2bN6nWfAQCgK74iBzrXpFkHHQXia665Jp555plYsmRJLFmyZN/rg4ODcfHFF8d1110Xd955Z2zZsiXuvffennUWAACAV5ke3r2OAvHo0aPjW9/6Vixbtizuu+++ePbZZ2Py5MkxZ86cmDZtWkREnH766fGnP/0pDjnkkJ52GAAAAHqho0C81zHHHBPHHHNMy/fGjx/fTdMjZlQEAACATnT0UC0AAACoO4EYAACAlLqaMg0AlMMtQ8WzjfvDdgbKpEIMAABASirEUBFGyPvDdgaqznkKoH9UiAEAAEhJIAYAACAlgRgAAICUBGIAAABSEogBAABIyVOmGZGbN60utP1Lj76y0PYBAAD+kQoxAAAAKQnEAAAApCQQAwAAkJJADAAAQEoCMQAAACl5ynTNeNozAGUq8nPIZxAA/SYQAwAAPaeQQx0IxB1ygDdbpgpIWfuyY4imsC8DQH0JxAB9kCk0ZVpXAKDePFQLAACAlFSIAQAAupTplrsmEYgBaARTtQFa+9mGiwttf3qhrUOxBGKAhhIQAQAOTiAGACrN4A5NYV+mKZo060AgptJ8cAAAAEURiAEAMAgNNeS47Z6vXQIAACAlgRgAAICUTJkGACAd3xnbXKYR0w4VYgAAAFISiAEAAEjJlGkAqCFTAgGgeyrEAAAApFT7CrERcgAAADpR+0AMAFAEg+79YTsDZRKIoSJcEAAQ4fMAoJ/cQwwAAEBKAjEAAAApmTINpGIqIgAAe6kQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQ0uiyO0A9/Mc/nVl2FwAAAHpKhRgAAICUVIgBAKDBbt60utD2Lz36ykLbhyIJxAB94LYDmsK+DECTCMQAAEDPGUBrribNOhCIa8aJhV6zTwEAkJVADC0IicWzjWkK+zIA1JdA3CEXQACtOT8CkJHPv3oSiAFoBBciAEC7BGKAhhIQaQr7MgBFEYgBAKAPDO5A9QjEAAAIa1BDjtvuCcRUmoMcAAAoikAMAEA6Bt2BCIEYAABoEIMdtEMgBoAacsFXPNsYoPkGy+4AAAAAlEGFGCpCJQIAAPpLhRgAAICUal8hVlUDAJok27VNtvUFqkWFGAAAgJQEYgAAAFISiAEAAEhJIAYAACCl2j9UC6AdHt4CAMBeKsQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkJJADAAAQEq+hxgAABrsP/7pzLK7AJUlEAMAADBiTRpkMWUaAACAlARiAAAAUhKIAQAASEkgBgAAICUP1aqbB75abPv//IVi2wcAAKgIgZiREcQBAICGEYghO4MdAAAkJRAD9IOBB5qiyH3ZfgxAnwnEADRDWYMOBjsAWnN+pAYE4k45wAFac34EoEyZPocyrWtBBGIAAIBuuaWklgRiaMUJrXhGNItnG9MU9mUACiIQA0AdCYkA0DWBGACglWyDDtnWtwy2MVTOYNkdAAAAgDKoEFNtmUZSM60rAJTN80KAUCEGAAAgKYEYAACAlEyZBgCokmy30GRbX6BSVIgBAABIqf4VYqOKAAAAdECFGAAAgJTqXyEGaIdZJQAA/J0KMQAAACkJxAAAAKQkEAMAAJCSe4gBAAAYuQY9k0WFGAAAgJQEYgAAAFISiAEAAEjJPcQAANBkDbrfE3pNhRgAAICUBGIAAABSEogBAABISSAGAAAgJYEYAACAlARiAAAAUhKIAQAASEkgBgAAICWBGAAAgJRGl90B2rP2tEcLbf+cQlsHAACoDhViAAAAUlIhBugDszuKZxv3R5Hb2TYGoN8EYkbEhSYAANA0AnGHBESawr4M3XEMAbSW7fxoBk09CcQAQKVlu6gGuuOcQTsEYmjBCB9N4IIAaIdzBpCRQAxAT7moBoD+8JnbPYEYKIUTOHTHMVQ82xig+QRiqAgXXgAA0F8CMQAA9IHBb6gegZhK88EBAAAURSAGACAd3ygBREQMlt0BAAAAKIMKMZCKafhA1WU7T2VbX6BaVIgBAABIqfYVYqOKAAAAdEKFGAAAgJQEYgAAAFKq/ZRpAAAA+qdJt60KxAAA0GBNCi/Qa6ZMAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASqPL7gDt+T//79eFtn/OrDMKbb8uitzOtjEUw/kRAGiXQNwhF140hX25P2xnmsKAIUBrzo/1JBAzIi7mAVpzfiyebQz1VNax65xBOwRiqAgnb3rNhUiz2c70WrZ9SjWPJsh23BZBIAZK4QQOQDY++6B6BGIAgBaEF4DmE4ipNBcjAABAUQRiAIAKyTYYnG19gWqpfSB2EgXa4ZwBAMBeg2V3AAAAAMpQ+woxAAAA/dOkGXcCMQAANFiTwgv0minTAAAApCQQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKo8ta8PDwcEREvPDCC121s+e//tqL7hzQgfpnuZZbp2Vabs7lZlpXy23GcjOtq+U2e7mZ1tVy+7PcTOt6sOV20sbe3HkgA8Nv9BMFefLJJ2PKlCllLBoAAIAEnnjiiZg8efIB3y8tEO/Zsyd27NgR48aNi4GBgTK6AAAAQAMNDw/Hiy++GJMmTYrBwQPfKVxaIAYAAIAyeagWAAAAKQnEAAAApCQQAwBvyB1WADSRQAwAr/HrX/86FixYEFOnTo1DDz00jj322Ljkkkvij3/8Y9ld64vHHnssBgYG4pZbbomIiJ07d8ZHP/rRuP/++8vtGAAUQCAGgL9buXJlzJ49O55++un44he/GOvXr4+lS5fGxo0b473vfW9s27at7C4W7m1ve1v84he/iDPOOCMiIrZu3Rq33XZb7Nmzp+SeAUDveco0AETEz3/+8/jABz4Qn/zkJ+Omm27a771nn302jj/++Jg4cWJs3ry5nA6WZOPGjTFv3ry47777Yu7cuWV3BwB6SoUYACLi+uuvj8MPPzyuu+6617131FFHxY033hhnn312vPTSS7F79+746le/GjNnzoxDDz00pk6dGkuXLo1XXnll3+9cdNFFcdppp8U3vvGNOPbYY+PQQw+NOXPmxKOPPhp33XVXzJw5Mw477LA44YQTYuvWrV3/3ty5c18XWDdu3BgDAwOxcePGiIi45ZZbYvTo0fGrX/0qZs+eHWPGjIm3v/3tccMNN+z7nddOmd4bhiMi5s2bF3Pnzo2VK1fGwMBAPProo/sta82aNTFq1Kh44oknOvwLAED/jS67AwBQtuHh4bjnnnvirLPOisMOO6zlz5x77rn7/v3xj388Vq9eHUuXLo2TTjoptmzZEp///Odjy5YtsWHDhhgYGIiIiAcffDB27NgRN954Y7z88suxaNGimD9/fgwMDMTVV18dY8eOjUsvvTQuuOCC+O1vf7uv/U5/byT27NkT5557bixZsiSuvfbaWLVqVVx55ZUxc+bMOPXUU/f72fe85z2xcuXKWLx4caxcuTLmzp0bkyZNiiuuuCJuu+22WLFixb6fvfXWW+ODH/xgTJkypa3+AECZBGIA0nvuuefilVdeiWnTpr3hz/7ud7+LVatWxRe+8IVYunRpRESccsopMWnSpFiwYEGsX78+5s+fHxERL774Ynzve9+Ld77znRERsWnTpvj6178eP/3pT+Pkk0+OiIjf//73ccUVV8TOnTvj8MMP7+r3RmJ4eDiWL18eCxcujIiIOXPmxA9+8IO46667XheIh4aGYsaMGRERMWPGjH3/Puecc2LNmjVx9dVXx8DAQDz55JNx7733xpo1a0bcDwCoAlOmAUhv9OhXx4d37979hj+7adOmiIg4//zz93v9vPPOi1GjRu2bnhwRMX78+H2hNiJi4sSJERFxwgkn7HttwoQJEfHq05y7/b2Rmj179r5/H3LIIXHUUUfFSy+9NOLfX7hwYTz22GP7njy9evXqGDduXJxzzjlt9wUAyiQQA5De+PHjY9y4cbF9+/YD/sxLL70Uzz//fPzlL3+JiIi3vvWt+70/evToOPLII/cLqENDQy3bGjt27EH70+nvjdQ/TgsfHBxs6ynSJ598ckybNi1Wr14dEa8G4vPOOy/GjBnTk/4BQL8IxAAQEaeeemrcd999+z0Y67W++c1vxpFHHrnv/0899dR+7//tb3+L5557br+f6aeBgYHXVbh37dpV2LIuuuiiWLt2bWzevDkeeeSR+NjHPlbIsgCgSAIxAETE5ZdfHn/+859j2bJlr3vvqaeeihtuuCFmzJixb1rwHXfcsd/PfPe7343du3fH+973vr709x8NDQ3Fk08+ud9rDzzwQNftjho1quXrF110UezcuTMuv/zymD59+n7TuQGgLjxUCwAiYtasWbFixYpYtmxZPPzww3HhhRfGkUceGb/5zW/i+uuvj5dffjl+8pOfxPTp0+PCCy+M5cuXx1//+td4//vfH1u3bo2rrroq5s2bF6eddlop/T/zzDNj3bp1sWTJkjjrrLPi/vvv3zeluRt7H9j1ox/9KMaPHx/vfve7IyJi6tSpccopp8Q999wTX/rSl7peDgCUQYUYAP7uM5/5TNx9990REfGpT30q5s+fH1/5ylfizDPPjK1bt8b06dMjImLVqlXxuc99Lm6//faYP39+rFy5Mi677LK4++67Y3CwnI/Wiy++OD796U/Hd77znZg/f348+OCDceedd3bd7rve9a44//zz4+abb44LLrhgv/fOOOOMGDVqVCxYsKDr5QBAGQaGh4eHy+4EAFA/p59+eowZMybWrl1bdlcAoCOmTAMAbVmxYkU88sgjcc899/TkPmUAKItADAC0Zd26dfGHP/whrr/++jjxxBPL7g4AdMyUaQAAAFLyUC0AAABSEogBAABISSAGAAAgJYEYAACAlARiAAAAUhKIAQAASEkgBgAAICWBGAAAgJQEYgAAAFL6b5xkRt21RnAbAAAAAElFTkSuQmCC",
      "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": 10,
   "metadata": {},
   "outputs": [],
   "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 = 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 = shuffle_matrix_rows(G.T).T\n",
    "    dg_m_shfls_col[rep] = np.mean(CalculateDs(ConvertToG(G_shuffled, relabus)))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkK0lEQVR4nO3de5CV9X348c+hy3VhQS5yF3dobIzZBE1UqkboKMEoFUlsUWMcLlZtm+AEW2f6i8o6NuNkUk0zqY7aOJVYXcMlUVGTRlOWRERMNJClK42GgCAgIKByvz2/Pwhb1xXZXZ7vnr28XjOO2XOe8+znwSdn9833nOcUsizLAgAAAMhdp2IPAAAAAO2V6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIpKdY3PnToUKxfvz569eoVhUKhWGMAAABAk2VZFu+9914MGTIkOnU6+np20aJ7/fr1MXz48GJ9ewBamcGDB8f1118f999/f2zYsKFIU5weEa9ExBkR8ZsizdBODI6I6yPi/ogo1n9OAGgBa9eujWHDhh31/kKWZVkLzlPnnXfeiT59+sTatWujrKysGCMAtC/LlkWMGROxaFHEqFHFnqaeO+889jaFwlvRteujsXfvVZFlA9MP9SE2buwUDz3UM6ZM2RGDBh0qygyNcl4j/kCLbOPOjfHQqw/FlFOnxKDSQcUeJzf/9Ll/KvYIALQS7777bgwfPjy2b98evXv3Pup2RVvpPvKS8rKyMtENkIeePf/v363sebVr18ZstTMiIrp06RkRxZm/S5cj/+7ZyJmLpLQ1D3dYl0OH/zC7dO8SXdvAvI3ldxYAPuhYb5d2ITUAAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACCRon1ONwB0dNVR2azHjc11CgAgJSvdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkEhJsQcAAGgrKqsrm/6YsU1/DADth5VuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCIlxR4AgHz99Kc/jW3/+7/HtY8rr7wyp2kAADo2K90AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREqKPQAAfJSamppc9lNRUZHLfgAAmsJKNwAAACQiugEAACARLy8HgBxUR2WxRwAAWiEr3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASKSk2AMAQEuoqak55jbbtnWLiI/Fa6+9Flu27Glwf0VFRYLJAID2zEo3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEjEhdQASKKqqqruf9fUfPKY23fvvi0+9rGI119/LXbv3pJyNACAFmOlGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEVcvB6CB9195HACA5rPSDQAAAIlY6QZoJY53dfmEP/whLsppFgAA8iG6AQASqqyubN7jxjbvcQC0Ll5eDgAAAImIbgAAAEhEdAMAAEAiohsAAAAScSE1AGikmpqao9731sC3Gr2fgScOzGMcAKANsNINAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIpKfYAAG1dVVVVsUcAAKCVstINAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBEfGQYAHzAH/rNjQ0D1ibb/1ub3jqux9fUbIqKioqcpgEAUrLSDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJFJS7AEAgKapffXEePvto98/dmyLjQIAHIOVbgAAAEhEdAMAAEAiXl4OdFhVVVXFHgEAgHbOSjcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkUlLsAQCApntr01tHva+mZlOj91NRUZHHOADAUVjpBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIiXFHgAAgIYqqyub/pixTX8MAGlZ6QYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASMTVywHIxbx5nyz2CAAArY6VbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBGf0w20SVVVVcUeAQAAjslKNwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEvGRYQDQgdXU1OSyn4qKilz2AwDtjZVuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBEGh3dd999d7z11lspZwEAAIB2pdHR/Y//+I+xZs2auq+zLIu/+7u/i3Xr1iUZDAAAANq6Rkd3lmX1vj506FDcf//9sWnTptyHAgAAgPbguN7T/cEQBwAAAP6PC6kBAABAIqIbAAAAEmlSdBcKhUbdBgAAAESUNGXj0aNHN7jts5/9bIPbCoVCHDhwoPlTAQAAQDvQ6OieNWtWyjkAAACg3RHdAAAAkEiTXl7+fmvXro2NGzdGoVCIYcOGxaBBg/KcCwAAANq8Jl1Ibd++fXHnnXfGiBEj4uSTT47Ro0fH2WefHUOHDo0/+7M/i+9+97tx8ODBVLMCAABAm9Lole49e/bEuHHjYvHixTF69OiYOnVqDBo0KLIsi3Xr1sVzzz0XX//61+Opp56KZ555Jjp37pxybgA4ptqB9zZp+z2dhkfE/0szDADQITU6ur/97W/Hr3/963j88cfj0ksvbXD/N7/5zXj66afjr//6r+Pf/u3f4utf/3qugwIArVdNTU29r7cd2BYREa+9/lpsKdnSqH1UVFTkPhcAFFujX14+d+7cmDlz5ocG9xGXXHJJzJgxI6qqqnIZDgAAANqyRkf3qlWr4vzzzz/mdmPGjImVK1ce11AAAADQHjQ6unfv3h29e/c+5nZ9+vSJnTt3HtdQAAAA0B40OrqzLItOnY69eaFQOK6BAAAAoL1o0keGCWoAAABovEZfvTwi4rLLLouuXbt+5DZ79+49roEAAACgvWh0dF9zzTVWugEAAKAJGh3dDz30UEREvPTSS7FmzZoYOXJknHHGGanmAgAAgDav0dG9ffv2mDBhQixZsqTutnPOOSceffTRGD58eJLhAAAAoC1r9IXUbrnllnjllVfi9ttvj6effjruuuuuWLlyZVx//fUp5wMAAIA2q9Er3QsWLIg777wzbrzxxoiIuOiii2Lo0KFx1VVXxc6dO6O0tDTZkAAAANAWNXqle+PGjfGZz3ym3m1jx46NgwcPxhtvvJH7YAAAANDWNTq69+/fH126dKl3W9++fSMiYs+ePflOBQAAAO1Ao6P7o2RZlsduAAAAoF1pUnQf7XO6fX43AAAANNToC6lFRPzt3/5tlJWV1X19ZIX7uuuui169etXdXigU4uc//3lOIwLtSVVVVbFHAACAFtPo6D7//POjUCg0eCn5mDFjIqL+S8y93BwAAACaEN3V1dUJxwAAAID2J5cLqQEAAAANiW4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACCRkmIPAABAPiqrK5v3uLHNexwAx2alGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBESoo9AAAcS+3Ae4s9AgBAs1jpBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkEhJsQcA2oaqqqpijwC0czU1Nbnsp6KiIpf9AEAerHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIpKTYAwAA+ap99cSj3veJUze14CQAgJVuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiZQUewAAAIqrsrqy6Y8Z2/THAHREVroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgkZJiDwBA2zBv3ieLPQIAQJtjpRsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImUFHsAADqW2oH3fujtJ5SUxMeif6zqOze2HTjQwlMBAKRhpRsAAAASEd0AAACQiOgGAACARLynGwBoV2pqanLZT0VFRS77AaBjs9INAAAAiYhuAAAASER0AwAAQCKiGwAAABJxITXoAKqqqoo9AgAAdEhWugEAACAR0Q0AAACJiG4AAABIRHQDAABAIqIbAAAAEnH1cgDoQGpfPfEj7//EqZtaaBIA6BisdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQSEmxBwCgbaodeG+xRwAAaPWsdAMAAEAiohsAAAASEd0AAACQiPd0AwB8iJqamuPeR0VFRQ6TANCWWekGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEhEdAMAAEAiohsAAAASEd0AAACQiOgGAACAREqKPQBwdFVVVcUeAQAAOA5WugEAACARK90AADRZZXVl8x43tnmPA2irrHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJiG4AAABIpKTYAwBQfPO2zDvmNrUDf9ECkwAAtC9WugEAACAR0Q0AAACJiG4AAABIRHQDAABAIi6kBgCQSE1NTS77qaioyGU/ALQ8K90AAACQiOgGAACAREQ3AAAAJCK6AQAAIBHRDQAAAImIbgAAAEjER4YBANBiKqsrm/e4sc17HECxWekGAACAREQ3AAAAJCK6AQAAIBHv6QagTu2rJxZ7BACAdsVKNwAAACQiugEAACAR0Q0AAACJeE83AEArV1NTc9z7qKioyGESAJrKSjcAAAAkIroBAAAgEdENAAAAiYhuAAAASER0AwAAQCKuXg6JVFVVFXsEAACgyEQ3QDvz3Pbn4g9bTij2GAAAhJeXAwAAQDKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENAAAAiYhuAAAASKSk2AMAAK1H7asnfuT9nzh1UwtNAgDtg5VuAAAASMRKNwAArV5ldWXTHzO26Y8ByJuVbgAAAEhEdAMAAEAiXl4OH1BVVVXsESCZY10kC2i/ampqctlPRUVFLvsB6CisdAMAAEAiohsAAAASEd0AAACQiOgGAACAREQ3AAAAJOLq5QCt1Lwt85q0ffn2bXFRolkAAGgeK90AAACQiJVuAADapcrqyuY9bmzzHgfwYax0AwAAQCJWugEAaLSamppc9lNRUZHLfgBaOyvdAAAAkIiVboB25g9/6Bu1OwYUewzaqdpXT2zUdntKChH9Dp+PGw4cPh8/ceqmlKMBQKtkpRsAAAASsdJNu1FVVVXsEQAAAOqx0g0AAACJWOkGAKDF5XEVdFdAB9oC0Q3QBPO2zGvyYy7vf3mCSQBIpbK6sumPGdv0xwAdg+gGSKw5oQ7t0Udd+dyVzWkOnxkOtAWim1bBRdAAAID2SHRzXMQyAEB+vxNdeeWVuewHaD1EN0Abc9SX6G4utOwgANRp7luJWvK6H96rDsUhujswq9S0B+3x/dIf9b5XaK+Odd57zzcp5fXe8Ob44M+xed879s+1xoS6FXNoPYoW3VmWRUTEu+++W6wROrxdu3YVewQ4bvt37y/2CLk7uG9fsx63f//+ePeP/27uPorp4B9/Lhzcvz8O7i/Of9cj3/dgG/0zbE0OZu/7s9x//H+W7fH/69Bcjfkd7sN+x967c2+Tv1cev6vPnTv3uPcREfFXf/VXuewH8nLk/x9H2vZoCtmxtkhk3bp1MXz48GJ8awAAAMjF2rVrY9iwYUe9v2jRfejQoVi/fn306tUrCgXvQwQAAKDtyLIs3nvvvRgyZEh06tTpqNsVLboBAACgvTt6jgMAAADHRXS3UXv27Inp06dHnz59YvDgwXHXXXcdddvf/OY3cfbZZ0ePHj3izDPPjJdffvlDt5s7d66X+ncgeZ5Dffr0iUKhUO+fHTt2pD4EiijP82fevHlxyimnRGlpaXz+85+PNWvWpB6fViCvc+iDzz1H/vnBD37QEodBkeR1/mRZFpWVlTFs2LA44YQTYvLkybF58+aWOASKLM9z6F/+5V+ivLw8+vTpE1OnTvU7EA1ltElf/epXs0996lPZyy+/nP3oRz/KevXqlc2dO7fBdjt27MgGDRqU3XTTTVltbW02Y8aMbODAgdmOHTvqbbdt27Zs0KBBmVOi48jrHFq3bl0WEdnvf//7bMOGDXX/HDp0qKUPiRaU1/mzePHirKSkJLvvvvuylStXZhMmTMhGjx7d0odDEeR1Dr3/eWfDhg3ZzTffnI0YMSLbvn17Sx8SLSiv8+e+++7Lhg0bllVXV2c1NTXZeeedl1166aUtfTgUQZ7nUGlpafboo49mK1asyMaNG5f95V/+ZUsfDq2cwmqDduzYkXXr1i1buHBh3W133HFHNmbMmAbbPvjgg1l5eXldAB06dCj70z/90+w//uM/6m137bXXZueee67o7iDyPIeeffbZbPDgwS0wNa1FnufPpEmTsilTptRtv2rVqmzEiBHZ5s2bUx4CRZbi51iWHT5/unXrlj377LOJJqc1yPP8ufTSS7Obbrqpbvsnn3wyKy0tTTk+rUCe59Bpp52W3XbbbXXbr1+/PisUCtnKlStTHgJtjJeXt0HLly+P/fv3xznnnFN323nnnRdLly6NQ4cO1dv2xRdfjPPOO6/uZeOFQiHOPffcWLJkSd02ixYtiurq6vjGN77RMgdA0eV5DtXW1sYpp5zScsNTdHmeP9XV1fHFL36xbvvy8vJYvXp19O/fvwWOhGLJ++fYEbfddltccMEFceGFF6Y9AIoqz/OnX79+8fTTT8ebb74Zu3fvjqqqqjj99NNb7mAoijzPoVWrVsXZZ59dt/3gwYNjwIABH/ocRcclutugDRs2RP/+/aNLly51tw0cODD27NkTb7/9doNthwwZUu+2gQMHxrp16yIiYu/evXHdddfFPffcE927d08/PK1CnufQq6++Grt27YqxY8fG4MGD4+KLL47f/e536Q+Cosnr/Nm+fXts27YtDhw4EOPHj49BgwbFxIkT480332yR46B48nwOOuKNN96IRx99NG699dZ0g9Mq5Hn+3HbbbVFSUhLDhg2LXr16xS9/+cuoqqpKfxAUVZ7n0MCBA+v93Nq5c2ds3bo1tmzZkvAIaGtEdxu0a9eu6Nq1a73bjny9d+/eRm17ZLs77rgjzjjjjPj85z+fcGJamzzPoZUrV8bWrVvjlltuiSeeeCK6d+8eF1xwQbz33nsJj4Biyuv8OXKhmRkzZsTVV18dCxYsiL1798aECRMarDTQvuT5HHTEgw8+GJ/97GfrrTjRPuV5/qxevTp69OgRCxYsiEWLFsWwYcNi2rRpCaenNcjzHJo8eXLceeed8eqrr8aePXti5syZERGxb9++VOPTBpUUewCarlu3bg2eEI583aNHj0Zt26NHj1ixYkU88MADUVNTk3ZgWp28zqGIiJ/+9Kexf//+6NmzZ0REPPLIIzF8+PBYsGBBXHXVVakOgSLK6/wpKTn8I+jaa6+Nr3zlKxFx+PwZOHBgvPjii/Ve9kf7kudz0BHz5s2LG264IcG0tDZ5nT9ZlsU111wT3/72t2PChAkRETFnzpwYMWJELF261F/gtGN5PgfdeuutsWrVqjjttNOic+fOcf3118eoUaOirKws4RHQ1ojuNmjo0KGxZcuWOHDgQN0vrRs3bozu3btHnz59Gmy7cePGerdt3LgxBg8eHPPnz4+tW7fGyJEjIyLi4MGDERHRs2fPuP/+++PLX/5y+oOhKPI6hyIO/23v+/8GuFu3blFeXu4lwu1YXudP//79o3PnzvHxj3+87r5+/fpFv379Yu3atcmPg+LJ8zkoImLt2rVRW1sbEydOTD47xZfX+bN58+ZYu3ZtfPrTn667b/jw4dG/f/9Ys2aN6G7H8nwOKi0tjTlz5sQ777wThUIhysrK4sQTT4yTTz65JQ6FNsLLy9ugUaNGRefOnePFF1+su+3555+PM888Mzp1qv+fdPTo0fHCCy9ElmURcfizBBcvXhyjR4+Or33ta7Fy5cpYtmxZLFu2LL7//e9HRMSyZcvi0ksvbbkDosXldQ5lWRYjR46Mhx56qG77nTt3xmuvvVYvpGhf8jp/SkpK4jOf+UwsX768bvstW7bEli1b/LLSzuV1Dh2xdOnSGD58eJx00kktcwAUVV7nT9++faNr165RW1tbt/2WLVvi7bffjvLy8pY5GIoiz+egm2++OWbPnh29e/eOsrKy+NWvfhXvvPOOV2tRXzEumc7xu/7667PTTjste+mll7If//jHWVlZWTZ//vwsyw5/ZumuXbuyLMuyd955JxswYEA2Y8aM7H/+53+yGTNmZIMGDWrwOd1ZlmULFy70kWEdSF7n0Ne+9rXspJNOyhYuXJitWLEimzRpUvbJT34yO3DgQNGOjfTyOn/mzp2blZaWZnPmzMlqa2uzCRMmZGeccYbPee8A8vw5NmvWrGz8+PFFOQ6KI6/z54YbbsjKy8uzRYsWZTU1Ndn48eOzc845x3NQB5DXOfSd73wnKy8vz1566aXs17/+dXbKKadkN954Y7EOi1ZKYbVRO3fuzK655pqstLQ0GzJkSPad73yn7r6IqPf5pUuXLs1OP/30rFu3btlZZ52VvfLKKx+6T9HdseR1Du3evTubOXNmNnjw4KxHjx7ZhAkTsjfeeKMFj4RiyPM56IEHHshGjBiRde/ePfvCF76QrV27toWOgmLK8xy64YYbsiuuuKKFJqc1yPNn2E033ZQNHTo069u3bzZ58uRs06ZNLXgkFEte59CBAweyG2+8MevXr182YMCAbObMmdn+/ftb8EhoCwpZ9sfXSgAAAAC58p5uAAAASER0AwAAQCKiGwAAABIR3QAAAJCI6AYAAIBERDcAAAAkIroBAAAgEdENADRLlmXFHgEAWj3RDQA5GjduXPTt2zf27dt31G0qKiri/PPPb9b+V69eHYVCIR566KHcH3PyySfHlClTIiKiuro6CoVCVFdXR0REZWVlFAqFum0XL14cl1xySROnB4COp6TYAwBAezJt2rR47rnn4plnnonLLruswf2vvPJKrFixImbPnt2s/Q8ePDiWLFkSI0eOPM5JG/rxj38cZWVlH3rftddeGxdddFHd1//+7/8etbW1uc8AAO2N6AaAHE2aNClOOOGEeOSRRz40umfPnh1lZWVx+eWXN2v/Xbt2jdGjRx/nlB/u9NNPP+p9w4YNi2HDhiX5vgDQnnl5OQDkqFu3bnHVVVfFU089Fe+++269+/bv3x9VVVVx5ZVXxq5du+Lv//7vY8SIEdGlS5fo27dvTJo0KVavXl23/dixY+Pqq6+Oyy+/PEpLS2PcuHEf+lLxX/ziFzF+/Pg44YQTokuXLlFeXh6VlZVx6NChet//zTffjAkTJkT37t1j+PDhMWvWrDh48GDd/e9/efkHvf/l5VOmTInZs2fHmjVr6mY588wz45xzzmnwuAsvvDDGjRvXxD9FAGg/RDcA5GzatGmxZ8+emD9/fr3bf/KTn8TmzZtj+vTpcckll8TPfvaz+Na3vhU/+9nPorKyMn7+85/HDTfcUO8xP/zhD6NXr17x5JNPxs0339zgey1fvjwuuOCC6N+/f/zwhz+MBQsWxOc+97m4/fbbY86cOfW2nTVrVgwcODAef/zxmDZtWnzzm9+Mf/iHf2jy8d16661x8cUXx6BBg2LJkiVxySWXxPTp02PJkiXx+uuv1223du3aWLhwYUydOrXJ3wMA2gsvLweAnJ1xxhkxatSoeOSRR+oF5+zZs+NTn/pUDB06NEpLS+Ouu+6K8847LyIOr2q//vrr8cADD9TbV5cuXeK+++6Lrl27RkTUWwmPiPjtb38b48aNi4cffjg6dTr8d+njxo2LJ598Mqqrq+OKK66o2/aiiy6KBx98MCIixo8fH++++27ce++9ceutt0bfvn0bfXwjR46MAQMG1Hup+5VXXhkzZ86Mhx9+OG6//faIiHj44YejV69eMWnSpEbvGwDaGyvdAJDAtGnTYuHChbF+/fqIiNi6dWs89dRTMX369BgyZEj893//d5x77rmxevXqePbZZ+N73/teLF68OPbu3VtvP6eeempdcH+Yr3zlK/HMM8/Evn374re//W3Mnz8/Zs2aFQcOHGiwr8mTJ9f7+otf/GLs27cvXnzxxeM+3t69e8eXvvSl+M///M+622bPnh2TJ0+O7t27H/f+AaCtEt0AkMCXv/zl6Ny5czz22GMREfHYY49FoVCIq6++OiIiHnnkkRgxYkSUl5fHFVdcEU888UT06NGjwX569uz5kd9n9+7dce2110bv3r1j1KhRcfPNN8fq1aujc+fODT5He9CgQfW+PvHEEyMiYtu2bc0+zvebPn16rFq1Kp5//vlYsmRJ/O53v/PScgA6PNENAAn07ds3Lrvssnj00UcjIuIHP/hBTJo0Kfr27RvPP/98XHPNNfGlL30p1q1bF2+//XY899xz8ed//udN/j433nhjzJs3L+bMmRPvvfde/P73v4+HH344Onfu3GDbrVu31vt648aNEfF/8X28xowZEyNHjow5c+bEnDlz4uMf/3iyK60DQFshugEgkWnTpsXLL78cixYtiqVLl8b06dMjIuKFF16IQ4cORWVlZQwdOjQiIg4ePBjPPvtsRESDq45/lOeffz7+4i/+IiZOnBilpaUREfHyyy/H5s2bG+zn6aefrvf1Y489Fj169Iizzz67ycf2J3/yJw1uKxQKMXXq1Hj88cfjySefPOqV0AGgI3EhNQBI5MILL4yTTjoprrvuuigvL48LLrggIiLOOuusiIj46le/GtOmTYutW7fGPffcE8uXL4+IiJ07d0avXr0a9T3OOuusmDNnTtx3331x6qmnxvLly+Of//mfo1AoxM6dO+ttO3/+/Bg6dGiMGzcu/uu//ivuv//+uOOOO6KsrKzJx9anT59466234ic/+UmMGjUqBg8eHBGHP05s1qxZEXH4/eYA0NFZ6QaARDp16hRTpkype2/zkc+5Hjt2bNxzzz3xwgsvxBe+8IWYOXNmnHTSSfGjH/0oIiJ++ctfNvp73H333TFp0qS45ZZbYsKECfH9738/brnllvibv/mbWLJkSb3P4f7ud78bv/rVr+Liiy+OOXPmxL/+67/GN77xjWYd29SpU+Pkk0+OiRMnxuzZs+tuHzp0aHz605+O8ePHx5AhQ5q1bwBoTwrZB6+yAgDQTG+++WaMGDEi5s+fHxMnTiz2OABQdKIbADhuy5YtiyeeeCLmz58fBw4ciBUrVtR9bjgAdGR+GgIAx23Pnj1x9913x4EDB+Kxxx4T3ADwR1a6AQAAIBF/DQ0AAACJiG4AAABIRHQDAABAIqIbAAAAEhHdAAAAkIjoBgAAgERENwAAACQiugEAACAR0Q0AAACJ/H/qPXGB/puLEQAAAABJRU5ErkJggg==",
      "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()"
   ]
  }
 ],
 "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
}