Hadza_total_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",
    "mm = relabus.shape[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30, 101)"
      ]
     },
     "execution_count": 13,
     "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_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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGI0lEQVR4nO3de3xV5Z32/2/Ap40yKvyKlUOwMBAsKGPSymMCRSXIgBHRTFtfBQRMIzAYxDHqiCdAsaBotS2kFiGjBVFrh0NRw2HKUU5OgaAcAqgcI4mIEgVa5qeQ5w9HFLs3ZGfnWvc6fN5/tQv8XncQMNdea913Sk1NTY0BAAAAAIB618D1AgAAAAAACCtKNwAAAAAAIpRuAAAAAABEKN0AAAAAAIhQugEAAAAAEKF0AwAAAAAgQukGAAAAAECE0g0AAAAAgMhZroJPnDhh+/fvt3PPPddSUlJcLQMAAAAAgITV1NTY4cOHrUWLFtagQfz72c5K9/79+61Vq1au4gEAAAAASNq+ffssLS0t7o87K93nnnuumX2xwPPOOy+pWeeff359LCmuTz75xHfZ6lyX2VH8ml1m8zX7IzuKv95h/ppdZvM1+yM7ir/eYf6aXWbzNfsjO4q/3mH+mk+XnYhPP/3UWrVqdbLbxpNSU1NTk3RaHXz66ad2/vnn2yeffJJ06VY/nn66XyJX2V48ku8qO4pfs8tsvmZ/ZEfx1zvMX7PLbL5mf2RH8dc7zF+zy2y+Zn9kR/HXO8xf8+myE1HbTstGagAAAAAAiFC6AQAAAAAQoXQDAAAAACBC6QYAAAAAQITSDQAAAACACKUbAAAAAAARSjcAAAAAACKUbgAAAAAARCjdAAAAAACIULoBAAAAABChdAMAAAAAIELpBgAAAABAhNINAAAAAIAIpRsAAAAAABFKNwAAAAAAIpRuAAAAAABEKN0AAAAAAIhQugEAAAAAEKF0AwAAAAAgQukGAAAAAECE0g0AAAAAgAilGwAAAAAAEUo3AAAAAAAilG4AAAAAAEQo3QAAAAAAiFC6AQAAAAAQoXQDAAAAACBC6QYAAAAAQITSDQAAAACACKUbAAAAAAARSjcAAAAAACJnuV4AEARj3hzjegkAAAAAAojSDdTCvy6/Rh/yf/URAAAAALzF4+UAAAAAAIhwpxuohRe/myrPKJInAAAAAPAapRuohZvuz9OHDN6nzwAAAADgKR4vBwAAAABAhNINAAAAAIAIpRsAAAAAABFKNwAAAAAAIpRuAAAAAABEKN0AAAAAAIhwZBhQCz8a9Dt5xm55AgAAAACvUbqBWth2diPXSwAAAAAQQDxeDgAAAACACKUbAAAAAAARSjcAAAAAACKUbgAAAAAARCjdAAAAAACIsHs5UAupo+/2IGWdBxkAAAAAvMSdbgAAAAAARCjdAAAAAACI8Hg5UAtLlkyUZ+TkyCMAAAAQML/ctMb1EpAkSjdQCzX2PddLAAAAQARtn/ScNmBKlnY+eLwcAAAAAAAV7nQHFI+ZAAAAAID/UboDSv6YiRmPmnzNwtlb5Bk9ctrKMwAAAAB4i9IN1MKfs2Z4kNLXgwwAAAAAXqJ0B9TPeo90vYRIuf3iX7peAgAAQOS9v3Ky6yV4bsqUH7peApJE6Q6osk//Js/oLk8AAAAAcDovbf8v6fx+Fw+Vzq+LsO1fRekGAAAAAPhG2PavonQDAAAAAHwjbI/UU7oD6o6PnvEgpcSDDAAAAADxVK69VxtwsXZ8XSyd01We0T1PHnESpTugRs7VPw5RXCSPCIz1FUflGfmd5REAAAAImHYLB2oDBpdr59dB2PavonQHVHGOFxseDPEgIxgaT/Lgrn/ek/oMAAAABMrCCu2dMA6t1aN0B9Srt0yUZ1wvTwiOERPbuF4CAAAAgACidAdU5aLz9SH+Oz0AAAAAQMgVDd7gQcrlHmR8gdIN1MahJR6EFHqQAQAAgCB5vWsL6fxi6fS6KbwyRZ5RvEIecRKlGwAAAADgG2Hbv4rSHVCDqtq7XgIAAAAA1LuKqWnyjLSx8oiTKN0AAAAA4FPbzm7kegmeK1s1Up6hr/VfoXQH1Jwtw+UZ/cx/Z/YBAAAACLewbRpN6Q6ofn+I3ideLlW2yZVnNJMnAAAAAPBaA9cLAAAAAAAgrLjTDQAAAAA+lTr6bnHCOvF8cKcbAAAAAAAR7nQH1JQb/k2eMexP8ggAAAAAp/F+yw+k81u+Lx1fJ2E7HpnSDQAAAAA+dXTpYtdL8FzYTmqidAfU6/kvyTOG2c3yDAAAAADxLV6sLd3t24frrrIfUboBAAAAwKfaX/zPrpfgubAdj0zpDqhnCt/Wh9yojwAAAAAQ38LZW6Tze+S0lc4HpTuwxlzbS54xTZ4AAAAA4HT+nDVDnNBXPB+U7oDq1/8+10sAAAAAILbuQJXrJSBJlG4AAAAA8Knfd5spnZ8vnV43YTsemdIdUAu+/aw8o4c9Ls8AAAAAEN+kBe9K5+d3vkg6H5TuwFq8a6c+pKs+AgAAAKit5U2XSef3s0Lp/Loo+dD1CrwXtuORKd0BxbsdAAAAAMLoQNlIfciN+ogvUboD6s2+r8ozsuUJAAAAQO11n3eHNuAe7fi6KK8ol87PtBzp/LqYNe1WfcjD+/QZ/4vSHVDbDn0qz8i2xvIMAAAAAPFtrd7legmeC9vxyJTugFJvqGDGpgoAAACAayN+2sH1EjwXtuORKd0BNSutueslRErlkTJ5RqY8AQAAINhe/G6qdH6RdDqiitINAAAAIBDeKxHfCBl8uXY+aiVsxyNTugEAAAAEwn3vPSJOGCKej9oI2/HIlG4AAAAAgfCjQb+Tzt8tnY7aCtvxyJRuAAAAAIFwf5v3XS8BHgjb8ciUbsDnfrlpjeslAAAA+MKgqvaulwAPhO14ZEp3QFWXVepD8qN3PIEfFVy2QBtwPEs7HwAAAIgwSndAlVeUyzMyLUeeAQAAANRW6ui7xQnrxPNRG+srjsoz8jvLI06idAO1MPuVDfKM3Ax5BAAAAOB7d1W7XkH9onQDAAAACIT3W34gnd+SfdogQOkGAAAAEAjzHrxfOn+4dDqiitIN+Fzv2zKl89kbHQAAANChdAdUTs9U10uARx7o1dD1EgAAAADUEaUb8LkfNNnjegkAAAAA6ojSDfjc8qbLpPP7WaF0PgAAABBllG4AAAAAgG9Ul1XqQ/I76DP+VwPPkgAAAAAAiBjudAMAAAAAfKO8olyekWk58owvUboBn6tce6824GLteAAAACDKKN0BVfJiiTzjgeEF8gyc2fZ/3CFOuFw8HwAAAIguSjdQC9npM10vAQAAAIiErdW7XC+hXlG6AQAAAAC+MeKn3u0s7gV2LwcAAAAAQIQ73QFV8PQA10sAAAAAAJwBpRsAAABAILS/+J9dLwFIGKUb8LnrnhmnDejWXzsfAACgniycvUU6v0dOW+l8RBOlG/C5vi/e6XoJAAAAvvDnrBnihL7i+YgiSjfgc4VXpkjnF6+QjgcAAAASUvJiiTzjgeEF8owvUboBnyvOGSpOGCKeDwAAAEQXpRsAAABAINx+8S9dL8Fz6ru+Xt7xjSpKN+Bzl7e5VDp/nXQ6AAAAEG2UbsDnxlw1yPUSAAAAfGHSgnel8/M7XySdj9opeHqA6yXUK0o34HO3T+konX/9BOl4AACAejMrrbnrJQAJo3QH1aElHoQUepABAAAAAOFF6QZ87v4277teAgAAAIA6auB6AQAAAAAAhBV3ugGfG1TV3vUSAAAAANQRpRvwuTlbhkvn97Ny6XwAAAAgyijdAbW1olKe0UyeAAAAAADhRukGAAAAAJ8K25nVUUTpDqiysmPyjJx8eQQAAAAAHxo6dKi78JAdj0zpDqickZyhDQAAAEBj4rTm+pAp+gg/oHQDAAAAgF/J7/rGvpnX+7ZMca7ZGnmCP1C6AZ/r94dGrpcAAAAAR9QbKMfbPHnNpA+kuVFC6QYAAAAAnOLE8o/kGQ2ukkf4AqUb8Lux68XzteMBAAAQPKWHO8kz+sS5HrbjkSndAAAAAOBTTboVOMn9QZM9TnLDiNINAAAAADhFiwumepASexO3sB2PTOkGfO79qWnS+S3HSscDAAAggF6qSZdn9ItzPWzHI1O6AQAAAARCdZn4Xd/8Dtr5iCRKN+BzR5cudr0EAAAAOFJ5pEw6X38aNyjdAAAAAIBTVK69Vx9ysT7CDyjdAAAAAAKhvKJcOj/TcqTzEU2U7oBSP2ZixqMmAAAA8JecnqmulxAZd3z0jAcpJR5kuEfpBmqhtHSVPKMgv03M6/sqGkpz27eXjgcAAEAA3da5hzxjijzBHyjdQC1UVm51vQQAAADAM9c9M04f0q2/PsMHKN2Azy2cvUU6v0dOW+l8AAAABE/fF+90vYTQoHQDPvcvo893vQQAAAA4MvuVDdL5uRmxr//l1W7SXDOzztfLI3yB0g0AAAAAOEVVp9c8SOkQ82rYNo2mdAMAAAAATlG5yIOnLYfqI/yA0g0AAAAAOMWgKo64qS+UbqAWhjw3wvUSAAAAAM+kjr7bg5R1HmS4R+kOqPJNlfKMeJsqRNF9FQfkGfnpLeQZAAAAALxF6U7S0gEXOsktL9eXbgAAAABAcijdSVo28l+l86+WTkdtzaz5sQcp73mQAQAAAJzZ+y0/kGe0fF8e4QsNXC8AAAAAAICw4k434HOz3imVzs/67lXS+QAAAECUUbqTNPSzC1wvAQAAAADgU5TuJL387hXS+UVdpeMRAIt37dQG8HsMAAAA3/DXZWWulxAalO4kvVci/s04+HLtfAAAAAD4hr37PpFnpKc3jXl99isb5NleHo9M6U5Scc5QccIQ8XzURqNzZjrLXtepr7NsAAAARFOjSys8SGnrQYZ7lO4ktf6f16Tzd0unIxB2vqmdf9lA7XwAp1g64ELXSwCAwCp5sUQ6/4HhBdL5iCZKN+BzH7w0XTq/ed5k6XwAp2o//CHXSwAAAB6idCfp/jYROdE94tLmfqYPyYp9ecE9m6Wx+dLpAL6JDTgBAIgWSneSBlW1d70EAECAFA1Wbw7DBpxAFFQ3GON6CfBIdrq7vYVQPyjdSUodfbc4YV3Mq/zhAwAAiK7et2XKM9bIE4Lj+LLxrpeAAKN0AwDgocIrU6Tzi1dIxwNAJP1qd095RtFV8gg4QulO1tj14vna8QAAAEBQFDw9wEnu9n/c4UEKrweFFaUb8Ln1FUel8/M7S8cDAADUn0NLxAGF4vnBMeudUnlG1nejcXuf0h1QpaWr5BkF+W3kGQAQNa93bSGdXyydDgCIisf/4VLXSwgNSneSdvTfLp3P3ugAEC7bzm7kegkAEFiVbXKl85tJpwfMzjf1GZcN1Gf4AKUb8LnVmyu1AXkdtPMBnMLVqRcAwmV8TjvXS0DY7XrL9QpCg9INAICX2IATABAAI9s8K8+YLE/wB0p3kioqKqTz27fnAfOom5XW3PUSANSjrJoF0vlrpdMBAFEhf9rSLDJPXFK6AQDw0PQBbFIJIHlVHcd6kPJHDzKA8KN0AwDgIZ6QAgAkQn1qEScW6VG6AQDw0ILZR6Tzc3Kk4wH4ROXae/UhF+sjgCigdAMAAAABUzR4gwcpl3uQAYQfpRvwueoy8SYW+dHYwALwi5p+2scEzfqK5wMAoJWdPtP1EuoVpTtJNfY910tAyJVXlEvnZxrPogJeWrxrpzagq3Y8AABIDKUb8Lmt1btcLwEAAPhM4ZUp8oziFfIIIBIo3Una+95BbUBOW+18AAAAAPAR9Y7tZt7u2k7pTlL5Wx+4XgJCbsRPeecaAACc6vWuLeQZxXGuL1m0Q54NhAmlG/C5khdLpPMfGF4gnQ/gVF0aPeJ6CQAAnFHJh65XEB6UbgAAAAC19vJ/Pimd373nFOl8wGuU7iT9OWuGOIGjXwAAqKtfblrjeglA6Lz9rVzXS4AHmndgM9/6QulO0uvzFmoDbtaOh/8VPD3A9RIAILD+OEW790rRJOl4AHBG/YqjWXRec6R0J2nBPZul8/Ol0wEACLfxOe1cLwEInTWT3G0kXHmkTDo/Uzo9WLjxU38o3QAAILSO/h/145GXiOcD/nNi+UfS+Q2uko4HPEfpBgAAAAJm29mNnGUv/1i751B36XTAe5RuAAAQWofTfy9O6COeD/jP8pZ/lM7vzhMkCBlKN+B3h5aIAwrF8wEAQH1LHX23BynrPMgAwo/SDfjc1opK6fxm0ukA4NZNr5+nDbhYOx4AnJHf+DGLys0fSndAVVZudb0EeKRJt2gcpQAACg2LrnC9BEDimWf03x8MHy6PgI/NfFV/TvddPeURvkDpBgAAAAKm/cX/7HoJAGqJ0g0AAEKr8MoU6fziFdLxAOBMzshoPPrtBUo34HOVR8qk8zOl0wF80+rN2n0aLK+Ddn7A/HpgjeslIOSGDh3qegmAhPp7ULPofB9K6QYAwEOlTfXvyOErWd+aJJ2/zihcUTdlyg+d5N6zf5Q8Y4NpjwYDooLSDQAAQmvMVYNcLwGQeH3eQn3IzfoIIAoo3QAAILSOHOqnDWitHQ//e+21FtL5ffrEvn7s2jeluQDqD6UbAAAPlbxYIp3/wHCOGfy6w38p1QZk8ng54KXZr2yQzs/NkI4PlPJN4j1ILDq/3pRuAAAQWoOq2rteAkKuT/sH1Qni+Ykb+tkFrpcABAqlGwDgRNaIua6XgAiYs2W4dH4/K5fOh/9V/WmydH6ze6Tj6+Tld6+Qzi/qKh0PeI7SnaT1FUel8/M7S8cDgDNrJn3gegmIgH5/aOR6CUDo3PHRM+IE7Ws4gNco3QE15LkRrpcAAElZOkd/K6N7njwCfjd2vXi+djz878XvpkrnF8W5Xl2mf9/W8jvoM3BGlZVbXS8BSaJ0J2n1ZvFfeHmx/7K7r+KANtfM8tO1u3EmasmiHa6X4AQbhiCsqjqO9SCFM2ajrm9GrnT+POl0BIGru77NO+wS5wKoL5TuJN1dpn6XK0c8HwDcqFx7rz7kYn1EogqeHuB6CZFyoGykNuBG7Xj4X8Mi7fvN8ahPQjCLfxrCyLlZ0tzieLf3ERlhu7tP6UZg3D9P/4dvTc90eQYAwDuzpt2qDXh4n3Y+fO8vr3aTzu98fezrLj/AK85RH5U3RDwf8BalO6BWtlriQUqGBxk4k+z0ma6XAEjoH8k08+NmPFs3af9MN7u8UDofwKk6r++oDbi+Rju/Diqmpknnp42Vjgc8R+kGAACoZ0OHqu8Ewi8+bzlFOp9v1oHg489xQO3ZWyzPaNfOXy/UjM9p53oJAJC0Jt1ivyMJjXkP3i+drz0FHEEwv+cn0vlxni6HxxaWnq8PifOmAqcWBR+lGwAAhNb69eIjwxB5t0/RPl5+/QTpeAAeaOB6AQAAAAAAhBV3uoFa2JP+tjwjLd4zRUBIuTpmB9Fy772PuV4CoHHIi011/bcxo6sNZu9v876TXIQDpRuBsbzlH+UZ3e0SeQYAwDt792nft01Pbxrz+sRpzaW5Zmam3b8LAFBPKN1J2lq9y/USAABAHI0urRAntI159fzjLcS5iLrKNrnyjGbyBPhZeXml6yWEBu90AwAAAAAgwp3uJI2Y2Mb1EiJj6GcXuF4CgHpUeGWKPKN4hTwCAOChso3aV0bwFVfvz5uF75i0ei/dn3/+uZ11Fl1e7eJFz+tDcvQRibjwzQ36kK76CAAAgGRVHimTZ2TKE+BnpaWr5BkF+bFvYN5XcUCenZ/u3WtACbfjIUOG2K9//Ws755xz/u7HNm7caLfccott3LixPtYWDPKdI/23ayS8xSe6AFB3s94plc7P+u5V0vkA4Epl5VbXSwiNhEv3zJkzbcWKFfbCCy9Y586dzczs+PHj9uijj9ovfvELa9myZb0vEjDjeCEgbN4s+I0HKUM8yAAQZVE8SurG8Y2l89dJpwPeS7h0b9iwwQYOHGhdu3a1hx56yPr06WMFBQW2adMmu/322+3RRx9VrBPg/U8gZA5u4/xkAMF37uKntQFDh2rn18GYqwa5XgJCbmUrL86hz/Ag4wsJl+7vf//7tnbtWnv88cdtzJgxNnbsWPv+979v//3f/22Zmbz5AZ1Jj2xxlr120Xx5RtcCXiVAtGw7u5HrJSACFu/aqQ1gLxBEUOWi87UB/vucAR7bs7dYntGuXZE840t12vHsvffesyVLltjx48etRYsWtmvXLluwYIH90z/9kzVs2LC+1+hr6jMSOR/xKyfe6SjPaHC1PCJhlfs3ixMyxPOB2FJH3+1BCg8pRt26Tn2d5C6do2/j3fPkEaiFfh22Ocmd/Yp+g9ncjNjXB1W1l2f7TRS/ZtSfhEv3I488YhMmTLDvfOc79tprr1nPnj1tzJgxNnr0aPvDH/5gU6dOPfmuN1Cf5vfUbyh2vTwhcWxigbDqO/dCeca8G+UR8Ludb2rnXzYw5uWqjmO1uWZm9kcPMgB/cXUzIvVEd3GumVmNBxlwIeHSPXbsWLv55pvtN7/5jTVu3NjMzMaPH295eXk2ePBg69Kli3322Wf1vU7Aek734MOc0foIAED4dZ93hz7kHn0E/MvlGcpRtKP/dnkG99LDK+HS/ac//cmuv/7v7wd27tzZNm7caA899FC9LAwAACBpu95yEntBw+c8SPmRBxmAv7h6AnBfhf4V2va07tBKuHTHKtxf+ta3vmWPP/54UgsCAAAIug3p+j0L/PYy35LZ6sd+/SmrZoF0/to410tLV0lzzcwK8tvEvK7fk8N/+3G8n7bUg5S2HmTAhTptpAZETXb/ia6XAACog9930z6Cmy+dHiwXf3eZBymXeJCRmFnTbtUGPLxPO78uxq4Xz9eOD5r7Kg5I5+ent4h5fchzI6S5UULpRmCw0zEQLs8Uvq0PuVEfAX+7evNRbYDfbjc79PK7V8gzinx4RNvRpYtdLwGAz1G6AQBAaFWXVWoD8jto5wfIHR8940FKiQcZiZk4Ufs03LRp06Tz62LJldoPGnKk01Fb6jvsZvHvsnf9n7nybC9RuhEc6keZzHicCfDQvAfvl2cMlyfA78oryqXzM6kHiKDt27U7eefk8Ocq6g4ur5ZnpF0rjziJ0o3A4KgGIFzWr/fggzQAELv33sec5LraxRvRsbLVEg9SMjzIcI/SjcCoqKiQZ7T34VkNbGKBsHL1jSoA1Kepv14tnT9xcl/p/KDh+yIEEaUbAHwga8Rc10uARyqPlEnnZ0qnB8/W6l2ulwAAgbRnb7E8o127InmGHzRwvQAAAAAAAMKKO90A4AOrfuLFu3k3eJCBM5n9ygbp/NwM6XgA31DTb5U4gcfLgaCjdAfUnnT9+bZp1k2eAQCA0oDWx1wvASG3eNdObYAPzyYHkBhKNwD4QIPmf/Ag5T4PMmpv775P5Bnp6U3lGUHBvgHeyqx6w4MUf50R3rDoCtdLcGJdJ+5EAzg9SjcAwIkF335WntHDHpdnBAWvMHir3c6W8ozd8gTUymLx32WXDYx5mV28geCgdAMAEAG/2t1TnlF0lTwCgA9cc81PXS8BHnhjRewPfOpTjxx5hC9QugEATjz+D5e6XkKk3PHRMx6klHiQAcA19etBvBrkD6Ou+LHrJYQGpRsAAKCeTbjpMtdL8FzhlSnyjOIV8gjUwsLZW6Tze+S0lc4HvEbpBhBTdYMxrpcAAEmb/Y72iLYH4lw//JdSaa6ZmWUO1WfAt+6rOCDPyE9vEfM6x6QBiaF0AwCA0Cp4eoDrJQCh8/ibq7UBHJN2ipWtlogTMsTzQelOUuWRMun8TOl0BIH6k+x4n2Kffzz2dQAAEG0frH1LOr+5dDrgPUp3kma/on1sLTdDOh6I68Tyj+QZDdjpGACAOtHf/TTjDqg/7NlbLJ3frl1RzOsHl1dLc83M0q6VR/gCpRsAAAAIGHURM4tfxnaP3yvN5U43wobSnaTs9Jmul4CQ4z0eAPWhYdEVrpcAyPxy0xrXS4iUbYc+lc7PtsbS+YDXKN0BtXbRfHlG14JCeQaACNv5pj7jsoH6jID4fKo+46whsa9TiKBWcNkCfcjxLH1GQGTMeFcb0Pki7XzAY5TuJJWWao9MKMhvI50P/3P1Hg+g9manR+UZ2fKE4Jjf8xN5xvVxrssL0enK0CH100J8QO0HvW/Tbz3LR0dfad5hl7NsVxvMAsmgdCepsnKr6yUAEss/1p+R2V2eAD8rfEm7EaWZ2YaHcuQZ8LeZr2rLwV09Y18fVNVemotTrfqJF9+P3eBBBs6E1+6iYU/62/KMNOsmz/gSpRtATFUdx3qQ8kcPMuBXs9LYKsdLlYvO14cM1UcEReoJLz5WrPEgIxiqzjpPnsH9T3+I4hOAFy96XhvA59NylO4kDXluhOslAEAglcwokWc8mv+kPAPAF4pzvPjUJfbGAS0u8GDTAl4lOOnQX8qk85sNl44HPEfpThLvlURD2fwqeUaXGzrKMwA/GfHTDq6XEClD9w/zIiXm1fOP89+yKKiYmibPSBsrjwiMN1boN4rsEecOaHmV9rUN/uuAsKF0IzBq7HvOsgszt3uQwrM9AIDk1dRMcZLbbOxDTnLNzJ5aM02eURTnNf0P1r4lzY33Is6tx3pIcwHUH0p3kmbW/Fic8J54fnAsnL1FntEjp608AwBc6JuRK8+YJ0+An7ncIf+Oj56RZ5vFfiVm9/i90lR2v4ArYdvMzCVKNwLjz1kzPEjR79idqEbnzHS9BECj8zHXK3AiO50/0win26foX5O6fkLs67d11t/1dfP8AIAwoHQjMNYd0L9XDQBqpaWrpPML8ttI5wN+dN0z4/Qh3frrMwKi/PAh10sAAoXSnSTuQnpnZJtn5RmT5QmJW7BggXR+VlaWdD6AU1VWenGeMACvbDv0qXR+tjWWzq+LgqcHuF4CECiU7iSlzf1MG0AfOqn5UwX6kLzV+gz4Wk3Nq66XEB2HlngQ4r8jfnaOu8n1EoDQye1+p+slAEBclG4EBp+qAuGytaJSntFMnpC4KG7A+dpr2iPD+vSRjq+TKG5cd3+b910vAV6Rf2ga/wPTKD5lqt7QLN5mZmsXzZfmmpl1LfDfh+MKlO6Ayu4/0fUSAABAHAfKRupDbtRHJGLoUP1rYPHOgne5czqA+he2wk/pRnBE9FFUeGvYMO3j5VOm+PCWHDwVxbs08NZTv/+BdH7R4Dg/MHa9NPeLjNiXjxzqp89uHfvy+oqj0tj8ztLxADxA6U6Sq8c9EB23HtMfgwIg/Fzeee3T/kFxsP8+zJpx8w+dZevPrI59XjWio7KN9vUJP74aBCSD0o3AUP8Fbxb/L/kpvY7Is4fJE1AbE6c11wZw0OtJTbp5sDkiTpo17VZ9yMP7Yl5+qSZdGnu6e5w5I3mCCahvlUfKpPMzpdMB71G6ERjF/6H/ZH3aeEoAANS3yrX3agMuPk025SAS8l59Vx+SmRbzcuNJ4u9P8p7Uzq+D8k3ajTBzM6TjAc9RuhEY2em8Bwm9849rd1nGV9RlyIxC9HVHly52lt1u4UBtwOBy7Xz4Xurouz1IWedBBoAwonQDABABEyfqT72YNm2aPAP+lVWzQJ6xNs71vnMvlGfPuzH29QGtj8mzAQQbpRuBUVIyVZ5RkD9engF/O7H8I+n8BldJxwNx3XvvY66X4ASPwQIAXKN0J0l9hlxUDoyH/8jfwTQ77XuYCD91GTKjEH3d3n36c4zT05vGvL6wokia21c6HTi98qpd0vkd4lxXn6Bjxik6QH2hdAOISX/kjJkfj51Z/rH22/fu0ulAfI0urfAgpW3Mq8U5Q8W5Q8TzEzf116vlGRMnx/776rbO2qMm4x3C4HKHfJdH4gHAmVC6k5TdX/+OHOCC+ps2M3+envVy07ek87vbJdL5QVJerr/TDX/4vKX2TzvfzPiDy836AMDP+O8UAmPnuJtcLwERcN0z47QB3fpr5wM+NL+n9tH266XT66am3yoPUmLf6Z7S7YgH2f7i8i47AJwJpRtATPLyaUYBjTiOAYyOI4f6aQNaa8cHzthh4vnq1wUSN+baXvIM9ub/Ck8qQS1sTxNTuhEYM2t+7EHKex5kBEPfF+90vQSEXGmp/k5gQX4beQbgNzzOD7Uofmi66cN39CHt9RFwg783AQBOVFZudb0EeOTwX0q1AZn+u/O6eNdOfUjX2Jf/o0Yb679fbUDP5Z9pBB+lG0BMhVemyDOKV8gjEsbRRkD9O3fx09qAof6rgesaXeh6CfBI+eFDrpfguZKSqdL5BfnjpfPr4vV5C/UhN8e+zBHFwUfpRmBcvOh5fUiOPgL+NumRLa6XAIROvw7bnGXz7inUCp4e4HoJ8MCCezbLM/LlCXCF0h1QZfOr5Bldbugoz8CZ7Ul/Wzo/zbpJ5wNA34xc6fx5p/mxKL57Co8dWiIO8N9dyCHPjXC9BCBQKN0A4APHX/uhPKOh9sl5AIikrRXapymaxbmufuTYjMeOgfpC6UZgqO/4mnHXF+7c1rmHPEO7n3HiuFMSHc8Uiv/+vlE7HjidsrJj0vk5Pnzm+L6KA9L5+ektpPMBr1G6AcRUnOPFxkRDPMgIBs5FR5ipz1Dm/GQAgJ9RugEAAEJk6P5h6gTxfPjdylbq99gzxPMBb4WidFc3GON6CZ4rzNzuQQpbeSN6Tryj3UCwwdXS8YGifjzRjEcU/eLeex9zll228RNn2a5MKZshna+u9AAQNqEo3QAAAPjC6/kvSecPi3eYMACJ7P4TXS8BSQpF6e59W6Z0/hrpdOD01LuTxtuZ9POW+m23QvEXEOpM/3iiGY8owpldb7leATySM5IdvgGcHt/zIjA4GsNb83vqH8m8Xp4AIOoq928WJ2SI5yfuQNlIbcCN2vEAEDahKN2rfrJVnHCDeD7gP7dP0b7bbGZ2/QR5RGD0ffFO10sAZKb+erV0/sTJfeP+WGWl+nsE/5k17VZtwMP7tPMDpvJImXR+vOc5eeQYYVY2v0qe0eUG/fe6XwpF6V7+cfz/2NaH7tLpgD9tO7uR6yUAQGD9vttMeUa845s5os1b5ZsqpfNzM2JfD1spAcIsFKUb0ZDapMj1EiJlzpbh8ox+Vi7PAODepSP1O9XjKy53iwcA/L1QlO6j/2eXOOES8XzUBsekAUAwbfrwHW1Ae+141M6+iobyjPZx/l3zQQMAPwtF6QaAoPt8qj7jrCH6DCCWxbt2agO6asfXxfqKo/KM/M7yiIQsnL1FntEjp608A/6WNvczbUCWdjyiidINAD6Q9a1J8ox1NlSeAf+a9U6pPCPru1fFvL6u0YXy7Hh2jrvJSe7qzdr3fM3MLK9DzMsfna8tv+kW+9+zS3v36U/cSE9vKs8A8IWwPeEaitLdp/2D6gTxfAB+oS6/FF8gGmalNXe9hEjhLru39uwtls5v1459fBAuoSjdAOpfvz9Ec/fyg9t4L9Ar6m/azPjG7ese/4dLXS8hUqrLPLjTnR/7TncU/cvo810vAQDionQjMKb0OiLPGCZPgN9xVBpCa/Gz+ozLBsa+vustfTac+3PWDA9StMfEBknYHr8FwozSnST1GYmcjwhnxq73IEMfkSj1UWnxjkkbc9UgaS7g8tzoKCqv0B+JmEkhAiKBvhF8lO4k6T9l5D+oXxr27sf6kNb6iERl95/oegmR0q/DNtdLACQmLXhXnpHf+SJ5BgD3Gp2j/xAPX6FvBB+lGwC+pm9GrnT+POn0YHljRZxHketRD76PgCM5PVNdL8Fz6w5o78YBQFBRuhEcK3+rz7hmgj4jQTxS5K0DZSO1ATdqxwPwh0N/KZNnNIvzNoz6eLh4R8MhOjIOsf8JkIhQlO6n1kyTzi9qLx0PAABEVrZaIk7IEM9HbeyofkeewYcNAOoqFKW73ULxI4qD9ZuhAH7z/tQ0eUbLsfIIAD6g3sSNDdyw6UN96TZuwgCoo1CUbpfUx1hxhBWAsBp1xY9dLyFSSj50l63exI0N3ACoZczQb0Zp/F0WWpRuADEdXbrY9RIi5fYp+nfrr/fZlgUHl1fLM9KulUcERuqR2R6kxN65zmXh37O3WDq/Xbsi6fwg+WCt/jz25nGuL961U55tXfURAMKJ0g0gpsWL9aW7fXue1fvS7gljPEi5zoMM4O+pz6zmvOpTPf7mam2AD8vnukYXul4CAMQVitLd98U7XS8BHpjTe4c8I0+eAMQxdr0HGfoIIJat1btcLwEInfKnfyPP6DAszhb5ABISitINAEGXVbNAnrFWnpCYPelvyzPSrJs8A2c24qcdnGWrz4OPdxb87Hc2SHPNzB6QJyRmwT2b5RlsmocoYg+p4KN0A4APzJp2qz7k4X36DCAG9ZnV8c6rNovmhn0j2zwrnT9ZOr2OdunfJweiJrUJe1bUF0o3gJiuueanrpcAoB41uc3d3eYoKnh6gLPs1ZsrtQF5/vu9pD6Wzoy77H6hfkqKJ6SgEIrSXXhlinR+8QrpeMCX9u77RJ6Rnt5UngEAThxa4kFIoQcZAIBkhaJ0F+cMFScMEc9PnPrdDjPe7/CLwszt4gR2/YUbaxfNl2d0LaCUwI2Zr+o3j7urpzwiIesrjsoz8jvHvs4ZygD8LBSlGwCCbtuMZfKMlvIEAHBj1uZSeQZH0yFq9Dd+zKJy8ycUpfvzllOk80PxiwQkaOHsLfKMHjlt5RlBwa83EC45I3nKwksud8h3ZcWCn8szoverCmjQJwHEdOnIA66XECn8egMaB5dXS+enXSsdj1pS75BvFn+X/PJy8cZ1AAKP0o2ELVm0w/US4IFf/fdCecag9h4ckxUQN//t266XgLBjY6/IKPnQ9Qq8V16lf4c+3l3f1KO/lWfHMuzdj/UhrfURQBRQugHADxZrz9U1M7PLBuoz4FtR3NgrqlKPzBYn+O8dzC6P5LpeAgDERelOkvxTxtaxL09ePl2ba2bDWt8jz4B/vT5Pf6fbbtZHBMXINvrSPVmekJjs/hNdLwEIJVdnshcs9eAxax+eEQ5vqU++4NQLKFC6AZ9THw/H0XD+cFe16xUACA35qwSxS0l5Rbk4lx3E/YK9EoDEhKJ0Z31rknT+OlOfAw74z+7xe+UZzeUJwbHmVf37tm3yuUMEeKXyiH5jr8w41zdcoP3b1ZcPckdxz4KVHrxLfs0EfQYQg/qmk5m3N55CUbqdUv+FF+cvu00t+2hzAXhqa7X+fVtEW1SPsNqT/rZ0fpp1k86vi/JN2se8czOk4+tka4X+0fZm8gTAX8JWfF2idAM+52rfAADh4vLOaxSpi6+ZP8svoiGKH2YByQhF6T647THXS4iU++dtlc5f0zNdOh+1M9ODbxiz2zWWZwTFD1L170Ii2lyWwNnvbJDmPiCdDgBAckJRul2a01t7ZnWedDoAv+C4G4RZwdMDXC8BIVdWdkyekZMvjwB8hdOS6k8oSveEmy5zvQQgdNhN21tbN82UZzS73F/v9JbNr5JndLmhozwDteBoN214i70pANSXsBX+UJTuw38p1QZksnv5143Paed6CUDoHMv8geslADKVbbRPcrDBFRAd2f0nul4CkLBQlO5BVe1dLwEAksKmSwgz9SZubOAGRIf6KSmekIJCKEr3nC3DpfP7GRscAVExa9qt2oCH92nnAz4UxSOsgDBbu2i+dH7XAl4ZQbiEonT367DN9RK8pz4f3CzuGeGIhuoy/Z1Xy++gz0jQthnLpPNbxrleXu7Brzcijd9j0dH48yed5I74qf/+TgcQTJta9nG9hHoVitINoP6VV+if8Mi0HHlGohbO3iKd3yOnbczrW1dcLc0FstP1m/XhKy4/5Cjb+ImzbERDapMiZ9mFmdvFCf773gTBR+kGEBO70ALhQhHzlssPOTasf0WckCGej9pQH1trFv/oWoovkBhKNwLjzrnr5Bl5WdfJM4KCxwQBoO5KS1fJMwry28S8PuS5EfJsAEDthaJ0TymbIZ0/TDo9eJa3/KN0fne7RDoftXPoL9rdhs3Mmmn3QAQAAACcC0Xpfj3/Jen8YXazdD7gR+VV+sfLuZf+lZ3jbnK9BAAAAAiEonQfKBupDbhROx7wo/LDh5xlZ42Y6yzblZWtlniQkuFBRjBE8fcYouO+igPS+fnpLaTzASBsQlG6AdS/gqcHuF6CE4/94xPihL4xr+7ZWyzONWvXzt1uswC886/L52oDumdo5wNAyISidM+adqs24OF92vkAIu+NFQPlGT3YDBZASOWMLHSWzW7x3prS64h0Pns5QSEUpRsAAPhX48+fdL0EAECCNrXs43oJoRGK0j3vwful89lgGQAQdJX7N3uQkhHzKmeER0PJiyXyjAeGF8S8XnlEf+JGZpzrHNEG4ExCUbrXr1/veglA+BzyYmMvd48DAkBYjbrix66XAKAeDXv3Y21Aa+14hKR033vvY66XAKAerfrJVnHCDeL5gP9UVqr/XPnT2kXzpfO7FvDhIQDg9EJRuvfu0z62lp7eVDofAAAg6FyeelG+qVKekZshjwAQUqEo3Qtnb5HO75HTVjofwKmWfxz7WK360l06HfAn3juF2tZNM+UZzS7nyQIAwROK0g2g/m2t0N81aCZPAOAH+k3cMsTz4Xfl5fr/ZgHw0Mrf6jOumaDP+F+UbgAxlZUdk2fk5MsjAoONj6B2X8UBeUZ+eouY16P6PnnUNOkWe2dxAIg6SjcSNvb/xv6mCiLqT/o8/JQP8IvxOe1cL8FzM2u8+GDnPQ8yEpPdf6KTXD5oQJhN6XVEOn+YdDrgPUo3ACByul/ixYZPGz3IqL1G5+jft41n57ibnGUDAOBaKEr3n7NmiBO0mzoBOFXZp3+TzvfjRmoHl1fLM9KulUcExrAP/l2eMaW9PCIwVrZaIk7IEM9PXBQ3rqs8UibPyJQnAED9C0XpBhAuRbs6ixNqxPMTtyf9bXlGmnWTZwTFlG7aRyOBKHJ5bFd2ursnOeCtycunS+cPa32PdD6iKRSle10n7kQDCLa1i+bLM7oWcNTOlz6fqs84a4g+IxELFiyQZ2RlZckz4F/sIA4vbGrZx/USgISFonTbzje18y8bqJ1fB3N675Bn5MkTAAAIp24pL3uQUuRBRu0deuc210twQv/v2l//ngEkLhylG0C9yxnJXVEACCL1HhHsDwEAiaF0I2GFV6ZI5xevkI5HLbncEOfzllOkuaf7i+/Nvq9Ks7Ol01FbWd+aJM9YZ0PlGUGxZ2+xdH67dvHvBJbNr5Jmd7mhY8zr6q/ZLP7Xrd4jgv0hTuXy9zcgpT621iwyR9eGo3Tvesv1CgCExLZDn0rnZ1tj6XzAj95YoX1Nq0eOdHydqL9ms/hft3qPCD/uD1G28RNn2S6P4wMQDOEo3QDq3exXNsgz4u1CC6iNuWqQ6yVEyqgrfux6CZGS3X+i6yV4rnL/Zg9SMjzIABBGlG4AACKA4gu1PuPdPXJQWbnVWTYAnEkoSvcHa7WPlzeXTgf8iTNPEWbXP//v+pCxnPWKaIniju0AUBuhKN27x++Vzqd0AwizrBFzXS/Bc30zcuUZ8+QJABBB6s29IrKxF7wVitINIFza7Wwpnb9bOh1B8EyhdndnMzO7UR8BAAD8LxSlm92GgXDZdnYj10tAyI25tpc8Y5o8AfAXl8ekAYCfhaJ0r684Kp2f31k6HsA3pI6+W5ywTjw/WMbntHO9BAAhENWjsxYsWCCdn5WVJZ0PxDOn9w55Rp48wR9CUboBtWHvfqwPaa2PSERp6Sp5RkF+m9g/MHa9NnisdjwAwFtDnhvhegkAEBelO0l3ztXeMcvLuk46H/CjKWUzpPOHnebHovjkTPe3LtUG5NVo5wPwhYxDvBoEALFQugHEFNUzTxtPKtEG5D2pnV8HFVPTpPPTxkrHAwAA+BqlGwAAAEk7uLxanpF2rTwCAOodpRuA7zyXcYF0/ukeL4+iG8c3ls5n2zp/oBBFR9n8Kun8Ljd0lM4HgLBt4kbpBmpj5W/1GddM0GckIKqb0gxofcz1Ejw35qpBrpfguX7973O9BM/tSdefTZ5m3WJeVxd+yv6pzhol3ozScmJedfl7DAD8jNINwHemD4izqzlQT95PW+pBSlsPMoJBXcYoYqeavHy6dP6w1vdI59fFfRUH5Bn56S1iXh91xY/l2QCCjdKNhBXnDBUnDBHPR224/AbGpfKqXdL5HaTTgfjWLpovz+haUCjPCIooFjGXv8dWtloizzbL8CADiEH9xKXPnrYMI0o3Etb6f16Tzt8tnY7acvkNzL6KhtLU9u2l4wF8g7qMUfb9IbVJkeslAIAvUboBn1NvJBFvE4k9e4uluWZm7dr57xu08sOHXC8BHtj04Tv6ED7c8YXCzO3ihNjvNwMA8KVQlO7Vmyu1AXk8EPp1q6b/qzZgwj7tfOA0CkY3d70EeODxf7jU9RI8l91/ouslADJR/aAYQDCEonTPSuObZC8dXbrY9RIAma0V2g/xmkmno9YWP6vPuGygPgMAAPheKEp3dZn4Tnc+d7oRPW+s0BeGHnGeylw4e4s4l12lASBMXP43C9FQXlEuz8jkdZXQCkXpVv8h4A8A4K0/Z80QJ/QVzw+W65//d23AWP8dLzSyjf5O92R5AuAv+vfnzeK9Qx/F3eIBBEcoSndOz1TXS/DcnXPXyTPysq6TZwB+06RbgesleO79qWnS+S3HSsfXiXwvEDP2A/ka3icHAERZKEr36ldnS+f/y/DofRMOuLSu0YWulxAp22f8Xjq/pXQ6AEQTN2CA4AhF6e7ySK7rJQChc+uxHu7Cd73lLLrySJl0fqZ0et0smH1EOj8nzhs6Lh9FBQAA8EooSjcA1JfyTdrHjnMzpOPr5NKRB1wvAQidg8ur5Rlp18ojEjKll/YDPDOzYXGuR/HXG0BwhKN0H1oiDigUzwfwdR+s1d7p9uMhg6lN3J3/mt7gJ86yoyht7mfagKzYl7tULdfmmplZRw8yAAAIlnCUbgAAAB/Zk/62PCPNuskzAADJo3QDwNeUl3uwq7XPzBQ/Up/drnHM6y4fRQXU1i6aL8/oWsCTeF/iQw4AfhaK0j3z1V3S+Xf1lI4HAKebijWeVKKNzXtSOx+10in1t/KMTTZcngEAQNCEonTnjOSTXgDB5vKu74DWx+TZsUxePl2eMaz1PfIMIBbOJgcAfCkUpXvJmlLp/Mx2nNMNREV2+kwnuRRQQEP9gRavMAAAziQUpRtA/Yvq+3Glpauk8wvy20jn10V5lfYVnQ7S6QAAQOHOuevkGXlZ18kz/IDSDQBfU1m51fUSImNTyz6ulwAAAHwobIWf0g3Ad3aP3yudf7pzuoc8N0KaDaifIvHjEySIBl6TAYDYKN0AgMgp+dD1CgAAQFRQupGwiooK6fz27dtL58P/th36VDo/2xpL56OWVuqPsLJrJsS8vOvdX8qjM+Mc0aY+v5mzmxFFnIsOwM9CUbrLyytdLwFAPVpfcVQ6P7+zdDwCoMudqa6XAAAAIiIUpdvVET8ANO6qdr0CQOeunGjs1AoAAL4QitINoP5F9VG9i/aJN6FKf087HwAAMfWmeWyYh7AJRemO4rm6ADQancOTM1Ew81Xt2eRmZnf1lEcERtn8Kun8Ljd0lM4HcCqOfAQSE4rSzbm6QLiUzCiRzn80/8m4P5ZxqJE0Gyh/+jfS+R2GDZfOBwAAiQlF6eZcXQAAgOjK7j/R9RIAIK5QlO5uKS+LE4rE8wH/4RsYAEBQqF9hMOM1BgB1F4rSvWdvsXR+u3aUbkQP38BER/nhQ66XAMiw4ZN3eM8X0JjTe4d0fp50OsxCUrrfWDFQOr9HjnQ8ADiVOrS56yUg5Aozt4sT+A813Dm4vFo6P+1a6XgAHghF6QYAAPCTSD4ttPK3+oxrJsS8rP9gx4wPdwDUFaUbqAX1Yz1mPNqD6OHPFYAw2JP+tnR+mnWTzgegR+kGEJPLuwYjftrBg2x/4V1Ib+WMLHS9BCB0+CANXuD9ZgQRpRuA/3Q+5noFCLkla0rlGZntCmJe7zD5Rnk2AO+sXTRfOr9rAR8SAkFH6QbgP4eWiAP4BgYOqd97jfPOq0tTeh2Rzh8mnV43vGMMAPgSpRsAAACoo9QmHC0L4PQo3QCAyNm17iF9iPY0S/hcp1T9Tt6bbLg8AwCQvAauFwAAAAAAQFhxpxsJ2/NeI20Ar6gBEPvhiO86y2bnXaD+8WQBAD+jdCNh5W994HoJAJCUPTOf1YeMHavPAHzkzrnr5Bl5WdfJMxKl3zSPuxFA0PF4OQAAAAAAItzpBgBEzqgrfux6CZEyefl06fxhre+RzgcAIBmUbiTszgM3ixM+Fc+H322tqJTObyadDgCIkiieQw8gMaEo3dyxQJip35Hz4/txxzJ/4HoJAOrRppZ9XC8BAABnQlG64a3ig0Ol8x+VTkcQzH5lg3R+boZ0fODIz6zmvGoAIcbrEwDOhNKNhPUZ/xPXSwAAAACAQKB0A/Cd7PSZrpcAAADiWSk+F/2aCXF/KIqv3SH4KN0AfKe0dJV0fkF+G+l8AP7B++QAANco3QB8p7Jyq7Psg8urpfPTrpWOr5Ob+x1zvQQAQIDM6b1DOj9POr1ucoYv8SCl0IMMuEDpBuA7O8fd5Cx7T/rb0vlp1k06HwAAAP5C6QbgOzNr1McAvhf3R9Yumi9N7lrgv0+xMw41cr0EzxU/f588465rX5NnAAC8sbWiUp7RTJ4AVyjdAPA1qU2KXC8BHuDfc4Q43PAJAAAzSjdQK+qdMs38t1vmlF5H5BnD5AmJK8zcLk7IEc8HAACAn1C6AfjOGysGSuf3OE3vLX/6N9LsDsOGS+cDQBSxSz3CjGPSgi8UpTuKuw0DYXbrsR6ulxApG/7/htL5XaTT6+Z32XfIMwptkzwDAAD4XyhKN4D6N+zdj/UhrfURAAAoRfEuO3degcSEonRzxA8QLvyZBgAEBpv1ATiDUJRuAALqbyLMfPmNxIoFP5fO7yCdDvgUpQQAEGGUbgC+E8WzsgEAABBOoSjdfIMOAHW36bLXpfO7WEfpfPjfnN47pPPzfJoND0X06SwAwRCK0p3df6LrJQCoR/yZ9tbk5dOl84e1vifmdfVGPGZsxgNAjw92AJxJKEp32fwq6fwuN3CXBogKVwUUAAAA4RSK0l2YuV2ckCOeD/iP+pN7s/if3vNBGgAAAMIiFKUbAAKP9xEhFsWzhKOKf9cA4C+UbgC+w9MrQLio39/n3X0AgJ9RugEA8JCzAsrTFAAAOEHpBgAfcPkOPRBmrj7k2LRrrzQX/sGTHNHQpFuB6yUgwCjdAHxnSq8j0vnDpNMBAEDYFP9HiTxj2niKfViFonR3StU+MrfJhkvnAwAAAADCKRSlGwDqi6tHQtWPJ5rxiCLgJV4ZAcLl0Du3uV4CAiwUpZv3pgAAAACo9BnPySeou1CUbgD1z+Wd12HvfqwNbq0dD/gRd14RZvz+BuBnoSjd6r9o+UsW8Jj6aKPTHGvE3ycA6oPTV0Yc/h0KhFW3lJc9SCnyIAMuhKJ0AwAAAIDKnr3F8ox27SjdYRWK0s35iEC4cLcZqH9s1ocw4/c31FJ3UYhRd6Eo3QDCxeUHaXyIByDo+OASqH+be/xVntFVngBXnJXumpoaMzP79NNPk5514n+0fwhOt0ZX2erc02UfPXrUSa7LrznM2XzN/sh2+eu9atseJ7lh/ppPlx3F32P8mfY2+44/rJDm9ujYK+b1MP8ec5nN1+yP7H5vVUpzzcw+/f/89esd5t9jp8uuy4wvu208KTVn+hkiFRUV1qpVKxfRAAAAAADUi3379llaWlrcH3dWuk+cOGH79++3c88911JSUlwsAQAAAACAOqmpqbHDhw9bixYtrEGDBnF/nrPSDQAAAABA2MWv4wAAAAAAICmUbgAAAAAARCjdAABAgjfYAACgdAMAkLR169bZwIED7aKLLrKzzz7b2rZta0OHDrVdu3a5Xpondu/ebSkpKfb888+bmVl1dbUNGjTI3njjDbcLAwDAByjdAAAkobi42LKzs+2DDz6wxx57zObPn2+jRo2yZcuW2eWXX25vvfWW6yXKNW/e3NasWWPXXXedmZlt3LjRZsyYYSdOnHC8MgAA3GP3cgAA6mjVqlV21VVX2YgRI+xXv/rVKT/24YcfWmZmpl144YW2fv16Nwt0ZNmyZda9e3dbunSpXX311a6XAwCAU9zpBgCgjp544glr3LixjR8//u9+7IILLrCnnnrKbrzxRjt69KgdP37cfvvb31qnTp3s7LPPtosuushGjRplx44dO/nP3HLLLda7d2979tlnrW3btnb22Wdb165dbceOHfbaa69Zp06d7JxzzrErrrjCNm7cmPQ/d/XVV/9dKV62bJmlpKTYsmXLzMzs+eeft7POOsvefPNNy87OttTUVPve975nTz755Ml/5uuPl39ZuM3MunfvbldffbUVFxdbSkqK7dix45SsF154wRo2bGj79u2r478BAAD87yzXCwAAIIhqamps4cKF1rdvXzvnnHNi/pybbrrp5P++9dZbbfr06TZq1Cjr1q2blZWV2cMPP2xlZWW2YMECS0lJMTOz1atX2/79++2pp56yv/3tbzZ8+HDLzc21lJQUe+SRR6xRo0Y2bNgwGzBggG3ZsuXk/Lr+c7Vx4sQJu+mmm6yoqMh+8YtfWElJid1zzz3WqVMn69Wr1yk/9wc/+IEVFxdbYWGhFRcX29VXX20tWrSwu+++22bMmGHjxo07+XN///vfW48ePaxVq1YJrQcAgCChdAMAUAcHDx60Y8eOWZs2bc74c7du3WolJSU2YcIEGzVqlJmZ9ezZ01q0aGEDBw60+fPnW25urpmZHT582F555RX7/ve/b2Zmy5cvt9/97ne2ePFiy8nJMTOzd955x+6++26rrq62xo0bJ/XP1UZNTY2NHj3aCgoKzMysa9euNnv2bHvttdf+rnSfd9551rFjRzMz69ix48n/nZeXZy+88II98sgjlpKSYhUVFbZkyRJ74YUXar0OAACCiMfLAQCog7PO+uJz6+PHj5/x5y5fvtzMzPr163fK9Z/97GfWsGHDk49ym5k1adLkZHE2M7vwwgvNzOyKK644ee073/mOmX2xS3iy/1xtZWdnn/zf3/72t+2CCy6wo0eP1vqfLygosN27d5/c0Xz69Ol27rnnWl5eXsJrAQAgSCjdAADUQZMmTezcc8+1PXv2xP05R48etUOHDtnHH39sZmbNmjU75cfPOussa9q06Skl+Lzzzos5q1GjRqddT13/udr65iP0DRo0SGh38pycHGvTpo1Nnz7dzL4o3T/72c8sNTW1XtYHAIBfUboBAKijXr162dKlS0/ZDO3rpk6dak2bNj35/6uqqk758c8++8wOHjx4ys/xUkpKyt/dqT9y5Igs65ZbbrE5c+bY+vXrbfv27Zafny/JAgDATyjdAADU0V133WUfffSRPfjgg3/3Y1VVVfbkk09ax44dTz5C/dJLL53yc15++WU7fvy4/ehHP/Jkvd903nnnWUVFxSnXVq5cmfTchg0bxrx+yy23WHV1td11113WoUOHUx59BwAgrNhIDQCAOsrKyrJx48bZgw8+aOXl5TZ48GBr2rSpbd682Z544gn729/+Zv/1X/9lHTp0sMGDB9vo0aPtr3/9q1155ZW2ceNGGzt2rHXv3t169+7tZP19+vSxefPmWVFRkfXt29feeOONk49/J+PLTdpef/11a9KkiV122WVmZnbRRRdZz549beHChfb4448nnQMAQBBwpxsAgCQ88MADVlpaamZm//Zv/2a5ubk2adIk69Onj23cuNE6dOhgZmYlJSU2ZswYmzlzpuXm5lpxcbHdcccdVlpaag0auPnP8c9//nO799577cUXX7Tc3FxbvXq1/ed//mfScy+55BLr16+fTZ482QYMGHDKj1133XXWsGFDGzhwYNI5AAAEQUpNTU2N60UAAIBouPbaay01NdXmzJnjeikAAHiCx8sBAIDcuHHjbPv27bZw4cJ6eW8cAICgoHQDAAC5efPm2bvvvmtPPPGEdenSxfVyAADwDI+XAwAAAAAgwkZqAAAAAACIULoBAAAAABChdAMAAAAAIELpBgAAAABAhNINAAAAAIAIpRsAAAAAABFKNwAAAAAAIpRuAAAAAABEKN0AAAAAAIj8P/ssiRv1a29yAAAAAElFTkSuQmCC",
      "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8QAAAHbCAYAAAANs1OQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdM0lEQVR4nO3de6ydZZk34Hvv1tBSu0MpiHbaykETW8PnhzGhpaItTgMUguIYQiEVBAM2mOBXIDaftSIIqCUMGSkaTRVKEWPQTgjSdgzQCuKYSNo6KgElgsWGk9JQECaG7vnDaccOa5d1etd7uK/rr3YfnvdZa7+n33M/77OGRkdHRwMAAACSGS67AwAAAFAGgRgAAICUBGIAAABSEogBAABISSAGAAAgJYEYAACAlARiAAAAUhKIAQAASGl8WRves2dP7Ny5MyZPnhxDQ0NldQMAAICGGR0djd27d8e0adNieHjsOnBpgXjnzp0xY8aMsjYPAABAw+3YsSOmT58+5vdLC8STJ0+OiL91cGRkpOt2Hn/8pn51qaVjjvl0pba7ect7Ct3u/A9ut92CtzvWNrPtU47d/mr1Pmd6rWVu1/tc723a7v/Idl5+x09+Wdg2f/eB/9Py69n2Kdut9zYPtN2yjttOvPjiizFjxox9uXMspQXivdOkR0ZGegrEf35hTb+61NJxI/+/UtudNKnYx77H+lvYbvHbzLZPOXb7q9X7nOm1lrld73O9t2m7/yPbeXl40psL22bV/ra2W//tVu21lnXcduONHs8tLRAD0EyTDr697C4AALRFIAaAHhgAoN/sUwCDIxADAACNYVCJTgjE0IITKQDQb9nuL7K9Xuqp2KewAQAAoKJUiKEijKICwOA8OOO+Alv/vwW23Tn3GDA2gZhKcwKH+tm4cWOh7c+ZM6fQ9gGAPARiAKghA4YAlKVJ1yCBGAAAoEdlhMQmBdOyCMRAKi4czeVvS7/ZpwCaTyCuGRdn+s0+NRjeZ6BdnsMHGJzaB+KybjLd3AKQkesfAE1S+0AM1JObappCNQ96U9b14Mk/rC6s7Xe8Y1lhbfPGnJfphEBMW4QXqB/HLU1hXwagKAIxAADQdyq11IFADC0UeQJ38gYABkkwhbEJxAD01fLj/6nsLgBdMDW9eIIpVM9w2R0AAACAMqgQAwAA9Mgjd/UkEENFmEYFAEAn3D/2TiAGgB54ZhoA6ksgBoAaUhVoLoMs0BvHEJ2wqBYAAAApqRAD0FfPb9lVaPvTTy20eQD6RKWWOqh9IDZlDAAAYHCalMFqH4hptiYdbECxnC8AgE4JxDXjhg96U9Yx5NiF+nHcDkam99kUYqgegRgAAKBHBjzqySrTAAAApKRCDEBfPfnOXxba/vQ4sdD2AYA8BOIuZXreBQAYHB9dRr+ZygtjE4ghOYM7QCeKPGc4XwD9YFCJTgjEAEClGbgDoCgCMbRgahEAQG9UaqkDq0wDAACQkgoxbTFdDerHcQsAzWZWY+8EYqgIJzQAABgsgRgoheolAABlE4gBoAcWjaHfnnznLwttf3qc2PLrBiqL53wB1SMQAwAA9KjIAQ+DHcURiAEAoMFUpum3Jq19IxADQA016WakqrzHUE9lPXZAPdU+ELtYAQAA0I3aB2KATpQ1iGbwDoBsVGqpA4EYgEYw6AD15NgFyjRcdgcAAACgDCrEAACkM+ng28vuAlABAnHNZJtWVNbr9TlyAADQfAIxVITPCARoLdtgMACDIxBDcm40AQDqSUGldwIxAAAMgI8hguoRiIFSqExDPTl2oX4E8cEo8n32HhdHIO6SGwKg6pynAAAOzOcQAwAAkJIKMQAAAG1r0mJeAjFtMfUSoDXP5gFAfQnEAAAY/AZS8gwxAAAAKakQAwC00KRn5ABoTSAGGIBMN9b//m8bCm1/3oWXFNp+XWTapwCgKAIxVISFeQAAYLA8QwwAAEBKKsTQQpHVWpXacplmCgDAXgIxlSa8AO1yvgAgG4/c9a72gdgNEAAAAN2ofSAGAGgSK7UDDI5FtQAAAEhJhRgAgNJ4/A0ok0AMAEA6GzduLKztOXPmFNY2VEGTFvMSiGvGKCr9Zp8CACArgRgAqDQDdwAURSAGAIABsIL4YHif6YRVpgEAAEhJhRgAAOg7lVrqQCAGANpW5PO8VXuWt0mrqALQmkAMMABurAEAqkcg7pIVLwEAgL2KnCJuenhxBGKAhjJw12xmHQBA7wRiKs0NH1B1Fo2hKQyiNZfzFIxNIAZK4cYLAICyCcRAKmYdABARsfz4fyq7C0AFDJfdAQAAACiDCjFtMb0VaNfcc75aynZV/wGATqkQAwAAkJIKMbTgc+QAAKD5BGKoCB+JAAAAgyUQAwAA1JCCSu9qH4gtogIAAEA3ah+IAQCKoPIC0HwCMSRnlgUAAJ1o0oChQFwzwgsA0CTubYAyCcQAQKVlC0xzz/lq2V2AWnMM0QmBGAAA6DvBlDoQiAEGoEnP2gAANMVw2R0AAACAMqgQA6mo1NJvpgQCQH0JxABQQwZ3AKB3AjEA0LYiV3yu2mrPADSfQAwAAAPgEQuoHotqAQAAkJIKcZeKnDIWYdrYXp6RA9rlfAG9cW8DvTEDoJ5UiAEAAEhJhRigoVR7AAAOTCCG5EwzBQAgK4EYAAAazLOtMLbaB2LVLQAAALpR+0CcjQEAmsLzrUC7XPsAKIpADABAaQx4AGXysUsAAACkpEIMLVh8AgCAqnPP2juBmLZ43hMAms3UZSAjgRgqwggfAAB10KT7VoEYAGhbkVXEqlUQm3TDB0BrAnGXTCsCOuHGGgCgegRiAPpq64anC23/hA/PLrR9IIfnt+wqrO3ppxbWNNBnAjEAANB3ZQ2QGpilEz6HGAAAgJQEYgAAAFISiAEAAEjJM8QADWU1fACAAxOIqTQfVQMAABRFIAZSMcgCVJ0VcgEGRyCG5MoKiKbzQm8M7kD9GOyA6hGIaYvwAlRdWQHRDS4AEcVeD1wLimOVaQAAAFJSIQYAaEH1H6D5ah+Isz1Dle31AgAAFKX2gRgAgPoy2A+USSAGACpNYAKgKBbVAgAAICUVYgAAaDALxDWXv23vBGIAoG2mLwPQJKZMAwAAkJIKMQAAqv9ASgJxl1w0gE54xgcAoHpMmQYAACAlFWJoochqnkoeg2ImCwBQhCbNfBOIgVKUFdaadAIHAKA3AjFUhKAGADTJJcc9WvAWTiq4fTIQiAEAgMYQxOmERbUAAABISYUYkjNVGwCArARi2mK1WgAAoGkEYipN9RIA1wIAiiIQA9AIQhMA0CmBGACgBSvVAjSfQAwANaQiDgC9E4hrxg0QAABAfwjEAPRVtmmm2V5vkQOzBmUZpCff+cvC2p4eJxbWNtVV7PWgWteCJql9IFYxBQDonXsqIKPaB2IAAOpLEAfKJBADAMAAZHvEAupguOwOAAAAQBlUiLtkeg8AUARVRPrNPgVjE4gBAABqyGBH70yZBgAAICWBGAAAgJRMmQZKke05/DKmNGV7jwEAOqVCDAAAQEoqxLRFpYmmsPgEAAB7qRADAACQkgoxAFBpZnYAVEuTzssqxAAAAKSkQgwtFDvqVa1KRJNG+ACgXf/+bxsKa3vehZcU1jbQXyrEAAAApCQQAwAAkJIp01ARpi4DUCbXISAjgbhmXKwAAAD6QyCm0gwAAO1yvgAAOiUQA0ANGQAAgN4JxAAAQN898s//Umj7sy5eWmj75CAQAwAAjSGI0wmBuEumqgGQUbHXP9c+AAar9oFYMAUAAKAbtQ/EAABFMO1yMBQ3gDINl90BAAAAKINADAAAQEqmTNMW05kAAICmUSEGAAAgJRViAPrKQkQAZFTk9c+1rzgCMVAK0/AByMaAIf1mn+qdQAzQUAYdAAAOTCAGAIAGU0WEsVlUCwAAgJQEYgAAAFIyZRoAemAqIv1mnwIYHIEYAIB05p7z1bK7AFSAKdMAAACkpEIMMACmQEL3HD8AFEUgBlJxYw0AwF6mTAMAAJCSQAwAAEBKpkxTaaa3Fs97DABAJ5p0/ygQQwtFHuQCIgAAVINADABAoyo+AO0SiKEi3IgAAMBgWVQLAACAlFSIa0YVEaA150cAoFMqxAAAAKRU+wqxigAADI5V+Ok393JAmWofiMvi5A1AmVyHAKB3pkwDAACQkkAMAABASgIxAAAAKXmGmLZ4Vo1+s08BAFA2gRgAoIWfbLyg0PZnFdo6AO0QiAEaShUeYGxbNzxdWNsnfHh2YW3XSVmDSgaz6IRniAEAAEhJIAYAACAlgRgAAICUBGIAAABS6mpRrT/84Q9jfm94eDje/OY3xyGHHNJtnwCoMYuZAAB10VUgPvLII2NoaOiAP3PooYfGpZdeGitWrOiqYwAA0CQGDOk3+1TvugrEt956a1x00UUxf/78WLx4cRxxxBHx7LPPxg9+8IO4++674/Of/3y89NJLcc0118TUqVNj6VIfzQEAADRXkeE0QzAtS1eB+I477oizzz47vvOd7+z39SVLlsTSpUvj4YcfjrvuuisOPfTQ+PrXvy4QAwC0ScUHYHC6WlRr8+bNcc4557T83kc/+tG49957IyLihBNOiMcff7z73gEAAEBBuqoQT506NbZv3x4LFy583fe2b98eIyMjERHx0ksvxaRJk3rrIakZJQcAAIrSVSA+99xzY+XKlfGmN70pPvaxj8Vb3vKWeOaZZ2L9+vVx5ZVXxqc+9al44YUX4sYbb4w5c+b0u88AXTPIAkA2rn0wtq4C8Ze+9KV49tlnY9myZbFs2bJ9Xx8eHo4LLrggrr322rjzzjtj69atcd999/WtswAAANAvXQXi8ePHx7e//e1YsWJF3H///fHcc8/F9OnTY968eXHUUUdFRMSpp54af/zjH+Oggw7qa4cBgFxUtwAoSleBeK+jjz46jj766JbfmzJlSi9NAzSKG3oAgOrpapVpAAAAqLueKsQAADSDmSxARgIxAAAAbWvSAJpADC0UeZAbIYdiNOniDAAMhmeIAQAASEmFuGZUQAAAAPpDIIaKMNgBdMI5AwB6JxBDcm6qAcjokuMeLbD1kwpsG+gnzxADAACQUu0rxGVVt1TVAIhwPQCAOqt9IAYAoL4MKgFlEogBgLb5nHYAmsQzxAAAAKQkEAMAAJCSKdO0xfM9UD+OWwCAAxOIAQBauGnL2kLbv/jIKwptH4A3JhADpVC9BACgbJ4hBgAAICWBGAAAgJRMmQYAAPqurOfwPf9PJ1SIAQAASEmFGAAABkDlEqpHhRgAAICUVIgB6CsVEAAyKvL653np4gjEVJqDHAAAKIpADABQIQaDAQZHIAZScaMJAMBeAjHAAAji0D3HDwBFsco0AAAAKQnEAAAApGTKNAAApqY3mL8tjE0gBgAgnUf++V8Ka3vWxUsLaxvoL1OmAQAASEmFGFoocmqRaUVQDFMCAWAwmnTNVSEGAAAgJRXimmnSaAwAAECZBOIuCaY0hX0Z6smxCwC9E4ihItzcAgDAYNU+EAsRAAAAdKP2gRgAgPpS3ADKJBDTFhcrACJ8LB0AzeJjlwAAAEhJIAYAACAlU6YBGsqjDgAAByYQA6UQ1gAAKJsp0wAAAKQkEAMAAJCSKdMA0APT/wGgvgRiAIAW/uMfTi+7CxToJxsvKKztWYW1DPSbKdMAAACkJBADAACQkkAMAABASgIxAAAAKVlUCwAA6LtsC9Nle71NIRADAACNIZjSCYEYAKg0N7c0hX0ZqkcgBlJxM1I87zEAUBcCMQAAQA0ZhO6dQAwAUCFucAEGRyCm0twUAAAARRGIAQbA4A5Qdc5TQEbDZXcAAAAAyiAQAwAAkJIp0wAA0GCmw8PYBGJowYUDAACaTyDuUlmBSVBrLn9b6I1jCADolEAMyQkRAAB0okn3jwIxbWnSTg/QBM7LANC72gdiNwQAAPXlXg4oU+0DMQCUKdvNfLbXC0Cz+RxiAAAAUlIhBmgolTyAsd20ZW1hbV985BWFtQ30l0AMlEJYAwCgbKZMAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkNL7sDgAAVNKDNxfb/j9eV2z7ALwhFWIAAABSEogBAABISSAGAAAgJYEYAACAlARiAAAAUrLKNABQbVZ7BqAgKsQAAACkpEIMAAD0X1mzO5q4XTNZCiMQA7mYegkANIX7mp4JxADQCzcjAFBbAjEAAAyCATSoHIGYanPhoCky7cuZXisAUGsCMQAA6fzHP5xedheAChCIaY+KT3P52wIQ4XoApCQQd6uJy7kfaLsAwGC41gMMzHDZHQAAAIAy1L9CbBSVpsi2L2d7vWXI9h5ne73QFI7d4nmP6bcG7VP1D8RQhCIPchcNoB+a+OiO8yMAA2bKNAAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkJJADAAAQErjy+4AAAAM3IM3F9f2P15XXNtAX6kQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASuPL7gAAQBWtP+WxQts/s9DWAWiHCjEAAAApqRADAAD0qMhZJWaUFEcgBkphKiLQLucLAIoiEAMAAH1X1mCWQTQ6UftAbIcfDO9zc/nbNpe/LVAHzlVAmWofiGk2F0n6zT5VPO8xAFAXAnGX3PABANAJ94/0m32qdwJxzdjpoZ4cu0C7nC8ABkcgBoAeCC8AUF8CMVSEm2oAyuQ6BGQkEAPQCG7mB6PI99l7DMCgCcTQghs+AABorUmD0AIxANRQk25GoAyZBr+dL2Bsw2V3AAAAAMogEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKQnEAAAApCQQAwAAkJJADAAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASgIxAAAAKY0vuwNATv/vX39RaPtnzjmt0PYBAKg/FWIAAABSUiEGAKA0ZgwBZap9IHYSBQCgU0XeQ7p/hPqofSAGoDUDhgAAB+YZYgAAAFJSIYbkVBEBAMhKIKYtQhMAANA0pkwDAACQkgoxAPTADBoAqC+BGACgBYMdQCd8lFc9CcRdcpEEgMFwzQU64ZxBJwRiKs0JDQCgntzHUQcCcc04sUBvyjqGMh27mV4rAJTJNbd3AjEA0DbPyEH3hBeoHoEYKsJFEoAI1wOAQfI5xAAAAKSkQgwtmBII9ZOtqpbt9VI8+xTQriadL1SIAQAASEmFGAAAGqxJ1TzoNxViAAAAUhKIAQAASEkgBgAAICWBGAAAgJQEYgAAAFISiAEAAEhJIAYAACAlgRgAAICUBGIAAABSEogBAABISSAGAAAgJYEYAACAlARiAAAAUhKIAQAASEkgBgAAICWBGAAAgJQEYgAAAFISiAEAAEhJIAYAACAlgRgAAICUBGIAAABSEogBAABISSAGAAAgJYEYAACAlARiAAAAUhKIAQAASEkgBgAAICWBGAAAgJQEYgAAAFISiAEAAEhJIAYAACAlgRgAAICUxpe14dHR0YiIePHFF3tqZ89//qUf3RnTWP2zXdut0zZtN+d2M71W223GdjO9Vttt9nYzvVbbHcx2M73WA223mzb25s6xDI2+0U8U5KmnnooZM2aUsWkAAAAS2LFjR0yfPn3M75cWiPfs2RM7d+6MyZMnx9DQUBldAAAAoIFGR0dj9+7dMW3atBgeHvtJ4dICMQAAAJTJoloAAACkJBADAACQkkAMALwhT1gB0EQCMQD8nV/84hexZMmSmDlzZkycODGOOeaYuOiii+L3v/992V0biCeeeCKGhobilltuiYiIXbt2xcc//vF44IEHyu0YABRAIAaA/7Z69eqYO3duPPPMM/HlL385NmzYEMuXL4/NmzfH+973vti+fXvZXSzc2972tvjZz34Wp512WkREbNu2LW677bbYs2dPyT0DgP6zyjQARMRPf/rT+OAHPxif/vSn48Ybb9zve88991wcd9xxccQRR8TDDz9cTgdLsnnz5liwYEHcf//9MX/+/LK7AwB9pUIMABGxatWqOOSQQ+Laa6993fcOP/zwuOGGG+IjH/lIvPzyy/Haa6/FzTffHMcee2xMnDgxZs6cGcuXL49XX3113++cf/75ccopp8Q3v/nNOOaYY2LixIkxb968eOyxx+Luu++OY489Ng4++OA4/vjjY9u2bT3/3vz5818XWDdv3hxDQ0OxefPmiIi45ZZbYvz48fHzn/885s6dGxMmTIi3v/3tcf311+/7nb+fMr03DEdELFiwIObPnx+rV6+OoaGheOyxx/bb1rp162LcuHGxY8eOLv8CADB448vuAACUbXR0NDZt2hRnnHFGHHzwwS1/5qyzztr3709+8pOxdu3aWL58eZx44omxdevW+OIXvxhbt26NjRs3xtDQUEREPPTQQ7Fz58644YYb4pVXXomlS5fGokWLYmhoKK666qqYNGlSXHzxxXHuuefGr3/9633td/t77dizZ0+cddZZsWzZsrjmmmtizZo1ccUVV8Sxxx4bJ5988n4/+973vjdWr14dl1xySaxevTrmz58f06ZNi8svvzxuu+22uPrqq/f97K233hof+tCHYsaMGR31BwDKJBADkN7zzz8fr776ahx11FFv+LO/+c1vYs2aNXHdddfF8uXLIyJi4cKFMW3atFiyZEls2LAhFi1aFBERu3fvju9///vxrne9KyIitmzZEt/4xjfi3nvvjZNOOikiIn7729/G5ZdfHrt27YpDDjmkp99rx+joaKxcuTIuvPDCiIiYN29e/PCHP4y77777dYF4ZGQkZs+eHRERs2fP3vfvM888M9atWxdXXXVVDA0NxVNPPRX33XdfrFu3ru1+AEAVmDINQHrjx/9tfPi11157w5/dsmVLREQsXrx4v6+fffbZMW7cuH3TkyMipkyZsi/URkQcccQRERFx/PHH7/va1KlTI+Jvqzn3+nvtmjt37r5/H3TQQXH44YfHyy+/3PbvX3jhhfHEE0/sW3l67dq1MXny5DjzzDM77gsAlEkgBiC9KVOmxOTJk+PJJ58c82defvnleOGFF+LPf/5zRES89a1v3e/748ePj8MOO2y/gDoyMtKyrUmTJh2wP93+Xrv+97Tw4eHhjlaRPumkk+Koo46KtWvXRsTfAvHZZ58dEyZM6Ev/AGBQBGIAiIiTTz457r///v0Wxvp73/rWt+Kwww7b9/+nn356v+//9a9/jeeff36/nxmkoaGh11W4X3rppcK2df7558f69evj4YcfjkcffTQ+8YlPFLItACiSQAwAEXHZZZfFn/70p1ixYsXrvvf000/H9ddfH7Nnz943LfiOO+7Y72e+973vxWuvvRbvf//7B9Lf/21kZCSeeuqp/b724IMP9tzuuHHjWn79/PPPj127dsVll10Ws2bN2m86NwDUhUW1ACAi5syZE1dffXWsWLEiHnnkkTjvvPPisMMOi1/96lexatWqeOWVV+LHP/5xzJo1K84777xYuXJl/OUvf4kPfOADsW3btrjyyitjwYIFccopp5TS/9NPPz3uuuuuWLZsWZxxxhnxwAMP7JvS3Iu9C3b96Ec/iilTpsR73vOeiIiYOXNmLFy4MDZt2hRf+cpXet4OAJRBhRgA/tvnPve5uOeeeyIi4jOf+UwsWrQovva1r8Xpp58e27Zti1mzZkVExJo1a+ILX/hC3H777bFo0aJYvXp1XHrppXHPPffE8HA5l9YLLrggPvvZz8Z3v/vdWLRoUTz00ENx55139tzuu9/97li8eHHcdNNNce655+73vdNOOy3GjRsXS5Ys6Xk7AFCGodHR0dGyOwEA1M+pp54aEyZMiPXr15fdFQDoiinTAEBHrr766nj00Udj06ZNfXlOGQDKIhADAB2566674ne/+12sWrUqTjjhhLK7AwBdM2UaAACAlCyqBQAAQEoCMQAAACkJxAAAAKQkEAMAAJCSQAwAAEBKAjEAAAApCcQAAACkJBADAACQkkAMAABASv8FKtIjZu66dW0AAAAASUVORK5CYII=",
      "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": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9sAAAHqCAYAAAAUMF39AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkJElEQVR4nO3de5CV9X348c/icr+IgC6XRWRIHf1FEjRRKRLdjhK8UJHEqWKMg2DVNoZMtHGaRsM6JpPJpNp0Wh20cSqxuIZLvKAmRi1LFJE0GgiU0FwcFFAQQYuCCLjP7w/L6orC2eXDHnb39ZphYM8+z3k+Z78+sG+fs+dUFEVRBAAAAJCmU7kHAAAAgPZGbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJKss14EbGhri5Zdfjt69e0dFRUW5xgAAAID9Kooi3nzzzRg8eHB06rT/69Zli+2XX345hg4dWq7DA3AIGjRoUFx11VVxxx13xCuvvFLucdq5EyPi+Yg4KSJ+U+ZZDnGDIuKqiLgjIvxnCdDhrV27Nqqrq/e7Xdliu3fv3hHx3qB9+vQp1xgAHEI2btwY9957byxatCiqqqrKPU6K732v3BN8tA0bOsXdd0dMmfLLGDiwodzjRETEU7H3F+tzY8swyIds2LYh7v7d3THlh1NiYM+BrXrsb37um616PEi3bFnEGWdELFoUMWpUuaeBA7J169YYOnRoY8vuT9lie89Tx/v06SO2AYiIiG3btkVERK9evdrNvw1du5Z7go/Wpcue33sdMjNWxt6DdO1ZhkE+pEvDe1+sLt27RNeerfvFai/nAR1Yr17v/+6/Z9qJUn8M2gukAQAAQDKxDQAAAMnENgAAACQr289sAwCwb7X1tS3br6Zl+wGQx5VtAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSVZZ7AADajrq6ugO+j8mTJydMAgBwaBPbAEC6+qgt9wgAUFaeRg4AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQrLLcAwBAc9XV1R3wfUyePDlhEgCAj+bKNgAAACRzZRuAVrWvq9I7duyIiIif//zn0a1bt9YaCQAgnSvbAAAAkExsAwAAQDJPIwfoIDJeVAwAgNK4sg0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkMz7bANAC9TWlnsCAOBQ5so2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJCsstwDAEA51NXVHdD+K1acECNHjkyaBgBob1zZBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAklWWewAA9q+urq7cIwAA0AyubAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJPM+2wAA7UxtfW3z96lp/j4AfDxXtgEAACCZ2AYAAIBkYhsAAACS+ZltAOBj1UdtuUcAgDbJlW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIFlluQcAAPg49fWlbVdTczCnAIDmc2UbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJJ5NXIAaKEVK1Yc8H2MHDkyYRIA4FDjyjYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQrLLcAwAAUH619bUt26+mZfsBtHeubAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAybz1F8BBVldXV+4RAABoZa5sAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkKzm2b7311ti4cePBnAUAAADahcpSN/zGN74RY8eOjaqqqoiIKIoivvKVr8Q//MM/RHV19UEbEADasxUrVhzwfYwcOTJhEgAgU8lXtouiaPJxQ0ND3HHHHfHqq6+mDwUAAABt2QH9zPaHAxwAAADwAmkAAACQTmwDAABAsmbFdkVFRUm3AQAAQEdW8quRR0SMHj16r9s++9nP7nVbRUVF7N69u+VTAQAAQBtWcmzPmDHjYM4BAAAA7YbYBgAAgGTNehr5B61duzY2bNgQFRUVUV1dHQMHDsycCwAAANqsZsX2zp0745ZbbomZM2fGunXrmnzuE5/4RPzt3/5tXHPNNXHYYYelDglQLnV1deUeAQCANqjk2N6xY0eMGzcuFi9eHKNHj47LL788Bg4cGEVRxLp16+KJJ56Ir3/96/Hwww/Ho48+Gp07dz6YcwMAAMAhq+TY/sEPfhC//vWv44EHHojzzz9/r89/97vfjUceeST+6q/+Kv71X/81vv71r6cOCgAAAG1Fye+zPXfu3Lj22ms/MrT3OO+882L69OmedgkAAECHVnJsv/DCC3H66afvd7szzjgjVq9efUBDAQAAQFtWcmy//fbbcfjhh+93u759+8a2bdsOaCgAAABoy0qO7aIoolOn/W9eUVFxQAMBAABAW1dybEcIaQAAAChFs95n+4ILLoiuXbvuc5t33nnngAYCAACAtq7k2L7ssstc2QYAAIASlBzbd999d0RE/OpXv4oXX3wxRowYESeddNLBmgsASFQftU0+fjMGRcRV8eu4I3rHK2WZCQDas5Jj+4033ogJEybEkiVLGm8bM2ZM3HvvvTF06NCDMhwAAAC0RSW/QNoNN9wQzz//fNx0003xyCOPxC233BKrV6+Oq6666mDOBwAAAG1OyVe2FyxYEN/73vfia1/7WkREnH322TFkyJC45JJLYtu2bdGzZ8+DNiQAAAC0JSVf2d6wYUN85jOfaXJbTU1NvPvuu/HSSy+lDwYAAABtVcmxvWvXrujSpUuT2/r16xcRETt27MidCgAAANqwkmN7X4qiyLgbAAAAaBeaFdsf9z7b3n8bAAAA3lfyC6RFRPzN3/xN9OnTp/HjPVe0r7zyyujdu3fj7RUVFfHkk08mjQgAAABtS8mxffrpp0dFRcVeTxk/44wzIqLpU8k9rRwAAICOrOTYrq+vP4hjAAAAQPuR8gJpAAAAwPvENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkK/mtvwAA4MNq62tbtl9Ny/YDaCtc2QYAAIBkYhsAAACSiW0AAABI5me2AaCNW7FixX632Vi1scnHO7Z0iYiILVs2x/bDNkbVUVUHZTYA6Khc2QYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSVZZ7AAA4lMybd0K5RwAA2gFXtgEAACCZ2AYAAIBkYhsAAACS+ZltoN2qq6sr9whAK6mvL227mpqDOQUAvM+VbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSVZZ7AACg/Da+uvGA76PqqKqESQCgfXBlGwAAAJKJbQAAAEgmtgEAACCZn9kGgDZkVdXt5R4BACiBK9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJCsstwDAADQ8dTW1zZ/n5rm7wNQLq5sAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAECyynIPAPBR6urqyj0CAAC0mNgGADqM+vrStqupOZhTANAReBo5AAAAJBPbAAAAkExsAwAAQDKxDQAAAMnENgAAACQT2wAAAJBMbAMAAEAysQ0AAADJxDYAAAAkqyz3AABA+7Dx1Y0HfB9VR1UlTAIA5efKNgAAACRzZRsAymRV1e3lHgEAOEhc2QYAAIBkrmwD0CHMm3dCuUcAADoQV7YBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZJXlHgBof+rq6so9AgAAlJXYBgCgTaitr23ZfjUt2w/gQHgaOQAAACRzZRsAOGRsfHXjAd9H1VFVCZMAwIFxZRsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwAAAJKJbQAAAEhWWe4BAADgYKqtr23ZfjUt2w8gQmwDwAFbVXV7uUcAAA4xnkYOAAAAycQ2AAAAJBPbAAAAkExsAwAAQDIvkAYA8CH19e//+c2IiIqIXz8X0ftD29XUtNpIALQxrmwDAABAMle2gSbq6urKPQLAAdn46sYDvo+qo6oSJgGgI3NlGwAAAJKJbQAAAEgmtgEAACCZ2AYAAIBkYhsAAACSiW0AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGSV5R4AAA7EvHknlHsE2qGNr25s/POOys0R/SO2bNkc23c3vU6xYsWrH3sfI0eOPGjzAXDoc2UbAAAAkoltAAAASOZp5AAA8BFq62ubv09N8/cB2iexDe1IXV1duUeANm1V1e3lHgEAaCfENgBAC6363VEf+7nNm9//c03NwZ8FgEOLn9kGAACAZK5sAwAcBB98+7B9vUXYvnj7MIC2y5VtAAAASCa2AQAAIJnYBgAAgGRiGwAAAJJ5gTQAAEhSW1/bsv1qWrYfcOgS2wAAh6gVK1Yc8H14RXOA8vA0cgAAAEgmtgEAACCZ2AYAAIBkYhsAAACSeYE0ANqlVVW3l3sEaLTqd0eVtN3/O/7VgzwJAK3FlW0AAABI5so2AIekefNOKPcIAAAtJrYBANox79UNUB5iGw4RdXV15R4BAABIIrYhgVAGAAA+SGwDcEjzquJQfp6KDtB8YhsA4BDhLcI6rtr62pbtV9Oy/YCDT2wD0Kr29Srj3bu/Hn/2Zy/Fk0+OaMWJoO1pi1Hu6jjQ0YhtAABoo1wRh0OX2KZN88JkANBxZFwdz9AerrC3JNIFOjRP2WK7KIqIiNi6dWu5RqDM5s6dW+4RoN14cPODLdpvYv+JyZN8vD0zrj7iyI/d5ojOnePPYlD8ofe98Xq3Xe/duLM1put43t21q/H3d3f6Iu/Lu8UHvla72tbXasXyviVtd9yxmw7uIO3M8796vtwjRETEJz/5yVY9Xou+b3/rrfd/930/bdyec2BPy+5PRVHqlsnWrVsXQ4cOLcehAQAAoEXWrl0b1dXV+92ubLHd0NAQL7/8cvTu3TsqKirKMQIAAACUpCiKePPNN2Pw4MHRqVOn/W5fttgGAACA9mr/OQ4AAAA0S4eM7R07dsS0adOib9++MWjQoLjllls+dtvf/OY3ceqpp0aPHj3i5JNPjueee+4jt5s7d+5eT4dvznEoTWut3f333x8VFRVNfl144YWpj6UjyVy3vn377rU2b/3fi68453K11ro53/Jlrt28efPi2GOPjZ49e8bnP//5ePHFF1t0HPavtdbNOZcra90+vCZ7fv34xz9u9nEoTWutnXMuV9a6FUURtbW1UV1dHUcccURcdNFFsWnT+y8WmXLOFR3QNddcU3zqU58qnnvuueKnP/1p0bt372Lu3Ll7bffWW28VAwcOLK677rpi1apVxfTp04uqqqrirbfearLd66+/XgwcOLD48Jez1ONQutZau+985zvFX/7lXxavvPJK46/XX3/9YD60di1r3datW1dERPGnP/2pydo0NDQ06ziUprXWzfmWL2vtFi9eXFRWVhYzZ84sVq9eXUyYMKEYPXp0s49DaVpr3ZxzubLW7YPr8corrxTXX399MWzYsOKNN95o1nEoXWutnXMuV9a6zZw5s6iuri7q6+uLFStWFGPHji3OP//8Zh9nXzpcbL/11ltFt27dioULFzbedvPNNxdnnHHGXtveddddxfDhwxu/IWxoaCg+8YlPFP/+7//eZLsrrriiOO2005oEW3OOQ2laa+2Koii+9KUvFd/85jezH0KHlLlujz/+eDFo0KADPg7711rrVhTOt2yZazdp0qRiypQpjdu/8MILxbBhw4pNmzY555K11roVhXMu08H43qQo3luzbt26FY8//nizj0NpWmvtisI5lylz3c4///ziuuuua9z+oYceKnr27Nns4+xLh3sa+fLly2PXrl0xZsyYxtvGjh0bS5cujYaGhibbPvvsszF27NjGpxhXVFTEaaedFkuWLGncZtGiRVFfXx/f+ta3WnwcStNaaxcRsWrVqjj22GMP0iPpWDLXbV/r4pzL1VrrVsrnaZ7Mtauvr48vfOELjdsPHz481qxZEwMGDHDOJWutdYtwzmXK/t5kj29/+9tx5plnxllnndXs41Ca1lq7COdcpsx169+/fzzyyCOxfv36ePvtt6Ouri5OPPHEZh9nXzpcbL/yyisxYMCA6NKlS+NtVVVVsWPHjti8efNe2w4ePLjJbVVVVbFu3bqIiHjnnXfiyiuvjNtuuy26d+/e4uNQmtZau6Io4n/+53/isccei2OPPTZGjBgRf//3fx87d+48SI+sfctct9/97nexffv2qKmpiUGDBsW5554bv//975t9HPavtdbN+ZYva+3eeOONeP3112P37t0xfvz4GDhwYEycODHWr1/f7OOwf621bs65XJl/V+7x0ksvxb333hs33nhji45DaVpr7ZxzuTLX7dvf/nZUVlZGdXV19O7dO5566qmoq6tr9nH2pcPF9vbt26Nr165Nbtvz8TvvvFPStnu2u/nmm+Okk06Kz3/+8wd0HErTWmv30ksvNe4/Z86c+Md//MeYPXt2fOMb38h8OB1G5rqtXr06tmzZEjfccEM8+OCD0b179zjzzDPjzTffdM4la611c77ly1q7PS9gN3369Lj00ktjwYIF8c4778SECROioaHBOZestdbNOZcr8+/KPe6666747Gc/G6eeemqLjkNpWmvtnHO5MtdtzZo10aNHj1iwYEEsWrQoqqurY+rUqc0+zr5UlrxlO9GtW7e9vkB7Pu7Ro0dJ2/bo0SNWrlwZd955Z6xYseKAj0NpWmvthg0bFps3b44jjjgiKioqYtSoUdHQ0BCXXnpp3HrrrXHYYYclPqr2L2vdIiJ+/vOfx65du6JXr14RETF79uwYOnRoLFiwwDmXrLXW7ZJLLnG+Jctau8rK975FuOKKK+LLX/5yRLy3dlVVVfHss88655K11rqNGTPGOZco8+/KPebNmxdXX311i49DaVpr7XxfmStr3YqiiMsuuyx+8IMfxIQJEyIiYs6cOTFs2LBYunRp2jnX4a5sDxkyJF577bXYvXt3420bNmyI7t27R9++fffadsOGDU1u27BhQwwaNCjmz58fW7ZsiREjRkSvXr3inHPOiYiIXr16xezZs5t1HErTWmsXEdGvX78mbwd2/PHHx44dO2LLli0H6dG1X1nrFvHe/1HcE2wR7/0lOnz48Fi/fr1zLllrrVuE8y1b1toNGDAgOnfuHMcdd1zj5/r37x/9+/ePtWvXOueStda6RTjnMmX+XRkRsXbt2li1alVMnDixxcehNK21dhHOuUxZ67Zp06ZYu3ZtfPrTn2783NChQ2PAgAHx4osvpp1zHS62R40aFZ07d45nn3228bann346Tj755OjUqemXY/To0fHMM89EURQR8d7PXCxevDhGjx4dX/3qV2P16tWxbNmyWLZsWfzoRz+KiIhly5bF+eef36zjUJrWWrvHHnss+vfvH9u3b2+8v2XLlkX//v3jyCOPbIVH2r5krVtRFDFixIi4++67G7fftm1b/OEPf4jjjjvOOZestdbN+ZYva+0qKyvjM5/5TCxfvrxx+9deey1ee+21OOaYY5xzyVpr3ZxzubLWbY+lS5fG0KFD4+ijj27xcShNa62dcy5X1rr169cvunbtGqtWrWrc/rXXXovNmzfH8OHD8865Zr12eTtx1VVXFZ/85CeLX/3qV8X9999f9OnTp5g/f35RFO+9T9727duLoiiK//3f/y2OPPLIYvr06cV///d/F9OnTy8GDhy413s1F0VRLFy4cK+3j9rXcWiZ1li7rVu3FkOGDCkmT55crF69unj00UeLwYMHF9///vdb50G2Q1nr9tWvfrU4+uiji4ULFxYrV64sJk2aVJxwwgnF7t2793scmq811s35dnBkrd3cuXOLnj17FnPmzClWrVpVTJgwoTjppJMa30bFOZerNdbNOZcv83uTGTNmFOPHj2/2cWiZ1lg751y+rHW7+uqri+HDhxeLFi0qVqxYUYwfP74YM2ZM6r9xHTK2t23bVlx22WVFz549i8GDBxf/9E//1Pi5iGjynnlLly4tTjzxxKJbt27FKaecUjz//PMfeZ8fFdv7Og4t01prt3LlyuKss84qevXqVQwaNKiora1tPPFovqx1e/vtt4trr722GDRoUNGjR49iwoQJxUsvvVTScWi+1lo351u+zL8r77zzzmLYsGFF9+7di3POOadYu3ZtSceh+Vpr3ZxzuTLX7eqrry4uvvjiZh+HlmmttXPO5cr8/uS6664rhgwZUvTr16+46KKLildffbWk45Sq4v+GAgAAAJL4IQ8AAABIJrYBAAAgmdgGAACAZGIbAAAAkoltAAAASCa2AQAAIJnYBgAAgGRiGwBolqIoyj0CABzyxDYAJBg3blz069cvdu7c+bHbjBw5Mk4//fQW3f+aNWuioqIi7r777vR9jjnmmJgyZUpERNTX10dFRUXU19dHRERtbW1UVFQ0brt48eI477zzmjk9AHQ8leUeAADag6lTp8YTTzwRjz76aFxwwQV7ff7555+PlStXxqxZs1p0/4MGDYolS5bEiBEjDnDSvd1///3Rp0+fj/zcFVdcEWeffXbjx//2b/8Wq1atSp8BANobsQ0ACSZNmhRHHHFEzJ49+yNje9asWdGnT5+48MILW3T/Xbt2jdGjRx/glB/txBNP/NjPVVdXR3V19UE5LgC0Z55GDgAJunXrFpdcckk8/PDDsXXr1iaf27VrV9TV1cXkyZNj+/bt8ZWvfCWGDRsWXbp0iX79+sWkSZNizZo1jdvX1NTEpZdeGhdeeGH07Nkzxo0b95FPCf/lL38Z48ePjyOOOCK6dOkSw4cPj9ra2mhoaGhy/PXr18eECROie/fuMXTo0JgxY0a8++67jZ//4NPIP+yDTyOfMmVKzJo1K1588cXGWU4++eQYM2bMXvudddZZMW7cuGZ+FQGg/RDbAJBk6tSpsWPHjpg/f36T23/2s5/Fpk2bYtq0aXHeeefFL37xi/j+978fv/jFL6K2tjaefPLJuPrqq5vs85Of/CR69+4dDz30UFx//fV7HWv58uVx5plnxoABA+InP/lJLFiwID73uc/FTTfdFHPmzGmy7YwZM6KqqioeeOCBmDp1anz3u9+Nv/u7v2v247vxxhvj3HPPjYEDB8aSJUvivPPOi2nTpsWSJUvij3/8Y+N2a9eujYULF8bll1/e7GMAQHvhaeQAkOSkk06KUaNGxezZs5uE5qxZs+JTn/pUDBkyJHr27Bm33HJLjB07NiLeu4r9xz/+Me68884m99WlS5eYOXNmdO3aNSKiyZXviIjf/va3MW7cuLjnnnuiU6f3/t/5uHHj4qGHHor6+vq4+OKLG7c9++yz46677oqIiPHjx8fWrVvj9ttvjxtvvDH69etX8uMbMWJEHHnkkU2e0j558uS49tpr45577ombbropIiLuueee6N27d0yaNKnk+waA9saVbQBINHXq1Fi4cGG8/PLLERGxZcuWePjhh2PatGkxePDg+M///M847bTTYs2aNfH444/Hv/zLv8TixYvjnXfeaXI/xx9/fGNof5Qvf/nL8eijj8bOnTvjt7/9bcyfPz9mzJgRu3fv3uu+LrrooiYff+ELX4idO3fGs88+e8CP9/DDD48vfvGL8R//8R+Nt82aNSsuuuii6N69+wHfPwC0VWIbABJ96Utfis6dO8d9990XERH33XdfVFRUxKWXXhoREbNnz45hw4bF8OHD4+KLL44HH3wwevTosdf99OrVa5/Hefvtt+OKK66Iww8/PEaNGhXXX399rFmzJjp37rzX+2APHDiwycdHHXVURES8/vrrLX6cHzRt2rR44YUX4umnn44lS5bE73//e08hB6DDE9sAkKhfv35xwQUXxL333hsRET/+8Y9j0qRJ0a9fv3j66afjsssuiy9+8Yuxbt262Lx5czzxxBPx53/+580+zte+9rWYN29ezJkzJ958883405/+FPfcc0907tx5r223bNnS5OMNGzZExPvRfaDOOOOMGDFiRMyZMyfmzJkTxx133EF75XQAaCvENgAkmzp1ajz33HOxaNGiWLp0aUybNi0iIp555ploaGiI2traGDJkSEREvPvuu/H4449HROz1KuL78vTTT8df/MVfxMSJE6Nnz54REfHcc8/Fpk2b9rqfRx55pMnH9913X/To0SNOPfXUZj+2ww47bK/bKioq4vLLL48HHnggHnrooY99ZXMA6Ei8QBoAJDvrrLPi6KOPjiuvvDKGDx8eZ555ZkREnHLKKRERcc0118TUqVNjy5Ytcdttt8Xy5csjImLbtm3Ru3fvko5xyimnxJw5c2LmzJlx/PHHx/Lly+M73/lOVFRUxLZt25psO3/+/BgyZEiMGzcuHnvssbjjjjvi5ptvjj59+jT7sfXt2zc2btwYP/vZz2LUqFExaNCgiHjvbcFmzJgREe/9PDkAdHSubANAsk6dOsWUKVMaf3Z5z/tU19TUxG233RbPPPNMnHPOOXHttdfG0UcfHT/96U8jIuKpp54q+Ri33nprTJo0KW644YaYMGFC/OhHP4obbrgh/vqv/zqWLFnS5H20//mf/zn+67/+K84999yYM2dO/PCHP4xvfetbLXpsl19+eRxzzDExceLEmDVrVuPtQ4YMiU9/+tMxfvz4GDx4cIvuGwDak4riw6+iAgDQTOvXr49hw4bF/PnzY+LEieUeBwDKTmwDAC22bNmyePDBB2P+/Pmxe/fuWLlyZeP7fgNAR+ZfQwCgxXbs2BG33npr7N69O+677z6hDQD/x5VtAAAASOZ/PwMAAEAysQ0AAADJxDYAAAAkE9sAAACQTGwDAABAMrENAAAAycQ2AAAAJBPbAAAAkExsAwAAQLL/D6PdVmpY2+2VAAAAAElFTkSuQmCC",
      "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
}