notebook/plot_genomic_structure_hff.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/liuzicheng/anaconda3/envs/evo-design/lib/python3.11/site-packages/cooltools/lib/numutils.py:1376: RuntimeWarning: invalid value encountered in divide\n",
      "  val_cur = ar_cur / armask_cur\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f104510ce90>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAGdCAYAAAAyiFt9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoBklEQVR4nO3df3BV9Z3/8de5F3MDksRQICESDRSERZukgmSywmrXlEA7FNzWQYZZYsZhdyjZkc3qOunUBGt3EpVlcFcGXFwWuiMF3R3s1O82lGYb2H4bBEIzRa2OOPBNEG4CzJiQsCR47/3+AZzLrYnm5uTe8/Hc56PzmSGH8+N9mLbvvN+fzznHikQiEQEAAGP53A4AAAB8PpI1AACGI1kDAGA4kjUAAIYjWQMAYDiSNQAAhiNZAwBgOJI1AACGG5PsC4bDYZ09e1YZGRmyLCvZlwcAOBCJRHTp0iXl5eXJ50tcvXflyhUNDAw4Pk9aWprS09NHISJ3JT1Znz17Vvn5+cm+LABgFHV0dGjq1KkJOfeVK1eUOzZL3XKerHNzc3Xq1KkvfcJOerLOyMiQJJ18skwZgaRf3l1pqTfr0LBgudshuOLtn2e7HYIrJgR73Q4h6cb3OE8oXyZXr/6v/vPAevv/yxNhYGBA3RrQJt/9GusgTf2vPlV18P9qYGCAZB2vG63vjMAYZabfkuzLuysFk3Vg/Di3Q3DFmEBq3vctaWG3Q0i6tFv8bofgimRMY47VGI21HKQpD335IsVKWwDAl4XPL/kc/E7gi0jyyO+PJGsAgJEsnyWfgwreinhnETPJGgBgJL9f8jvIt34PtcFTbxIVAIAvGSprAICRfA7b4D7a4AAAJJbPNwoLzDyCNjgAAIajsgYAGMnnpw1+A8kaAGAkn3WtFT7i4z3yjLVEGxwAAONRWQMAjOTzW/I5WGHmpIVuGpI1AMBIft+1MeLjRy8U19EGBwDAcFTWAAAj+XwO2+DyThucyhoAYCTLf/3LWyMc1gj74Fu2bFFBQYHS09NVUlKiI0eODLnv9u3btXDhQmVnZys7O1tlZWWf2f+xxx6TZVkxY/HixXHFRLIGABjJ54tW1yMb8V9z7969qq6uVl1dnY4fP66ioiKVl5erq6tr0P2bm5u1cuVK/frXv1ZLS4vy8/O1aNEiffzxxzH7LV68WOfOnbPHT3/60/j+LeK/FQAAvGnTpk1as2aNKisrNWfOHG3btk3jxo3Tjh07Bt3/tdde0/e//30VFxdr9uzZevXVVxUOh9XU1BSzXyAQUG5urj2ys7PjiotkDQAw0rXK2tmQpJ6enpjR398/6PUGBgbU2tqqsrKym2LwqaysTC0tLcOK+fLly7p69aomTJgQs725uVmTJ0/WrFmztHbtWl28eDG+f4u49gYAIEl8Puvas9YjHdcXp+Xn5ysrK8se9fX1g17vwoULCoVCysnJidmek5OjYDA4rJiffvpp5eXlxST8xYsX6yc/+Ymampr0/PPP6+DBg1qyZIlCodCw/y1YDQ4A8LSOjg5lZmbaPwcCgYRcp6GhQXv27FFzc7PS09Pt7Y8++qj956997WsqLCzUV7/6VTU3N+uhhx4a1rmprAEARrrxUhQnQ5IyMzNjxlDJeuLEifL7/ers7IzZ3tnZqdzc3M+NdePGjWpoaNAvf/lLFRYWfu6+06dP18SJE3Xy5Mlh/1uQrAEARnLUAr8+4pGWlqa5c+fGLA67sVistLR0yONeeOEFPffcc2psbNS8efO+8DpnzpzRxYsXNWXKlGHHRrIGAOC66upqbd++Xbt27dIf/vAHrV27Vn19faqsrJQkrV69WjU1Nfb+zz//vJ555hnt2LFDBQUFCgaDCgaD6u3tlST19vbqqaee0uHDh3X69Gk1NTVp2bJlmjFjhsrLy4cd14iSdTwPjAMAMBKjtRo8HitWrNDGjRtVW1ur4uJitbW1qbGx0V501t7ernPnztn7b926VQMDA/re976nKVOm2GPjxo2SJL/fr9///vf6zne+o7vuukuPP/645s6dq//5n/+Ja+487gVmNx4Y37Ztm0pKSrR582aVl5frgw8+0OTJk+M9HQAAg3L8utHIyI6tqqpSVVXVoH/X3Nwc8/Pp06c/91xjx47V/v37RxTHzeL+vSPeB8YBAIAzcVXWNx4Yv7lf/0UPjPf398c8gN7T0zPCUAEAqeTGu8FHfPzoheK6uCrrkTwwXl9fH/Mwen5+/sijBQCkDGfvBXfWQjdNwleD19TUqLu72x4dHR2JviQAwAOcfHHL57AqN01cbfCRPDAeCAQS9rYYAABSQVyV9UgfGAcAIF5+y5Lf52BY3mmDx/3oVnV1tSoqKjRv3jzNnz9fmzdvjnlgHACA0TDSZ6VvPt4r4k7WK1as0Pnz51VbW6tgMKji4uKYB8YBAMDoGtFXtz7vgXEAAEbDSN7vHXP8CF+KYiI+kQkAMBJt8CgP3QoAAN5EZQ0AMJLli8jyRRwd7xUkawCAkSzfteHkeK/w0K0AAOBNVNYAACNZVkSW5aAN7uBY05CsAQBGog0eRbIGABjJ8kXkY4GZJOasAQAwHpU1AMBIluWwDe6dF5iRrAEAZuI56yja4AAAGI7KGgBgJFaDR5GsAQBG4kMeUR66FQAAvInKGgBgJN5gFkWyBgAYiTnrKA/dCgAA3kRlDQAwEs9ZR5GsAQBGog0eRbIGABjp2qNbI6+OeXQLAAAkDZU1AMBItMGjSNYAACNZcvictVhg5lya79pIIUd3XnU7hKTbENjndgiu+LulD7sdgitO/Dzb7RCSbsynYbdDSKqrV1Pv/8dMQGUNADASbfAokjUAwEg8Zx3lod87AADwJiprAICRLIefyKQNDgBAgtEGj/LQ7x0AAHgTlTUAwEisBo8iWQMAjGRZDl+K4uBY05CsAQBG8jlcYMaHPAAAQNJQWQMAjMRq8CiSNQDASCwwi/LQrQAA4E1U1gAAM/ksyW85O94jSNYAACNZPkuWg4Tr5FjT0AYHAMBwVNYAADP5fdeGk+M9gmQNADCTz3I270wbHAAAJAuVNQDASJZfshysBrf8oxiMy0jWAAAz0Qa3kawBAGbyO3zO2smxhmHOGgAAw1FZAwCMZFkOX4pieaeyJlkDAMzEc9Y279wJAAAeFXeyPnTokJYuXaq8vDxZlqU333wzAWEBAFLdjXeDOxleEXey7uvrU1FRkbZs2ZKIeAAAuObGanAnYwS2bNmigoICpaenq6SkREeOHBly3+3bt2vhwoXKzs5Wdna2ysrKPrN/JBJRbW2tpkyZorFjx6qsrEwffvhhXDHFnayXLFmiH//4x3r44YfjPRQAAKPt3btX1dXVqqur0/Hjx1VUVKTy8nJ1dXUNun9zc7NWrlypX//612ppaVF+fr4WLVqkjz/+2N7nhRde0D/90z9p27Ztevvtt3XrrbeqvLxcV65cGXZcCZ+z7u/vV09PT8wAAOALuVBZb9q0SWvWrFFlZaXmzJmjbdu2ady4cdqxY8eg+7/22mv6/ve/r+LiYs2ePVuvvvqqwuGwmpqaJF2rqjdv3qwf/vCHWrZsmQoLC/WTn/xEZ8+ejWsaOeHJur6+XllZWfbIz89P9CUBAB6Q7DnrgYEBtba2qqyszN7m8/lUVlamlpaWYZ3j8uXLunr1qiZMmCBJOnXqlILBYMw5s7KyVFJSMuxzSklI1jU1Neru7rZHR0dHoi8JAIDtj7u7/f39g+534cIFhUIh5eTkxGzPyclRMBgc1rWefvpp5eXl2cn5xnFOziklIVkHAgFlZmbGDAAAvtAotcHz8/NjOrz19fUJCbehoUF79uzRvn37lJ6ePqrn5qUoAAAzWT7J56CmtK4d29HREVMoBgKBQXefOHGi/H6/Ojs7Y7Z3dnYqNzf3cy+1ceNGNTQ06Fe/+pUKCwvt7TeO6+zs1JQpU2LOWVxcPOxbiftfobe3V21tbWpra5N0rR/f1tam9vb2eE8FAMCQLL/leEj6THd3qGSdlpamuXPn2ovDJNmLxUpLS4eM84UXXtBzzz2nxsZGzZs3L+bvpk2bptzc3Jhz9vT06O233/7cc/6xuCvrY8eO6Rvf+Ib9c3V1tSSpoqJCO3fujPd0AAAYo7q6WhUVFZo3b57mz5+vzZs3q6+vT5WVlZKk1atX6/bbb7db6c8//7xqa2u1e/duFRQU2PPQ48eP1/jx42VZltavX68f//jHmjlzpqZNm6ZnnnlGeXl5Wr58+bDjijtZP/jgg4pEIvEeBgBAfFz4nvWKFSt0/vx51dbWKhgMqri4WI2NjfYCsfb2dvluas1v3bpVAwMD+t73vhdznrq6Om3YsEGS9Pd///fq6+vTX/3VX+mTTz7RggUL1NjYGNe8NnPWAAAzufQ966qqKlVVVQ36d83NzTE/nz59+gvPZ1mWfvSjH+lHP/rRiOKR+JAHAADGo7IGABjJ6cc4vPQhD5I1AMBMfM/a5p07AQDAo6isAQBm8svhArNRi8R1JGsAgJEsy+GcteWdOWva4AAAGI7KGgBgJpeeszYRyRoAYCYX3mBmKpI1AMBIN3+MY6THewVz1gAAGI7KGgBgJp/D71k7OdYwJGsAgJmYs7Z559cOAAA8isoaAGAm2uA2kjUAwEwka5t37gQAAI+isgYAmMlyuMDMQ+8GJ1kDAMxEG9zmnTsBAMCjqKwBAGaisraRrAEAZuKlKDbXknXDguUKjB/n1uVdsSGwz+0Qku7/vTXgdgiuqD32726H4Ip/+du/dDuEpHurebLbISTVp1fSpTeSdDGf5bCy9k6y9k6PAAAAj6INDgAwE3PWNpI1AMBMzFnbvPNrBwAAHkVlDQAwE21wG8kaAGAmkrXNO3cCAIBHUVkDAIxkWZYsBx/jcHKsaUjWAAAzWQ7b4JZ3msfeuRMAADyKyhoAYCYWmNlI1gAAM/FSFBvJGgBgJiprm3fuBAAAj6KyBgCYiTa4jWQNADAT37O20QYHAMBwVNYAADOxwMxGsgYAmIk5a5t3fu0AAMCjqKwBAGaiDW4jWQMAzMSHPGwkawCAmZiztnnn1w4AADyKyhoAYCbmrG0kawCAmSyfs3lnD81Ze+dOAADwqLiSdX19ve677z5lZGRo8uTJWr58uT744INExQYASGU3KmsnwyPiupODBw9q3bp1Onz4sA4cOKCrV69q0aJF6uvrS1R8AIBUdWPO2snwiLjmrBsbG2N+3rlzpyZPnqzW1lb92Z/92agGBgAArnG0wKy7u1uSNGHChCH36e/vV39/v/1zT0+Pk0sCAFKFZTlcYMZz1gqHw1q/fr3uv/9+3XPPPUPuV19fr6ysLHvk5+eP9JIAgFTCnLVtxHeybt06vfPOO9qzZ8/n7ldTU6Pu7m57dHR0jPSSAACkpBG1wauqqvTWW2/p0KFDmjp16ufuGwgEFAgERhQcACCF8Zy1La5kHYlE9Dd/8zfat2+fmpubNW3atETFBQBIdbzBzBZXsl63bp12796tn/3sZ8rIyFAwGJQkZWVlaezYsQkJEACQoqisbXHdydatW9Xd3a0HH3xQU6ZMscfevXsTFR8AACkv7jY4AABJQWVt40MeAAAzkaxt3rkTAABGwZYtW1RQUKD09HSVlJToyJEjQ+777rvv6rvf/a4KCgpkWZY2b978mX02bNggy7JixuzZs+OKiWQNADCU0xeixJ/i9u7dq+rqatXV1en48eMqKipSeXm5urq6Bt3/8uXLmj59uhoaGpSbmzvkee+++26dO3fOHr/5zW/iiotkDQAwkwsf8ti0aZPWrFmjyspKzZkzR9u2bdO4ceO0Y8eOQfe/77779OKLL+rRRx/93HeKjBkzRrm5ufaYOHFiXHGRrAEAntbT0xMzbv5exc0GBgbU2tqqsrIye5vP51NZWZlaWlocxfDhhx8qLy9P06dP16pVq9Te3h7X8SRrAICRLMsny/I7GNdSXH5+fsw3Kurr6we93oULFxQKhZSTkxOzPScnx36vyEiUlJRo586damxs1NatW3Xq1CktXLhQly5dGvY5WA0OADDTKK0G7+joUGZmpr052a/AXrJkif3nwsJClZSU6M4779Trr7+uxx9/fFjnIFkDADwtMzMzJlkPZeLEifL7/ers7IzZ3tnZ+bmLx+J122236a677tLJkyeHfQxtcACAmZL8icy0tDTNnTtXTU1N9rZwOKympiaVlpaO2m319vbqo48+0pQpU4Z9DJU1AMBMLrwUpbq6WhUVFZo3b57mz5+vzZs3q6+vT5WVlZKk1atX6/bbb7fnvQcGBvTee+/Zf/7444/V1tam8ePHa8aMGZKkJ598UkuXLtWdd96ps2fPqq6uTn6/XytXrhx2XCRrAICZXPjq1ooVK3T+/HnV1tYqGAyquLhYjY2N9qKz9vZ2+W4679mzZ/X1r3/d/nnjxo3auHGjHnjgATU3N0uSzpw5o5UrV+rixYuaNGmSFixYoMOHD2vSpEnDjotkDQDATaqqqlRVVTXo391IwDcUFBR84Xcz9uzZ4zgmkjUAwEy8G9xGsgYAmIlkbfPOnQAA4FFU1gAAM1FZ20jWAAAzubAa3FTeuRMAADyKyhoAYCbLctgGt0YvFpeRrAEAZmLO2uadOwEAwKOorAEAZqKytpGsAQBmIlnbXEvWb/88W2MC49y6vCv+bunDboeQdLXH/t3tEFzxf1676nYIrlh7yy63Q0i6S39d4XYISdXfd1mtybqYdX04Od4jvPNrBwAAHkUbHABgpEgk8oVftPqi472CZA0AMFJEYUUUdnS8V9AGBwDAcFTWAAAjRa7/x8nxXkGyBgAYKRIJKxJx0AZ3cKxpaIMDAGA4KmsAgJFog0eRrAEARrr26JaTNrh3kjVtcAAADEdlDQAwEm3wKJI1AMBIvBQlimQNADASj25FMWcNAIDhqKwBAEZizjqKZA0AMBJt8Cja4AAAGI7KGgBgJNrgUSRrAICReHQrijY4AACGo7IGABjp2rvBHbTBPfRucJI1AMBItMGjaIMDAGA4KmsAgJFYDR5FsgYAGImXokSRrAEARorIWXXsnbo6zjnrrVu3qrCwUJmZmcrMzFRpaal+8YtfJCo2AACgOCvrqVOnqqGhQTNnzlQkEtGuXbu0bNky/e53v9Pdd9+dqBgBAKnIYRtcqdoGX7p0aczP//AP/6CtW7fq8OHDJGsAwKhigVnUiOesQ6GQ3njjDfX19am0tHTI/fr7+9Xf32//3NPTM9JLAgCQkuJO1idOnFBpaamuXLmi8ePHa9++fZozZ86Q+9fX1+vZZ591FCQAIPXwUpSouF+KMmvWLLW1tentt9/W2rVrVVFRoffee2/I/WtqatTd3W2Pjo4ORwEDAFLDjdeNOhleEXdlnZaWphkzZkiS5s6dq6NHj+qll17SK6+8Muj+gUBAgUDAWZQAAKQwx89Zh8PhmDlpAABGA23wqLiSdU1NjZYsWaI77rhDly5d0u7du9Xc3Kz9+/cnKj4AQIriq1tRcSXrrq4urV69WufOnVNWVpYKCwu1f/9+ffOb30xUfAAApLy4kvW//uu/JioOAABi0AaP4t3gAAAjhSPXhpPjvYJkDQAwUihiKRSxHB3vFXE/Zw0AAJKLyhoAYCTa4FEkawCAkcIRS2EHrWwnx5qGNjgAAIajsgYAGCkckUK0wSWRrAEAhqINHkUbHAAAw1FZAwCMxGrwKJI1AMBIvBQlijY4AACGo7IGABgpfH04Od4rqKwBAEYKy7JXhI9oaGRt8C1btqigoEDp6ekqKSnRkSNHhtz33Xff1Xe/+10VFBTIsixt3rzZ8TkHQ7IGABjpxgIzJyNee/fuVXV1terq6nT8+HEVFRWpvLxcXV1dg+5/+fJlTZ8+XQ0NDcrNzR2Vcw6GZA0AwHWbNm3SmjVrVFlZqTlz5mjbtm0aN26cduzYMej+9913n1588UU9+uijCgQCo3LOwZCsAQBGurEa3MmQpJ6enpjR398/6PUGBgbU2tqqsrIye5vP51NZWZlaWlpGdA+jdU6SNQDASBGHLfDI9TZ4fn6+srKy7FFfXz/o9S5cuKBQKKScnJyY7Tk5OQoGgyO6h9E6J6vBAQCe1tHRoczMTPvnodrVJiNZAwCMNFrvBs/MzIxJ1kOZOHGi/H6/Ojs7Y7Z3dnYOuXgsWeekDQ4AMFIo4nzEIy0tTXPnzlVTU5O9LRwOq6mpSaWlpSO6h9E6J5U1AADXVVdXq6KiQvPmzdP8+fO1efNm9fX1qbKyUpK0evVq3X777fa898DAgN577z37zx9//LHa2to0fvx4zZgxY1jnHA6SNQDASG58InPFihU6f/68amtrFQwGVVxcrMbGRnuBWHt7u3y+aFP67Nmz+vrXv27/vHHjRm3cuFEPPPCAmpubh3XO4SBZAwCM5NZXt6qqqlRVVTXo391IwDcUFBQoEvniC33eOYfDtWQ9IdirW9K89ObWL3bi59luh5B0//K3f+l2CK5Ye8sut0Nwxdu/TL1lMM9O3+t2CEnVc+WqXnI7iBREZQ0AMNJIFon98fFeQbIGABjJjTlrU5GsAQBGcmvO2kSpN8EEAMCXDJU1AMBI4Zs+xjHS472CZA0AMBJt8Cja4AAAGI7KGgBgJFaDR5GsAQBG4jnrKNrgAAAYjsoaAGCkawvMnLTBRzEYl5GsAQBGYjV4FG1wAAAMR2UNADASC8yiSNYAACOF5DBZj1ok7iNZAwCMFHE4Zx3xUGXNnDUAAIajsgYAGIk56yiSNQDASCTrKNrgAAAYjsoaAGAkKusokjUAwEi8wSyKNjgAAIajsgYAGIk2eBTJGgBgpLDDZE0bHAAAJI2jZN3Q0CDLsrR+/fpRCgcAgGtuLDBzMrxixG3wo0eP6pVXXlFhYeFoxgMAgCTmrG82osq6t7dXq1at0vbt25WdnT3aMQEAYCdrJ8MrRpSs161bp29/+9sqKyv7wn37+/vV09MTMwAAwPDF3Qbfs2ePjh8/rqNHjw5r//r6ej377LNxBwYASG2hiKVQxHJ0vFfEVVl3dHToiSee0Guvvab09PRhHVNTU6Pu7m57dHR0jChQAEBqYYFZVFyVdWtrq7q6unTvvffa20KhkA4dOqSXX35Z/f398vv9MccEAgEFAoHRiRYAgBQUV7J+6KGHdOLEiZhtlZWVmj17tp5++unPJGoAAEaK1eBRcSXrjIwM3XPPPTHbbr31Vn3lK1/5zHYAAJzgDWZRvMEMAADDOX43eHNz8yiEAQBArHD42nByvFfwIQ8AgJGYs46iDQ4AgOGorAEARgrJYWU9apG4j2QNADCS0xebeGk1OMkaAGAk5qyjmLMGAMBwVNYAACNRWUeRrAEARuINZlG0wQEAMByVNQDASKwGjyJZAwCMxJx1FG1wAAAMR2UNADBSOGwpHLYcHe8VJGsAgJHCIUvhkINk7eBY09AGBwDAcFTWAAAj0QaPIlkDAIxEso4iWQMAjESyjmLOGgAAw1FZAwCMFAk7Ww0e8VBlTbIGABiJNngUbXAAAAxHZQ0AMBKVdZRryXp8z4DSbvG7dXlXjPk07HYISfdW82S3Q3DFpb+ucDsEVzw7fa/bISRd+0G3I0iu3k+Td61w+NpwcvxIbNmyRS+++KKCwaCKior0z//8z5o/f/6Q+7/xxht65plndPr0ac2cOVPPP/+8vvWtb9l//9hjj2nXrl0xx5SXl6uxsXHYMdEGBwDgur1796q6ulp1dXU6fvy4ioqKVF5erq6urkH3/+1vf6uVK1fq8ccf1+9+9zstX75cy5cv1zvvvBOz3+LFi3Xu3Dl7/PSnP40rLpI1AMBIN94N7mTEa9OmTVqzZo0qKys1Z84cbdu2TePGjdOOHTsG3f+ll17S4sWL9dRTT+lP/uRP9Nxzz+nee+/Vyy+/HLNfIBBQbm6uPbKzs+OKi2QNADDSjTlrJyMeAwMDam1tVVlZmb3N5/OprKxMLS0tgx7T0tISs790rcX9x/s3Nzdr8uTJmjVrltauXauLFy/GFRsLzAAAntbT0xPzcyAQUCAQ+Mx+Fy5cUCgUUk5OTsz2nJwcvf/++4OeOxgMDrp/MBi0f168eLH+4i/+QtOmTdNHH32kH/zgB1qyZIlaWlrk9w9v7RbJGgBgpNFaDZ6fnx+zva6uThs2bHASWlweffRR+89f+9rXVFhYqK9+9atqbm7WQw89NKxzkKwBAEYKhSz5HLzBLHT92I6ODmVmZtrbB6uqJWnixIny+/3q7OyM2d7Z2anc3NxBj8nNzY1rf0maPn26Jk6cqJMnTw47WTNnDQAwUjjicM46ci1ZZ2ZmxoyhknVaWprmzp2rpqamaAzhsJqamlRaWjroMaWlpTH7S9KBAweG3F+Szpw5o4sXL2rKlCnD/rcgWQMAcF11dbW2b9+uXbt26Q9/+IPWrl2rvr4+VVZWSpJWr16tmpoae/8nnnhCjY2N+sd//Ee9//772rBhg44dO6aqqipJUm9vr5566ikdPnxYp0+fVlNTk5YtW6YZM2aovLx82HHRBgcAGCnicM56JB/yWLFihc6fP6/a2loFg0EVFxersbHRXkTW3t4uny9a5/7pn/6pdu/erR/+8If6wQ9+oJkzZ+rNN9/UPffcI0ny+/36/e9/r127dumTTz5RXl6eFi1apOeee27ICn8wJGsAgJHcet1oVVWVXRn/sebm5s9se+SRR/TII48Muv/YsWO1f//+EcVxM9rgAAAYjsoaAGCkkb6F7ObjvYJkDQAwEl/diqINDgCA4aisAQBGuvaJTCeV9SgG4zKSNQDASLTBo2iDAwBgOCprAICRWA0eRbIGABjpxrvBnRzvFSRrAICRmLOOYs4aAADDUVkDAIwUcThnHWHOGgCAxAqHLVm0wSXRBgcAwHhU1gAAM4Uj14aT4z2CZA0AMJIvHJHPQcKNeChZx9UG37BhgyzLihmzZ89OVGwAAEAjqKzvvvtu/epXv4qeYAzFOQBg9FmhiKzQyKtjJ8eaJu5MO2bMGOXm5iYiFgAAbJbDNng4VdvgkvThhx8qLy9P06dP16pVq9Te3v65+/f396unpydmAACA4YsrWZeUlGjnzp1qbGzU1q1bderUKS1cuFCXLl0a8pj6+nplZWXZIz8/33HQAADv80Ui9iKzEY1IilbWS5Ys0SOPPKLCwkKVl5frv/7rv/TJJ5/o9ddfH/KYmpoadXd326Ojo8Nx0AAA77PCEcfDKxytDrvtttt011136eTJk0PuEwgEFAgEnFwGAJCCfCHJ52CRmC80isG4zNEbzHp7e/XRRx9pypQpoxUPAAD4I3El6yeffFIHDx7U6dOn9dvf/lYPP/yw/H6/Vq5cmaj4AAApytF8tcOV5KaJqw1+5swZrVy5UhcvXtSkSZO0YMECHT58WJMmTUpUfACAFOU04aZsst6zZ0+i4gAAAEPg9WMAACM5XdHNanAAABKMNngU37MGAMBwVNYAACP5QhGHz1l7p7ImWQMAjHTjdaNOjvcK2uAAABiOyhoAYCan7/f20AIzkjUAwEisBo8iWQMAjMQCsyjmrAEAMByVNQDASLzBLIpkDQAwki8cli8cdnS8V9AGBwDAcFTWAAAjsRo8imQNADASq8GjaIMDAGA4KmsAgJGsiMPV4B56NzjJGgBgJOaso2iDAwBgOCprAICRqKyjSNYAACP5QhH5fKwGl0jWAABThR2+MtQ7LzBjzhoAANMlvbKOXF9Kf/Xq/yb70q67evWq2yEk3adX0t0OwRX9fZfdDsEVPVdS77/jvZ+6HUFy9X567YYjSXgsKjRwWZ86qKxDn3onz1iRZPyL3+TMmTPKz89P5iUBAKOso6NDU6dOTci5r1y5omnTpikYDDo+V25urk6dOqX09C934ZD0ZB0Oh3X27FllZGTIsqykXbenp0f5+fnq6OhQZmZm0q7rNu47de47Fe9ZSs37dvOeI5GILl26pLy8PPl8iZtJvXLligYGBhyfJy0t7UufqCUX2uA+ny9hv40NR2ZmZsr8D/pm3HfqSMV7llLzvt2656ysrIRfIz093RNJdrSwwAwAAMORrAEAMFzKJOtAIKC6ujoFAgG3Q0kq7jt17jsV71lKzftOxXtOdUlfYAYAAOKTMpU1AABfViRrAAAMR7IGAMBwJGsAAAyXMsl6y5YtKigoUHp6ukpKSnTkyBG3Q0qoQ4cOaenSpcrLy5NlWXrzzTfdDinh6uvrdd999ykjI0OTJ0/W8uXL9cEHH7gdVsJt3bpVhYWF9gsySktL9Ytf/MLtsJKqoaFBlmVp/fr1boeSUBs2bJBlWTFj9uzZboeFJEiJZL13715VV1errq5Ox48fV1FRkcrLy9XV1eV2aAnT19enoqIibdmyxe1QkubgwYNat26dDh8+rAMHDujq1atatGiR+vr63A4toaZOnaqGhga1trbq2LFj+vM//3MtW7ZM7777rtuhJcXRo0f1yiuvqLCw0O1QkuLuu+/WuXPn7PGb3/zG7ZCQDJEUMH/+/Mi6devsn0OhUCQvLy9SX1/vYlTJIymyb98+t8NIuq6uroikyMGDB90OJemys7Mjr776qtthJNylS5ciM2fOjBw4cCDywAMPRJ544gm3Q0qourq6SFFRkdthwAWer6wHBgbU2tqqsrIye5vP51NZWZlaWlpcjAyJ1t3dLUmaMGGCy5EkTygU0p49e9TX16fS0lK3w0m4devW6dvf/nbM/7697sMPP1ReXp6mT5+uVatWqb293e2QkARJ/5BHsl24cEGhUEg5OTkx23NycvT++++7FBUSLRwOa/369br//vt1zz33uB1Owp04cUKlpaW6cuWKxo8fr3379mnOnDluh5VQe/bs0fHjx3X06FG3Q0makpIS7dy5U7NmzdK5c+f07LPPauHChXrnnXeUkZHhdnhIIM8na6SmdevW6Z133kmZ+bxZs2apra1N3d3d+o//+A9VVFTo4MGDnk3YHR0deuKJJ3TgwIGU+jLTkiVL7D8XFhaqpKREd955p15//XU9/vjjLkaGRPN8sp44caL8fr86Oztjtnd2dio3N9elqJBIVVVVeuutt3To0CFXP8eaTGlpaZoxY4Ykae7cuTp69KheeuklvfLKKy5Hlhitra3q6urSvffea28LhUI6dOiQXn75ZfX398vv97sYYXLcdtttuuuuu3Ty5Em3Q0GCeX7OOi0tTXPnzlVTU5O9LRwOq6mpKSXm9FJJJBJRVVWV9u3bp//+7//WtGnT3A7JNeFwWP39/W6HkTAPPfSQTpw4oba2NnvMmzdPq1atUltbW0okaknq7e3VRx99pClTprgdChLM85W1JFVXV6uiokLz5s3T/PnztXnzZvX19amystLt0BKmt7c35rftU6dOqa2tTRMmTNAdd9zhYmSJs27dOu3evVs/+9nPlJGRoWAwKEnKysrS2LFjXY4ucWpqarRkyRLdcccdunTpknbv3q3m5mbt37/f7dASJiMj4zNrEW699VZ95Stf8fQahSeffFJLly7VnXfeqbNnz6qurk5+v18rV650OzQkWEok6xUrVuj8+fOqra1VMBhUcXGxGhsbP7PozEuOHTumb3zjG/bP1dXVkqSKigrt3LnTpagSa+vWrZKkBx98MGb7v/3bv+mxxx5LfkBJ0tXVpdWrV+vcuXPKyspSYWGh9u/fr29+85tuh4ZRdubMGa1cuVIXL17UpEmTtGDBAh0+fFiTJk1yOzQkGJ/IBADAcJ6fswYA4MuOZA0AgOFI1gAAGI5kDQCA4UjWAAAYjmQNAIDhSNYAABiOZA0AgOFI1gAAGI5kDQCA4UjWAAAYjmQNAIDh/j9pxukyeVJGtAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "orca_path='genomic_structure/train/log/hff_orca/weight/model_hff_orcahff_a_orca20240514_063101.checkpoint'\n",
    "hyena_path='genomic_structure/train/log/hyena/weight/model_hyenahff20240513_133816.checkpoint'\n",
    "mamba_path='genomic_structure/train/log/mamba/weight/model_mambahff_a20240513_133536.checkpoint'\n",
    "import sys\n",
    "\n",
    "from selene_utils2 import *\n",
    "from matplotlib import pyplot as plt\n",
    "batch_size=5\n",
    "max_length=6000\n",
    "dest_path='data/genomic_structure/resources/resources'\n",
    "t = Genomic2DFeatures(\n",
    "            [dest_path + \"/4DNFI643OYP9.rebinned.mcool::/resolutions/1000\"],\n",
    "            [\"r1000\"],\n",
    "            (int(max_length/1000), int(max_length/1000)),\n",
    "            cg=True,\n",
    "        )\n",
    "\n",
    "sampler = RandomPositionsSamplerHiC(\n",
    "            reference_sequence=MemmapGenome(\n",
    "                input_path=dest_path + \"/Homo_sapiens.GRCh38.dna.primary_assembly.fa\",\n",
    "                memmapfile=dest_path + \"/Homo_sapiens.GRCh38.dna.primary_assembly.fa.mmap\",\n",
    "                blacklist_regions=\"hg38\",\n",
    "            ),\n",
    "            target=t,\n",
    "            # target_1d=MultibinGenomicFeatures(\n",
    "            #     self.dest_path + \"/h1esc/h1esc.hg38.bed.sorted.gz\",\n",
    "            #     np.loadtxt(self.dest_path + \"/h1esc/h1esc.hg38.bed.sorted.features\", str),\n",
    "            #     1000,\n",
    "            #     1000,\n",
    "            #     (32, 10),\n",
    "            #     mode=\"any\",\n",
    "            # ),\n",
    "            features=[\"r1000\"],\n",
    "            test_holdout=[\"chr9\", \"chr10\"],\n",
    "            validation_holdout=[\"chr8\"],\n",
    "            sequence_length=6000,\n",
    "            position_resolution=1000,\n",
    "            random_shift=100,\n",
    "            random_strand=False,\n",
    "            cross_chromosome=False,\n",
    "        )\n",
    "sampler.mode = \"validate\"\n",
    "for i in range(10):\n",
    "    sampled_data = sampler.sample(batch_size=batch_size)\n",
    "\n",
    "    while np.isnan(sampled_data[1]).any():\n",
    "                sampled_data = sampler.sample(batch_size=batch_size)\n",
    "\n",
    "# print(sampled_data)\n",
    "plt.imshow(sampled_data[1][0, :, :],cmap = 'Spectral_r')   \n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequence=sampled_data[0]\n",
    "target=sampled_data[1]\n",
    "target=torch.from_numpy(target).float().cuda()\n",
    "sequence_b=np.zeros((sequence.shape[0],sequence.shape[1]))\n",
    "mask=np.where(sequence==0.25)\n",
    "sequence_b=sequence.argmax(axis=2)\n",
    "sequence_b[[mask[0],mask[1]]]=4\n",
    "sequence=sequence_b\n",
    "map_to_genomic_sequence = {0: 'A', 1: 'C', 2: 'G', 3: 'T', 4: 'N'}\n",
    "#map sequence to genome \n",
    "sequence = \"\".join(map(lambda x: map_to_genomic_sequence[x], sequence.flatten().tolist()))\n",
    "#reshape list\n",
    "# sequence = np.array(list(sequence)).reshape(batch_size,-1)\n",
    "from scipy.stats import pearsonr\n",
    "def genomic_structure_corr(logit,y,ignore_index=-100):\n",
    "\n",
    "\n",
    "    pred=logit\n",
    "    target=y\n",
    "    normmat_bydist = np.exp(\n",
    "        np.load(\"orca/resources/resources/4DNFI643OYP9.rebinned.mcool.expected.res1000.npy\")        \n",
    "\n",
    "    )[:10]\n",
    "    normmat = normmat_bydist[np.abs(np.arange(6)[:, None] - np.arange(6)[None, :])]\n",
    "    normmat_r = torch.from_numpy(normmat).float().cuda()\n",
    "    eps = torch.min(normmat_r)\n",
    "    target = torch.log(((target + eps) / (normmat_r + eps)))\n",
    "    corr=[]\n",
    "    #convert to numpy\n",
    "    target=target.cpu().numpy().reshape((pred.shape[0], -1))\n",
    "    pred=pred.detach().cpu().numpy().reshape((pred.shape[0], -1))\n",
    "    for j in range(pred.shape[0]):\n",
    "            if np.mean(np.isnan(target[j, :])) < 0.7:\n",
    "                corr.append(\n",
    "                    pearsonr(\n",
    "                        pred[j, ~np.isnan(target[j, :])],\n",
    "                        target[j, ~np.isnan(target[j, :])],\n",
    "                    )[0]\n",
    "                )\n",
    "            else:\n",
    "                corr.append(np.nan)\n",
    "    corr=np.nanmean(corr)\n",
    "    return corr\n",
    "def genomic_structure_mse(logit,y,ignore_index=-100):\n",
    "\n",
    "    pred=logit\n",
    "    target=y\n",
    "    if torch.isnan(target).all():\n",
    "        target=torch.rand_like(target)\n",
    "\n",
    "    normmat_bydist = np.exp(\n",
    "        np.load(\"orca/resources/resources/4DNFI643OYP9.rebinned.mcool.expected.res1000.npy\")\n",
    "\n",
    "    )[:8]\n",
    "    normmat = normmat_bydist[np.abs(np.arange(6)[:, None] - np.arange(6)[None, :])]\n",
    "    normmat_r = torch.from_numpy(normmat).float().cuda()\n",
    "    eps = torch.min(normmat_r)\n",
    "    target = torch.nanmean(\n",
    "                torch.nanmean(torch.reshape(target, (target.shape[0], 6, 1, 6, 1)), axis=4),\n",
    "                axis=2,\n",
    "            )\n",
    "    target_r = torch.log(((target + eps) / (normmat_r + eps)))\n",
    "    target_cuda = target_r\n",
    "\n",
    "    loss = (\n",
    "        (\n",
    "            pred[~torch.isnan(target)]\n",
    "            - target_cuda[~torch.isnan(target)]\n",
    "        )\n",
    "        ** 2\n",
    "    ).mean()\n",
    "\n",
    "\n",
    "\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/liuzicheng/anaconda3/envs/evo-design/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n",
      "/liuzicheng/anaconda3/envs/evo-design/lib/python3.11/site-packages/torch/utils/checkpoint.py:460: UserWarning: torch.utils.checkpoint: please pass in use_reentrant=True or use_reentrant=False explicitly. The default value of use_reentrant will be updated to be False in the future. To maintain current behavior, pass use_reentrant=True. It is recommended that you use use_reentrant=False. Refer to docs for more details on the differences between the two variants.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.42186854902787496\n",
      "tensor(0.0905, device='cuda:0', grad_fn=<MeanBackward0>)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f0fd9639690>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGdCAYAAAAv9mXmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUVklEQVR4nO3dXYyVhb3v8f8wk1n4Moyi8jJlQI0VtpKhEYXMsbZWqWa24WgvGkI4KYeanrQZGggxaeam6EUz3OwGo4SSvnFTgrY7aGKqlNICaSoVh5CATY0YGsfwps3uzDAeF+yZdS5O9nSzBTcL+K/HNfP5JE/SWXkWz2+lyNe1HmZsqFQqlQCAq2xS0QMAGJ8EBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFI01fqCo6Ojcfz48WhpaYmGhoZaXx6AK1CpVGJoaCja2tpi0qRPf49S88AcP3482tvba31ZAK6i/v7+mDVr1qeeU/PAtLS0RETEg9P+dzRNaq715Qv13r9MLXpCzVUOTil6QiEW/vNbRU8oxPGPWoueUHNt1w4UPaGmzg2fi3/9ny+O/Vn+aWoemP/4WKxpUvOEC0zjtaWiJ9RcpTS56AmFaL5+Yv3e/g9NDRPv93jzdRPz/+tLucXhJj8AKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKS4rMBs2rQpbr311pg8eXIsXrw43njjjau9C4A6V3VgXnjhhVi3bl2sX78+Dh48GAsWLIhHH300Tp8+nbEPgDpVdWB++MMfxre+9a1YtWpV3HXXXfGjH/0orr322vjZz36WsQ+AOlVVYM6ePRt9fX2xZMmSf/wCkybFkiVL4vXXX7/gc8rlcgwODp53ADD+VRWYDz/8MEZGRmL69OnnPT59+vQ4efLkBZ/T29sbra2tY0d7e/vlrwWgbqT/LbKenp4YGBgYO/r7+7MvCcBnQFM1J998883R2NgYp06dOu/xU6dOxYwZMy74nFKpFKVS6fIXAlCXqnoH09zcHAsXLozdu3ePPTY6Ohq7d++Ozs7Oqz4OgPpV1TuYiIh169bFypUr4957741FixbFxo0bY3h4OFatWpWxD4A6VXVgli1bFh988EF8//vfj5MnT8YXvvCFeO211z5x4x+Aia3qwERErF69OlavXn21twAwjvhZZACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEjRVNSF3/uXqdF4bamoyxei/X8dK3pCzTVce03REwqx5465RU8oxP3/dLToCTW3oKW/6Ak19XHDv8f2SzzXOxgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACmqDsy+ffti6dKl0dbWFg0NDfHSSy8lzAKg3lUdmOHh4ViwYEFs2rQpYw8A40RTtU/o6uqKrq6ujC0AjCNVB6Za5XI5yuXy2NeDg4PZlwTgMyD9Jn9vb2+0traOHe3t7dmXBOAzID0wPT09MTAwMHb09/dnXxKAz4D0j8hKpVKUSqXsywDwGeP7YABIUfU7mDNnzsTRo0fHvj527FgcOnQopk6dGrNnz76q4wCoX1UH5s0334yvfOUrY1+vW7cuIiJWrlwZW7duvWrDAKhvVQfmwQcfjEqlkrEFgHHEPRgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEjRVNSFKwenRKU0uajLF6Lh2muKnlBzIx/+regJhTjWtbvoCYV47t/mFD2BZI1RueRzvYMBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIUVVgent747777ouWlpaYNm1aPPHEE/H2229nbQOgjlUVmL1790Z3d3fs378/du3aFefOnYtHHnkkhoeHs/YBUKeaqjn5tddeO+/rrVu3xrRp06Kvry++9KUvXdVhANS3qgLzXw0MDERExNSpUy96TrlcjnK5PPb14ODglVwSgDpx2Tf5R0dHY+3atXH//ffH/PnzL3peb29vtLa2jh3t7e2Xe0kA6shlB6a7uzuOHDkS27dv/9Tzenp6YmBgYOzo7++/3EsCUEcu6yOy1atXxyuvvBL79u2LWbNmfeq5pVIpSqXSZY0DoH5VFZhKpRLf/e53Y8eOHbFnz5647bbbsnYBUOeqCkx3d3ds27YtXn755WhpaYmTJ09GRERra2tcc801KQMBqE9V3YPZvHlzDAwMxIMPPhgzZ84cO1544YWsfQDUqao/IgOAS+FnkQGQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIqmoi688J/fiubrm4u6fCH23DG36Ak1d6xrd9ETCnHbS/+n6AmFKH1Q2B8phZn34LtFT6ipc8NnI2LvJZ3rHQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBRVBWbz5s3R0dERU6ZMiSlTpkRnZ2e8+uqrWdsAqGNVBWbWrFmxYcOG6OvrizfffDMeeuihePzxx+Ott97K2gdAnWqq5uSlS5ee9/UPfvCD2Lx5c+zfvz/uvvvuqzoMgPpWVWD+s5GRkfjlL38Zw8PD0dnZedHzyuVylMvlsa8HBwcv95IA1JGqb/IfPnw4rr/++iiVSvHtb387duzYEXfddddFz+/t7Y3W1taxo729/YoGA1Afqg7M3Llz49ChQ/GnP/0pvvOd78TKlSvjz3/+80XP7+npiYGBgbGjv7//igYDUB+q/oisubk57rjjjoiIWLhwYRw4cCCeffbZ2LJlywXPL5VKUSqVrmwlAHXnir8PZnR09Lx7LAAQUeU7mJ6enujq6orZs2fH0NBQbNu2Lfbs2RM7d+7M2gdAnaoqMKdPn45vfOMbceLEiWhtbY2Ojo7YuXNnfPWrX83aB0CdqiowP/3pT7N2ADDO+FlkAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASNFU1IWPf9QaTQ2loi5fiPv/6WjRE2ruuX+bU/SEQpQ+KOwfrUJtXPHToifUXM/GJ4ueUFMjZz++5HO9gwEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkOKKArNhw4ZoaGiItWvXXqU5AIwXlx2YAwcOxJYtW6Kjo+Nq7gFgnLiswJw5cyZWrFgRP/7xj+PGG2+82psAGAcuKzDd3d3x2GOPxZIlS/7bc8vlcgwODp53ADD+NVX7hO3bt8fBgwfjwIEDl3R+b29vPPPMM1UPA6C+VfUOpr+/P9asWRO/+MUvYvLkyZf0nJ6enhgYGBg7+vv7L2soAPWlqncwfX19cfr06bjnnnvGHhsZGYl9+/bF888/H+VyORobG897TqlUilKpdHXWAlA3qgrMww8/HIcPHz7vsVWrVsW8efPie9/73ifiAsDEVVVgWlpaYv78+ec9dt1118VNN930iccBmNh8Jz8AKar+W2T/1Z49e67CDADGG+9gAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEgRVNRF267diCar2su6vKFWNDSX/QEamTeg+8WPaEQPRufLHpCzQ39j4+KnlBTox99HLHl0s71DgaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIqqAvP0009HQ0PDece8efOytgFQx5qqfcLdd98dv/3tb//xCzRV/UsAMAFUXYempqaYMWNGxhYAxpGq78G888470dbWFrfffnusWLEi3nvvvU89v1wux+Dg4HkHAONfVYFZvHhxbN26NV577bXYvHlzHDt2LB544IEYGhq66HN6e3ujtbV17Ghvb7/i0QB89lUVmK6urvj6178eHR0d8eijj8avf/3r+Pvf/x4vvvjiRZ/T09MTAwMDY0d/f/8Vjwbgs++K7tDfcMMNceedd8bRo0cvek6pVIpSqXQllwGgDl3R98GcOXMm3n333Zg5c+bV2gPAOFFVYJ566qnYu3dv/PWvf40//vGP8bWvfS0aGxtj+fLlWfsAqFNVfUT2/vvvx/Lly+Nvf/tb3HLLLfHFL34x9u/fH7fcckvWPgDqVFWB2b59e9YOAMYZP4sMgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSNNX6gpVKJSIizg2fq/WlC/dxw78XPaHmGqNS9IRCnBs+W/SEQoyc/bjoCTU3+tHEes2j/7ccEf/4s/zTNFQu5ayr6P3334/29vZaXhKAq6y/vz9mzZr1qefUPDCjo6Nx/PjxaGlpiYaGhppdd3BwMNrb26O/vz+mTJlSs+sWzeueOK97Ir7miIn5uot8zZVKJYaGhqKtrS0mTfr0uyw1/4hs0qRJ/231Mk2ZMmXC/Cb8z7zuiWMivuaIifm6i3rNra2tl3Sem/wApBAYAFJMmMCUSqVYv359lEqloqfUlNc9cV73RHzNERPzddfLa675TX4AJoYJ8w4GgNoSGABSCAwAKQQGgBQTJjCbNm2KW2+9NSZPnhyLFy+ON954o+hJqfbt2xdLly6Ntra2aGhoiJdeeqnoSel6e3vjvvvui5aWlpg2bVo88cQT8fbbbxc9K93mzZujo6Nj7JvuOjs749VXXy16Vk1t2LAhGhoaYu3atUVPSfX0009HQ0PDece8efOKnnVREyIwL7zwQqxbty7Wr18fBw8ejAULFsSjjz4ap0+fLnpamuHh4ViwYEFs2rSp6Ck1s3fv3uju7o79+/fHrl274ty5c/HII4/E8PBw0dNSzZo1KzZs2BB9fX3x5ptvxkMPPRSPP/54vPXWW0VPq4kDBw7Eli1boqOjo+gpNXH33XfHiRMnxo4//OEPRU+6uMoEsGjRokp3d/fY1yMjI5W2trZKb29vgatqJyIqO3bsKHpGzZ0+fboSEZW9e/cWPaXmbrzxxspPfvKTomekGxoaqnz+85+v7Nq1q/LlL3+5smbNmqInpVq/fn1lwYIFRc+4ZOP+HczZs2ejr68vlixZMvbYpEmTYsmSJfH6668XuIxsAwMDERExderUgpfUzsjISGzfvj2Gh4ejs7Oz6Dnpuru747HHHjvvn+/x7p133om2tra4/fbbY8WKFfHee+8VPemiav7DLmvtww8/jJGRkZg+ffp5j0+fPj3+8pe/FLSKbKOjo7F27dq4//77Y/78+UXPSXf48OHo7OyMjz/+OK6//vrYsWNH3HXXXUXPSrV9+/Y4ePBgHDhwoOgpNbN48eLYunVrzJ07N06cOBHPPPNMPPDAA3HkyJFoaWkpet4njPvAMDF1d3fHkSNHPtufT19Fc+fOjUOHDsXAwED86le/ipUrV8bevXvHbWT6+/tjzZo1sWvXrpg8eXLRc2qmq6tr7H93dHTE4sWLY86cOfHiiy/Gk08+WeCyCxv3gbn55pujsbExTp06dd7jp06dihkzZhS0ikyrV6+OV155Jfbt21fofxqilpqbm+OOO+6IiIiFCxfGgQMH4tlnn40tW7YUvCxHX19fnD59Ou65556xx0ZGRmLfvn3x/PPPR7lcjsbGxgIX1sYNN9wQd955Zxw9erToKRc07u/BNDc3x8KFC2P37t1jj42Ojsbu3bsnxGfUE0mlUonVq1fHjh074ne/+13cdtttRU8qzOjoaJTL5aJnpHn44Yfj8OHDcejQobHj3nvvjRUrVsShQ4cmRFwiIs6cORPvvvtuzJw5s+gpFzTu38FERKxbty5WrlwZ9957byxatCg2btwYw8PDsWrVqqKnpTlz5sx5/1Zz7NixOHToUEydOjVmz55d4LI83d3dsW3btnj55ZejpaUlTp48GRH//z+OdM011xS8Lk9PT090dXXF7NmzY2hoKLZt2xZ79uyJnTt3Fj0tTUtLyyfurV133XVx0003jet7bk899VQsXbo05syZE8ePH4/169dHY2NjLF++vOhpF1b0X2Orleeee64ye/bsSnNzc2XRokWV/fv3Fz0p1e9///tKRHziWLlyZdHT0lzo9UZE5ec//3nR01J985vfrMyZM6fS3NxcueWWWyoPP/xw5Te/+U3Rs2puIvw15WXLllVmzpxZaW5urnzuc5+rLFu2rHL06NGiZ12UH9cPQIpxfw8GgGIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0CK/wdVMXz7KGeWlQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sys\n",
    "\n",
    "from src.models.sequence.orca import Orca\n",
    "orca_model=Orca().to('cuda')\n",
    "seq_ids =torch.from_numpy(sampled_data[0]).squeeze(0)\n",
    "seq_ids = seq_ids.float().to('cuda')\n",
    "target = torch.from_numpy(sampled_data[1]).squeeze(0)\n",
    "target = target.float().to('cuda')\n",
    "\n",
    "checkpoint=torch.load(orca_path)\n",
    "\n",
    "torch.nn.modules.utils.consume_prefix_in_state_dict_if_present(\n",
    "            checkpoint, \"module.\"\n",
    "        )\n",
    "\n",
    "orca_model.load_state_dict(checkpoint,strict=False)\n",
    "\n",
    "prediction,_=orca_model(seq_ids)\n",
    "\n",
    "corr=genomic_structure_corr(prediction,target)\n",
    "print(corr)\n",
    "mse=genomic_structure_mse(prediction,target)\n",
    "print(mse)\n",
    "\n",
    "#plot prediction as image\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.imshow(prediction.cpu().detach().numpy()[0])\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.31035648576963776\n",
      "tensor(0.1013, device='cuda:0')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f0fcfb54910>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGdCAYAAAAPGjobAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1CUlEQVR4nO3df3xU1Z3/8fdMQhKUzMRAyJgSRLZoggh8DRJGaVEJBqVWlEXkwRaIecBju4lFo27BUqntdmO3VUGhoC1KXU2hWKULtXHTIFAl4UeQKhRYtWoiOAlIk0CAJGbm+wcydUoImXszM5eZ17OP86C5OWfu50blk8+5555r8/l8PgEAgKhhj3QAAACgZ5HcAQCIMiR3AACiDMkdAIAoQ3IHACDKkNwBAIgyJHcAAKIMyR0AgCgTH+4Ter1eHTp0SMnJybLZbOE+PQDABJ/Pp2PHjikjI0N2e+jqw1OnTqmtrc305yQkJCgpKakHIrqwhD25Hzp0SJmZmeE+LQCgB9XV1WnAgAEh+exTp07J1dupJplP7i6XSx9++GHMJfiwJ/fk5GRJ0p7pY5WcEPbTR1T1g9MjHULYpfVuj3QIEfH/aj+IdAiREYu7WX/WGOkIwqr5RJsum13m/7s8FNra2tSkNj1hv169TaSpk/pcJZ631NbWRnIPtTNT8ckJ8XLEWHK/KLl3pEMIuz69Y+uf8RmOPomRDiEyvDGY3E8mRDqCiAjHbdXeildvm4m/Q2LwX8czYvNvXgCA5dnjJLuJ3yHsPkneHgvngkJyBwBYks1uk93EDIHNF7uLtknuAABLiouT4kzk57gYnpbnOXcAAKIMlTsAwJLsJqfl7UzLAwBgLXZ7Dyyoi1FMywMAEGWo3AEAlmSPY1reKJI7AMCS7LbTU/OGx8foM+4S0/IAAEQdKncAgCXZ42yym1hRZ2ZK/0JHcgcAWFKc/XQzPL7nQrngMC0PAECUoXIHAFiS3W5yWl5MywMAYCm2OHOr5WP4ljvJHQBgTXa7TFbusSuWrx0AgKhEcgcAWNLpyt1cM2LZsmUaNGiQkpKSlJubq+3bt3fZf+3atcrKylJSUpKuvvpqvfbaa2f12bdvn775zW/K6XTq4osv1rXXXqva2lpjAXYDyR0AYEl2u+30s+5Gm4Ep/TVr1qikpESLFi3Srl27NGLECOXn56uhoaHT/lu3btX06dNVWFiot99+W5MnT9bkyZO1Z88ef58PPvhAY8eOVVZWljZt2qR33nlH3//+95WUlGT4Z3M+Np/PF9b35jQ3N8vpdOrjWTfIkRBbt/zf/N7MSIcQdum92yMdQkTkfPxepEOIDG8MvobryN8iHUFYNZ9o0yV3rVJTU5McDkdozvFFnlg/eIIujutl+HNaOtp1218rgoo1NzdX1157rZYuXSpJ8nq9yszM1L333qv58+ef1X/atGlqaWnRhg0b/MfGjBmjkSNHasWKFZKku+++W7169dJ///d/G76WYFG5AwAs6cwmNmaadPqXhS+31tbWTs/X1tammpoa5eXl+Y/Z7Xbl5eWpqqqq0zFVVVUB/SUpPz/f39/r9er3v/+9rrjiCuXn56t///7Kzc3VunXrzP+AukByBwBYkqkp+S+aJGVmZsrpdPpbaWlpp+c7cuSIOjo6lJ6eHnA8PT1dHo+n0zEej6fL/g0NDTp+/Lgee+wxTZw4Uf/7v/+rO+64Q3feeac2b95s9kd0TrE1Lw4AiDl1dXUB0/KJiYlhO7fXe/rVdLfffrvuv/9+SdLIkSO1detWrVixQuPGjQvJeQ1V7sGuJAQAIFg9tVre4XAEtHMl9379+ikuLk719fUBx+vr6+VyuTod43K5uuzfr18/xcfHa+jQoQF9srOzrbVaPtiVhAAAGHFm+1kzLRgJCQnKyclRZWWl/5jX61VlZaXcbnenY9xud0B/SaqoqPD3T0hI0LXXXqsDBw4E9Pm///s/XXbZZUHFF4ygk/sTTzyhOXPmqKCgQEOHDtWKFSt00UUX6bnnngtFfAAAhE1JSYl+8Ytf6Fe/+pX27dunb3/722ppaVFBQYEkaebMmVqwYIG//7x581ReXq7HH39c+/fv1w9+8APt3LlTxcXF/j4PPfSQ1qxZo1/84hd6//33tXTpUq1fv17/9m//FrLrCOqe+5mVhF++sPOtJGxtbQ1Ymdjc3GwwVABALLHFSXYT7201snHttGnTdPjwYT3yyCPyeDwaOXKkysvL/YvmamtrZf/S7jjXXXedysrKtHDhQj388MMaMmSI1q1bp2HDhvn73HHHHVqxYoVKS0v1ne98R1deeaV++9vfauzYscYv7jyCSu5drSTcv39/p2NKS0v16KOPGo8QABCTTL8VzmdsbHFxcUDl/WWbNm0669jUqVM1derULj/znnvu0T333GMoHiNC/ijcggUL1NTU5G91dXWhPiUAIArY48y3WBVU5W5kJWFiYmJYHzsAACDWBVW5G1lJCACAEXE2m+LsJloMv9A96E1sSkpKNGvWLI0aNUqjR4/W4sWLA1YSAgDQE8y82e3M+FgVdHI/30pCAAAQWYa2n+1qJSEAAD3hy/vDGxpvcLV8NGBveQCAJTEtb1wMXzoAANGJyh0AYEk2u082u8/U+FhFcgcAWJLNfrqZGR+rYvjSAQCITlTuAABLstl8stlMTMubGHuhI7kDACyJaXnjSO4AAEuy2X2ys6DOkBj+vQYAgOhE5Q4AsCSbzeS0fOxuUEdyBwBYE8+5G8e0PAAAUYbKHQBgSayWN47kDgCwJF4cY1wMXzoAANGJyh0AYEnsUGccyR0AYEncczcuhi8dAIDoROUOALAknnM3juQOALAkpuWNI7kDACzp9KNwxqtvHoUDAABRg8odAGBJTMsbR3IHAFiSTSafcxcL6sKu+sHpuii5d6ROHxFfvzQp0iGE3b13RTqCyFi16upIhxARvoPvRzqEsLMNviLSIYSV7dhJSasiHQbOg8odAGBJTMsbF8OXDgCwsjPPuZtpRixbtkyDBg1SUlKScnNztX379i77r127VllZWUpKStLVV1+t11577Zx9//Vf/1U2m02LFy82FFt3kdwBAPjCmjVrVFJSokWLFmnXrl0aMWKE8vPz1dDQ0Gn/rVu3avr06SosLNTbb7+tyZMna/LkydqzZ89ZfV999VVVV1crIyMj1JdBcgcAWJPN/vfXvhppRqbln3jiCc2ZM0cFBQUaOnSoVqxYoYsuukjPPfdcp/2XLFmiiRMn6qGHHlJ2drZ+9KMf6ZprrtHSpUsD+h08eFD33nuvXnrpJfXq1cvIjyMoJHcAgCX11LR8c3NzQGttbe30fG1tbaqpqVFeXp7/mN1uV15enqqqqjodU1VVFdBfkvLz8wP6e71efetb39JDDz2kq666yuyPpVtI7gCAqJaZmSmn0+lvpaWlnfY7cuSIOjo6lJ6eHnA8PT1dHo+n0zEej+e8/X/yk58oPj5e3/nOd0xeSfexWh4AYEk9tVq+rq5ODofDfzwxMdFkZN1XU1OjJUuWaNeuXbLZbGE7L5U7AMCSbDaf6SZJDocjoJ0ruffr109xcXGqr68POF5fXy+Xy9XpGJfL1WX/P/3pT2poaNDAgQMVHx+v+Ph4ffzxx3rggQc0aNAgkz+hcyO5AwAsycxiujMtGAkJCcrJyVFlZaX/mNfrVWVlpdxud6dj3G53QH9Jqqio8Pf/1re+pXfeeUe7d+/2t4yMDD300EN6/fXXgwswCEzLAwDwhZKSEs2aNUujRo3S6NGjtXjxYrW0tKigoECSNHPmTH3lK1/x37efN2+exo0bp8cff1yTJk3S6tWrtXPnTj377LOSpL59+6pv374B5+jVq5dcLpeuvPLKkF0HyR0AYElmNqI5Mz5Y06ZN0+HDh/XII4/I4/Fo5MiRKi8v9y+aq62tlf1LUwLXXXedysrKtHDhQj388MMaMmSI1q1bp2HDhhmOuyeQ3AEAlhSp7WeLi4tVXFzc6fc2bdp01rGpU6dq6tSp3f78jz76yFhgQeCeOwAAUYbKHQBgTXabFGfi8TF7+B49sxqSOwDAkmx2m2wmErSZsRc6puUBAIgyVO4AAGuKs59uZsbHKJI7AMCa7DZz982ZlgcAANGCyh0AYEm2OMlmYrW8La4Hg7nAkNwBANbEtLxhJHcAgDXFmXzO3czYCxz33AEAiDJU7gAAS7LZTG5iY4vdyp3kDgCwJp5zNyx2rxwAgCgVdHLfsmWLbrvtNmVkZMhms2ndunUhCAsAEOvO7C1vpsWqoJN7S0uLRowYoWXLloUiHgAATjuzWt5Mi1FB33O/5ZZbdMstt4QiFgAA0ANCvqCutbVVra2t/q+bm5tDfUoAQDTgOXfDQr6grrS0VE6n098yMzNDfUoAQBTgnrtxIU/uCxYsUFNTk7/V1dWF+pQAAMS0kE/LJyYmKjExMdSnAQBEG6blDWMTGwCANdnskt3EBLMtdrdyCTq5Hz9+XO+//77/6w8//FC7d+9WamqqBg4c2KPBAQBily3OZvKVr1Tu3bZz507deOON/q9LSkokSbNmzdKqVat6LDAAAGBM0Mn9hhtukM/nC0UsAAD8He9zN4x77gAAa2JBnWGxu9oAAIAoReUOALAksxvRxPImNiR3AIA18T53w2L3ygEAiFJU7gAAa4qTyQV1PRbJBYfkDgCwJJvN5D13W+zec2daHgCAL1m2bJkGDRqkpKQk5ebmavv27V32X7t2rbKyspSUlKSrr75ar732mv977e3t+u53v6urr75aF198sTIyMjRz5kwdOnQopNdAcgcAWNOZ59zNtCCtWbNGJSUlWrRokXbt2qURI0YoPz9fDQ0NnfbfunWrpk+frsLCQr399tuaPHmyJk+erD179kiSTpw4oV27dun73/++du3apVdeeUUHDhzQN7/5TVM/mvOx+cK83Vxzc7OcTqfWvPuMLkruHc5TR9zXL02KdAhhd+9dpyIdQkSsWhWbN/t8B98/f6coY+s/INIhhFXzsZNK+adiNTU1yeFwhOYcX+SJo6umy3FRgvHPOdGm1Nm/DirW3NxcXXvttVq6dKkkyev1KjMzU/fee6/mz59/Vv9p06appaVFGzZs8B8bM2aMRo4cqRUrVnR6jh07dmj06NH6+OOPQ/ZOFip3AIAlnXlxjJkWjLa2NtXU1CgvL89/zG63Ky8vT1VVVZ2OqaqqCugvSfn5+efsL0lNTU2y2WxKSUkJKr5gsKAOABDVmpubA75OTExUYmLiWf2OHDmijo4OpaenBxxPT0/X/v37O/1sj8fTaX+Px9Np/1OnTum73/2upk+fHrKZD4nKHQBgVXa7+SYpMzNTTqfT30pLSyNyOe3t7brrrrvk8/m0fPnykJ6Lyh0AYE099Fa4urq6gCq5s6pdkvr166e4uDjV19cHHK+vr5fL5ep0jMvl6lb/M4n9448/1saNG0NatUtU7gCAKOdwOALauZJ7QkKCcnJyVFlZ6T/m9XpVWVkpt9vd6Ri32x3QX5IqKioC+p9J7O+9957++Mc/qm/fvj1wVV2jcgcAWNOXptYNjw9SSUmJZs2apVGjRmn06NFavHixWlpaVFBQIEmaOXOmvvKVr/in9ufNm6dx48bp8ccf16RJk7R69Wrt3LlTzz77rKTTif2f//mftWvXLm3YsEEdHR3++/GpqalKSDD+NEBXSO4AAGuKQHKfNm2aDh8+rEceeUQej0cjR45UeXm5f9FcbW2t7F/63Ouuu05lZWVauHChHn74YQ0ZMkTr1q3TsGHDJEkHDx7U//zP/0iSRo4cGXCuN954QzfccIOxazsPkjsAAF9SXFys4uLiTr+3adOms45NnTpVU6dO7bT/oEGDFObtZCSR3AEAVmUzuaAuhveWJ7kDAKwpAtPy0SJ2rxwAgChF5Q4AsCYqd8NI7gAAa+qhTWxiUcSSe1rvdvXpHVu/W9x7V6QjCL+HnmuNdAgRUT1uXaRDiIjcJwZHOoTwS0mNdATh1XYyfOey20xW7rGb3GN3zgIAgCgVW6UzAODCwT13w0juAABr4p67YbH7aw0AAFGKyh0AYE1MyxtGcgcAWBPJ3bDYvXIAAKIUlTsAwJJsNptsJl7+YmbshY7kDgCwJpvJaXlb7E5Ox+6VAwAQpajcAQDWxII6w0juAABrYhMbw0juAABronI3LHavHACAKEXlDgCwJqblDSO5AwCsife5G8a0PAAAUYbKHQBgTSyoM4zkDgCwJu65Gxa7v9YAABClqNwBANbEtLxhJHcAgDXx4hjDSO4AAGvinrthsftrDQAAUYrKHQBgTdxzN4zkDgCwJpvd3H3zGL7nHrtXDgBAlAoquZeWluraa69VcnKy+vfvr8mTJ+vAgQOhig0AEMvOVO5mmgHLli3ToEGDlJSUpNzcXG3fvr3L/mvXrlVWVpaSkpJ09dVX67XXXgv4vs/n0yOPPKJLL71UvXv3Vl5ent577z1DsXVXUFe+efNmFRUVqbq6WhUVFWpvb9fNN9+slpaWUMUHAIhVZ+65m2lBWrNmjUpKSrRo0SLt2rVLI0aMUH5+vhoaGjrtv3XrVk2fPl2FhYV6++23NXnyZE2ePFl79uzx9/mv//ovPfXUU1qxYoW2bdumiy++WPn5+Tp16pThH835BHXl5eXlmj17tq666iqNGDFCq1atUm1trWpqakIVHwAAYfPEE09ozpw5Kigo0NChQ7VixQpddNFFeu655zrtv2TJEk2cOFEPPfSQsrOz9aMf/UjXXHONli5dKul01b548WItXLhQt99+u4YPH64XXnhBhw4d0rp160J2HabuuTc1NUmSUlNTz9mntbVVzc3NAQ0AgPOy2UxOy59+zv0fc1Bra2unp2tra1NNTY3y8vL8x+x2u/Ly8lRVVdXpmKqqqoD+kpSfn+/v/+GHH8rj8QT0cTqdys3NPedn9gTDyd3r9eq+++7T9ddfr2HDhp2zX2lpqZxOp79lZmYaPSUAIJb00D33zMzMgDxUWlra6emOHDmijo4OpaenBxxPT0+Xx+PpdIzH4+my/5k/g/nMnmD4UbiioiLt2bNHb775Zpf9FixYoJKSEv/Xzc3NJHgAQNjU1dXJ4XD4v05MTIxgNOFhKLkXFxdrw4YN2rJliwYMGNBl38TExJj4QQIAelgPPefucDgCkvu59OvXT3Fxcaqvrw84Xl9fL5fL1ekYl8vVZf8zf9bX1+vSSy8N6DNy5MhuX0qwgvqp+Xw+FRcX69VXX9XGjRt1+eWXhyouAECsC/Nq+YSEBOXk5KiystJ/zOv1qrKyUm63u9Mxbrc7oL8kVVRU+PtffvnlcrlcAX2am5u1bdu2c35mTwiqci8qKlJZWZl+97vfKTk52X+/wOl0qnfv3iEJEAAQoyKwQ11JSYlmzZqlUaNGafTo0Vq8eLFaWlpUUFAgSZo5c6a+8pWv+O/bz5s3T+PGjdPjjz+uSZMmafXq1dq5c6eeffbZ0yHYbLrvvvv0H//xHxoyZIguv/xyff/731dGRoYmT55s/NrOI6jkvnz5cknSDTfcEHD8+eef1+zZs3sqJgAAImLatGk6fPiwHnnkEXk8Ho0cOVLl5eX+BXG1tbWyf2lG4LrrrlNZWZkWLlyohx9+WEOGDNG6desCFpr/+7//u1paWjR37lw1NjZq7NixKi8vV1JSUsiuw+bz+Xwh+/RONDc3y+l0auP7S9UnObaq/aVze0U6hLB76LnOHzmJdsduXBfpECIi94nBkQ4h7GxZV0Q6hLBqPnZSKdn/rqampm7dxzZ0ji/yRONfl8phIk80HzuplMHFIY3VqnhxDADAmnhxjGGxe+UAAEQpKncAgEWZrNxjuH4luQMArMngy18Cxseo2L1yAACiFJU7AMCSbDa7bLY4U+NjFckdAGBNrJY3LHavHACAKEXlDgCwJip3w0juAABrIrkbRnIHAFgTj8IZFrtXDgBAlKJyBwBYE9PyhpHcAQDWRHI3LHavHACAKEXlDgCwJip3w0juAABrYrW8YbF75QAARCkqdwCANdlsJqflbT0XywWG5A4AsCbuuRsWu1cOAECUonIHAFgTlbthJHcAgDWR3A2LWHL/f7UfyNEnMVKnj4hVq66OdAhhVz1uXaRDiIiR3/RFOoSIsI25IdIhhN1ncccjHUJYHevVK3wns33RzIyPUbH7aw0AAFGKaXkAgCX5fD75fMZnwcyMvdCR3AEAluSTVz55TY2PVUzLAwAQZajcAQCW5Pvif2bGxyqSOwDAknw+r3w+E9PyJsZe6JiWBwAgylC5AwAsiWl546jcAQCWdPpROK+JFtrkfvToUc2YMUMOh0MpKSkqLCzU8eNdb2p06tQpFRUVqW/fvurTp4+mTJmi+vp6//f//Oc/a/r06crMzFTv3r2VnZ2tJUuWBB0byR0AAANmzJihvXv3qqKiQhs2bNCWLVs0d+7cLsfcf//9Wr9+vdauXavNmzfr0KFDuvPOO/3fr6mpUf/+/fXiiy9q7969+t73vqcFCxZo6dKlQcXGtDwAwJKsPC2/b98+lZeXa8eOHRo1apQk6emnn9att96qn/3sZ8rIyDhrTFNTk1auXKmysjLddNNNkqTnn39e2dnZqq6u1pgxY3TPPfcEjBk8eLCqqqr0yiuvqLi4uNvxUbkDACzpzCY2ZpokNTc3B7TW1lbTsVVVVSklJcWf2CUpLy9Pdrtd27Zt63RMTU2N2tvblZeX5z+WlZWlgQMHqqqq6pznampqUmpqalDxkdwBAJZk7n773x+jy8zMlNPp9LfS0lLTsXk8HvXv3z/gWHx8vFJTU+XxeM45JiEhQSkpKQHH09PTzzlm69atWrNmzXmn+/8R0/IAgKhWV1cnh8Ph/zox8dxvJJ0/f75+8pOfdPl5+/bt67HYurJnzx7dfvvtWrRokW6++eagxpLcAQCW1FP33B0OR0By78oDDzyg2bNnd9ln8ODBcrlcamhoCDj++eef6+jRo3K5XJ2Oc7lcamtrU2NjY0D1Xl9ff9aYv/zlLxo/frzmzp2rhQsXdiv2LyO5AwAsKRI71KWlpSktLe28/dxutxobG1VTU6OcnBxJ0saNG+X1epWbm9vpmJycHPXq1UuVlZWaMmWKJOnAgQOqra2V2+3299u7d69uuukmzZo1Sz/+8Y+DvgaJe+4AAAQtOztbEydO1Jw5c7R9+3a99dZbKi4u1t133+1fKX/w4EFlZWVp+/btkiSn06nCwkKVlJTojTfeUE1NjQoKCuR2uzVmzBhJp6fib7zxRt18880qKSmRx+ORx+PR4cOHg4qPyh0AYElWfhROkl566SUVFxdr/PjxstvtmjJlip566in/99vb23XgwAGdOHHCf+zJJ5/0921tbVV+fr5+/vOf+7//8ssv6/Dhw3rxxRf14osv+o9fdtll+uijj7odm80X5rfZNzc3y+l06m8b75ejz7kXNUQj25CrIx1C2FWPK4t0CBEx8puxue1l0vzCSIcQdp/Fdb0jWbQ51nxSg9P/VU1NTd2+jx2sM3nikyOr5HBcZOJzTmhAv9khjdWqmJYHACDKMC0PALCk03vLm5iWD+/EtKWQ3AEAlvTlXeaMjo9VTMsDABBlqNwBAJZk9dXyVkZyBwBYUiQ2sYkWJHcAgCX5ZK76jt26Pch77suXL9fw4cP9+/S63W794Q9/CFVsAADAgKAq9wEDBuixxx7TkCFD5PP59Ktf/Uq333673n77bV111VWhihEAEItMTsuLafnuue222wK+/vGPf6zly5erurqa5A4A6FEsqDPO8D33jo4OrV27Vi0tLQFvs/lHra2tam1t9X/d3Nxs9JQAAKAbgk7u7777rtxut06dOqU+ffro1Vdf1dChQ8/Zv7S0VI8++qipIAEAsYdNbIwLehObK6+8Urt379a2bdv07W9/W7NmzdJf/vKXc/ZfsGCBmpqa/K2urs5UwACA2HBm+1kzLVYFXbknJCToq1/9qqTTL57fsWOHlixZomeeeabT/omJiUpMjK23vwEAEEmmn3P3er0B99QBAOgJTMsbF1RyX7BggW655RYNHDhQx44dU1lZmTZt2qTXX389VPEBAGIUb4UzLqjk3tDQoJkzZ+rTTz+V0+nU8OHD9frrr2vChAmhig8AAAQpqOS+cuXKUMUBAEAApuWNY295AIAleX2nm5nxsYrkDgCwpA6fTR0+m6nxsSro59wBAIC1UbkDACyJaXnjSO4AAEvy+mzymphaNzP2Qse0PAAAUYbKHQBgSV6f1MG0vCEkdwCAJTEtbxzT8gAARBkqdwCAJbFa3jiSOwDAktjExjim5QEAiDJU7gAAS/J+0cyMj1UkdwCAJXllcrW8mJYHAMBSziyoM9NC6ejRo5oxY4YcDodSUlJUWFio48ePdznm1KlTKioqUt++fdWnTx9NmTJF9fX1nfb97LPPNGDAANlsNjU2NgYVG8kdAAADZsyYob1796qiokIbNmzQli1bNHfu3C7H3H///Vq/fr3Wrl2rzZs369ChQ7rzzjs77VtYWKjhw4cbio3kDgCwpDOr5c20UNm3b5/Ky8v1y1/+Urm5uRo7dqyefvpprV69WocOHep0TFNTk1auXKknnnhCN910k3JycvT8889r69atqq6uDui7fPlyNTY26sEHHzQUH8kdAGBJPpNT8r4vpuWbm5sDWmtrq+nYqqqqlJKSolGjRvmP5eXlyW63a9u2bZ2OqampUXt7u/Ly8vzHsrKyNHDgQFVVVfmP/eUvf9EPf/hDvfDCC7LbjaVpkjsAIKplZmbK6XT6W2lpqenP9Hg86t+/f8Cx+Ph4paamyuPxnHNMQkKCUlJSAo6np6f7x7S2tmr69On66U9/qoEDBxqOj9XyAABL6qm95evq6uRwOPzHExMTzzlm/vz5+slPftLl5+7bt89wTOezYMECZWdn61/+5V9MfQ7JHQBgSR0m3wp3ZqzD4QhI7l154IEHNHv27C77DB48WC6XSw0NDQHHP//8cx09elQul6vTcS6XS21tbWpsbAyo3uvr6/1jNm7cqHfffVcvv/yyJMn3xb2Ffv366Xvf+54effTRbl0HyR0AgC+kpaUpLS3tvP3cbrcaGxtVU1OjnJwcSacTs9frVW5ubqdjcnJy1KtXL1VWVmrKlCmSpAMHDqi2tlZut1uS9Nvf/lYnT570j9mxY4fuuece/elPf9I//dM/dfs6SO4AAEuy8itfs7OzNXHiRM2ZM0crVqxQe3u7iouLdffddysjI0OSdPDgQY0fP14vvPCCRo8eLafTqcLCQpWUlCg1NVUOh0P33nuv3G63xowZI0lnJfAjR474z/eP9+q7QnIHAFiS1d8K99JLL6m4uFjjx4+X3W7XlClT9NRTT/m/397ergMHDujEiRP+Y08++aS/b2trq/Lz8/Xzn/+8x2Oz+c5M6IdJc3OznE6n/lZ5nxwXn3tRQ1Tqc1GkIwg/z5FIRxARtjE3RDqEiKi+7rlIhxB2zk2TIx1CWB1vPqncQd9RU1NTt+9jB+tMntiwf7kuTu5t+HNajp3UN7K+HdJYrYrKHQBgST21oC4WkdwBAJZk5XvuVkdyBwBYktXvuVsZO9QBABBlqNwBAJbkNfnyF6blAQCwGKbljWNaHgCAKEPlDgCwJFbLG0dyBwBYEs+5G8e0PAAAUYbKHQBgSacX1JmZlu/BYC4wJHcAgCWxWt44puUBAIgyVO4AAEtiQZ1xJHcAgCV1yGRy77FILjwkdwCAJflM3nP3xXDlzj13AACiDJU7AMCSuOduHMkdAGBJJHfjmJYHACDKULkDACyJyt04kjsAwJLYoc44puUBAIgyVO4AAEtiWt44kjsAwJK8JpM70/IAACBqmErujz32mGw2m+67774eCgcAgNPOLKgz02KV4Wn5HTt26JlnntHw4cN7Mh4AACRxz90MQ5X78ePHNWPGDP3iF7/QJZdc0tMxAQDgT+5mWqwylNyLioo0adIk5eXlnbdva2urmpubAxoAAAidoKflV69erV27dmnHjh3d6l9aWqpHH3006MAAALGtw2dTh89manysCqpyr6ur07x58/TSSy8pKSmpW2MWLFigpqYmf6urqzMUKAAgtrCgzrigkntNTY0aGhp0zTXXKD4+XvHx8dq8ebOeeuopxcfHq6Oj46wxiYmJcjgcAQ0AgAvd0aNHNWPGDDkcDqWkpKiwsFDHjx/vcsypU6dUVFSkvn37qk+fPpoyZYrq6+vP6rdq1SoNHz5cSUlJ6t+/v4qKioKKLahp+fHjx+vdd98NOFZQUKCsrCx997vfVVxcXFAnBwDgXKy+Wn7GjBn69NNPVVFRofb2dhUUFGju3LkqKys755j7779fv//977V27Vo5nU4VFxfrzjvv1FtvveXv88QTT+jxxx/XT3/6U+Xm5qqlpUUfffRRULEFldyTk5M1bNiwgGMXX3yx+vbte9ZxAADMsPIOdfv27VN5ebl27NihUaNGSZKefvpp3XrrrfrZz36mjIyMs8Y0NTVp5cqVKisr00033SRJev7555Wdna3q6mqNGTNGf/vb37Rw4UKtX79e48eP948N9rFzdqgDACBIVVVVSklJ8Sd2ScrLy5Pdbte2bds6HVNTU6P29vaAJ82ysrI0cOBAVVVVSZIqKirk9Xp18OBBZWdna8CAAbrrrruCXq9mem/5TZs2mf0IAADO4vWebmbGSzrrEezExEQlJiaaiEzyeDzq379/wLH4+HilpqbK4/Gcc0xCQoJSUlICjqenp/vH/PWvf5XX69V//ud/asmSJXI6nVq4cKEmTJigd955RwkJCd2Kj8odAGBJPbWJTWZmppxOp7+Vlpae85zz58+XzWbrsu3fvz9k1+z1etXe3q6nnnpK+fn5GjNmjH7961/rvffe0xtvvNHtz+GtcACAqFZXVxfwpFZXVfsDDzyg2bNnd/l5gwcPlsvlUkNDQ8Dxzz//XEePHpXL5ep0nMvlUltbmxobGwOq9/r6ev+YSy+9VJI0dOhQ//fT0tLUr18/1dbWdhnXl5HcAQCW1CGTq+W/+DOYx7DT0tKUlpZ23n5ut1uNjY2qqalRTk6OJGnjxo3yer3Kzc3tdExOTo569eqlyspKTZkyRZJ04MAB1dbWyu12S5Kuv/56//EBAwZIOv3I3ZEjR3TZZZd16xokpuUBABZl5U1ssrOzNXHiRM2ZM0fbt2/XW2+9peLiYt19993+lfIHDx5UVlaWtm/fLklyOp0qLCxUSUmJ3njjDdXU1KigoEBut1tjxoyRJF1xxRW6/fbbNW/ePG3dulV79uzRrFmzlJWVpRtvvLHb8ZHcAQCWZPUXx7z00kvKysrS+PHjdeutt2rs2LF69tln/d9vb2/XgQMHdOLECf+xJ598Ut/4xjc0ZcoUff3rX5fL5dIrr7wS8LkvvPCCcnNzNWnSJI0bN069evVSeXm5evXq1e3YmJYHAMCA1NTULjesGTRokHy+wN8wkpKStGzZMi1btuyc4xwOh1auXKmVK1cajo3kDgCwJKvvUGdlJHcAgCVZeYc6q+OeOwAAUYbKHQBgSWZXvMdy5U5yBwBYEvfcjWNaHgCAKEPlDgCwJK/XJq/XZmp8rCK5AwAsydthk7fDRHI3MfZCx7Q8AABRhsodAGBJTMsbR3IHAFgSyd04kjsAwJJI7sZxzx0AgChD5Q4AsCSf19xqeV8MV+4kdwCAJTEtbxzT8gAARBkqdwCAJVG5Gxe55P5Zo3QyIWKnjwTb4CsiHUL4paRGOoKI+CzueKRDiAjnpsmRDiHstvX9daRDCKuTvs/Ddi6v93QzMz5WMS0PAECUYVoeAGBJ7C1vHMkdAGBJ3HM3jml5AACiDJU7AMCSqNyNI7kDACypo8Mmu4n75h3ccwcAwFq8PpOVuy92kzv33AEAiDJU7gAAS/KZvOfOi2MAALAYFtQZx7Q8AABRhsodAGBJ7FBnHMkdAGBJTMsbx7Q8AABRhsodAGBJp1/5aqZy78FgLjBU7gAASzozLW+mhdLRo0c1Y8YMORwOpaSkqLCwUMePH+9yzKlTp1RUVKS+ffuqT58+mjJliurr6wP67NixQ+PHj1dKSoouueQS5efn689//nNQsZHcAQAwYMaMGdq7d68qKiq0YcMGbdmyRXPnzu1yzP3336/169dr7dq12rx5sw4dOqQ777zT//3jx49r4sSJGjhwoLZt26Y333xTycnJys/PV3t7e7djY1oeAGBJVl4tv2/fPpWXl2vHjh0aNWqUJOnpp5/Wrbfeqp/97GfKyMg4a0xTU5NWrlypsrIy3XTTTZKk559/XtnZ2aqurtaYMWO0f/9+HT16VD/84Q+VmZkpSVq0aJGGDx+ujz/+WF/96le7FR+VOwDAks7sLW+4fbG3fHNzc0BrbW01HVtVVZVSUlL8iV2S8vLyZLfbtW3btk7H1NTUqL29XXl5ef5jWVlZGjhwoKqqqiRJV155pfr27auVK1eqra1NJ0+e1MqVK5Wdna1BgwZ1Oz6SOwDAknrqnntmZqacTqe/lZaWmo7N4/Gof//+Acfi4+OVmpoqj8dzzjEJCQlKSUkJOJ6enu4fk5ycrE2bNunFF19U79691adPH5WXl+sPf/iD4uO7P9lOcgcARLW6ujo1NTX524IFC87Zd/78+bLZbF22/fv3hyzWkydPqrCwUNdff72qq6v11ltvadiwYZo0aZJOnjzZ7c/hnjsAwJJ8Ju+5+74Y63A45HA4ujXmgQce0OzZs7vsM3jwYLlcLjU0NAQc//zzz3X06FG5XK5Ox7lcLrW1tamxsTGgeq+vr/ePKSsr00cffaSqqirZ7Xb/sUsuuUS/+93vdPfdd3frOkjuAABL8nptsoV5h7q0tDSlpaWdt5/b7VZjY6NqamqUk5MjSdq4caO8Xq9yc3M7HZOTk6NevXqpsrJSU6ZMkSQdOHBAtbW1crvdkqQTJ07IbrfLZvt77Ge+9gbx4D7T8gAABCk7O1sTJ07UnDlztH37dr311lsqLi7W3Xff7V8pf/DgQWVlZWn79u2SJKfTqcLCQpWUlOiNN95QTU2NCgoK5Ha7NWbMGEnShAkT9Le//U1FRUXat2+f9u7dq4KCAsXHx+vGG2/sdnwkdwCANXl95lsIvfTSS8rKytL48eN16623auzYsXr22Wf9329vb9eBAwd04sQJ/7Enn3xS3/jGNzRlyhR9/etfl8vl0iuvvOL/flZWltavX6933nlHbrdbX/va13To0CGVl5fr0ksv7XZsTMsDACzJ7vXJbiJB+0Kc3FNTU1VWVnbO7w8aNEg+X2AMSUlJWrZsmZYtW3bOcRMmTNCECRNMxRZU5f6DH/zgrFWDWVlZpgIAAAA9K+jK/aqrrtIf//jHv39AEM/dAQDQXbYOn2wdxqtvM2MvdEFn5vj4+HMu8wcAoKfYTE7Le0M8LW9lQS+oe++995SRkaHBgwdrxowZqq2t7bJ/a2vrWVv/AQCA0Akquefm5mrVqlUqLy/X8uXL9eGHH+prX/uajh07ds4xpaWlAdv+ndkIHwCArth9Pv+iOkPNR+XeLbfccoumTp2q4cOHKz8/X6+99poaGxv1m9/85pxjFixYELDtX11dnemgAQDRz+b1mW6xytRquJSUFF1xxRV6//33z9knMTFRiYmJZk4DAIhB9g7JbmJRnL2jB4O5wJjaxOb48eP64IMPgnqwHgAAhFZQyf3BBx/U5s2b9dFHH2nr1q264447FBcXp+nTp4cqPgBAjDJ1v93kSvsLXVDT8p988ommT5+uzz77TGlpaRo7dqyqq6u7tck+AADBMJugSe7dtHr16lDFAQAAegjbywEALMnsindWywMAYDFMyxvHK18BAIgyVO4AAEuyd/hMPuceu5U7yR0AYElntp81Mz5WMS0PAECUoXIHAFiT2f3hY3hBHckdAGBJrJY3juQOALAkFtQZxz13AACiDJU7AMCS2KHOOJI7AMCS7F6v7F6vqfGximl5AACiDJU7AMCSWC1vHMkdAGBJrJY3jml5AACiDJU7AMCSbD6Tq+VjeG95kjsAwJK4524c0/IAAEQZKncAgCVRuRtHcgcAWJK9wye7ndXyRpDcAQDW5DW5hWzsblDHPXcAAKJN2Ct33xePJjSfaAv3qSPOduxkpEMIv7YYvGZJx3r1inQIEXH8ZOz9d33S93mkQwirkzp9vb4wPGbW0XZCn5uo3Ds+j82/fyTJ5gvHP6Ev+eSTT5SZmRnOUwIAelhdXZ0GDBgQks8+deqULr/8cnk8HtOf5XK59OGHHyopKakHIrtwhD25e71eHTp0SMnJybLZbGE7b3NzszIzM1VXVyeHwxG280Ya1x071x2L1yzF5nVH8pp9Pp+OHTumjIwM2e2hu7N76tQptbWZnwlKSEiIucQuRWBa3m63h+y3ve5wOBwx8xfAl3HdsSMWr1mKzeuO1DU7nc6QnyMpKSkmk3JPYUEdAABRhuQOAECUiZnknpiYqEWLFikxMTHSoYQV1x071x2L1yzF5nXH4jUjOGFfUAcAAEIrZip3AABiBckdAIAoQ3IHACDKkNwBAIgyMZPcly1bpkGDBikpKUm5ubnavn17pEMKqS1btui2225TRkaGbDab1q1bF+mQQq60tFTXXnutkpOT1b9/f02ePFkHDhyIdFght3z5cg0fPty/oYnb7dYf/vCHSIcVVo899phsNpvuu+++SIcSUj/4wQ9ks9kCWlZWVqTDggXFRHJfs2aNSkpKtGjRIu3atUsjRoxQfn6+GhoaIh1ayLS0tGjEiBFatmxZpEMJm82bN6uoqEjV1dWqqKhQe3u7br75ZrW0tEQ6tJAaMGCAHnvsMdXU1Gjnzp266aabdPvtt2vv3r2RDi0sduzYoWeeeUbDhw+PdChhcdVVV+nTTz/1tzfffDPSIcGKfDFg9OjRvqKiIv/XHR0dvoyMDF9paWkEowofSb5XX3010mGEXUNDg0+Sb/PmzZEOJewuueQS3y9/+ctIhxFyx44d8w0ZMsRXUVHhGzdunG/evHmRDimkFi1a5BsxYkSkw8AFIOor97a2NtXU1CgvL89/zG63Ky8vT1VVVRGMDKHW1NQkSUpNTY1wJOHT0dGh1atXq6WlRW63O9LhhFxRUZEmTZoU8N93tHvvvfeUkZGhwYMHa8aMGaqtrY10SLCgsL84JtyOHDmijo4OpaenBxxPT0/X/v37IxQVQs3r9eq+++7T9ddfr2HDhkU6nJB799135Xa7derUKfXp00evvvqqhg4dGumwQmr16tXatWuXduzYEelQwiY3N1erVq3SlVdeqU8//VSPPvqovva1r2nPnj1KTk6OdHiwkKhP7ohNRUVF2rNnT8zcj7zyyiu1e/duNTU16eWXX9asWbO0efPmqE3wdXV1mjdvnioqKmLqzWG33HKL//8PHz5cubm5uuyyy/Sb3/xGhYWFEYwMVhP1yb1fv36Ki4tTfX19wPH6+nq5XK4IRYVQKi4u1oYNG7Rly5aIvl44nBISEvTVr35VkpSTk6MdO3ZoyZIleuaZZyIcWWjU1NSooaFB11xzjf9YR0eHtmzZoqVLl6q1tVVxcXERjDA8UlJSdMUVV+j999+PdCiwmKi/556QkKCcnBxVVlb6j3m9XlVWVsbEPclY4vP5VFxcrFdffVUbN27U5ZdfHumQIsbr9aq1tTXSYYTM+PHj9e6772r37t3+NmrUKM2YMUO7d++OicQuScePH9cHH3ygSy+9NNKhwGKivnKXpJKSEs2aNUujRo3S6NGjtXjxYrW0tKigoCDSoYXM8ePHA36b//DDD7V7926lpqZq4MCBEYwsdIqKilRWVqbf/e53Sk5OlsfjkSQ5nU717t07wtGFzoIFC3TLLbdo4MCBOnbsmMrKyrRp0ya9/vrrkQ4tZJKTk89aS3HxxRerb9++Ub3G4sEHH9Rtt92myy67TIcOHdKiRYsUFxen6dOnRzo0WExMJPdp06bp8OHDeuSRR+TxeDRy5EiVl5eftcgumuzcuVM33nij/+uSkhJJ0qxZs7Rq1aoIRRVay5cvlyTdcMMNAceff/55zZ49O/wBhUlDQ4NmzpypTz/9VE6nU8OHD9frr7+uCRMmRDo09LBPPvlE06dP12effaa0tDSNHTtW1dXVSktLi3RosBhe+QoAQJSJ+nvuAADEGpI7AABRhuQOAECUIbkDABBlSO4AAEQZkjsAAFGG5A4AQJQhuQMAEGVI7gAARBmSOwAAUYbkDgBAlCG5AwAQZf4/WOAzlnWI/QoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#import autotokenizer\n",
    "from transformers import AutoTokenizer,AutoModel\n",
    "import torch\n",
    "from torch import nn\n",
    "\n",
    "from model.hyena import Hyena\n",
    "\n",
    "with torch.no_grad():\n",
    "    hyena_tokenizer=AutoTokenizer.from_pretrained(\"weight/hyenadna/hyenadna-large-1m-seqlen\",trust_remote_code=True)\n",
    "    hyena_model=Hyena(pretrained_weight='weight/hyenadna/hyenadna-large-1m-seqlen',max_length=6000).to('cuda')\n",
    "    full_sequence=[]\n",
    "    checkpoint=torch.load(hyena_path)\n",
    "    # torch.nn.modules.utils.consume_prefix_in_state_dict_if_present(\n",
    "    #         checkpoint, \"model.backbone.\"\n",
    "    #     )\n",
    "    # torch.nn.modules.utils.consume_prefix_in_state_dict_if_present(\n",
    "    #         checkpoint, \"decoder.0.length_transform_2.\"\n",
    "    #     )\n",
    "\n",
    "    #edit key name in hyena_decoder\n",
    "    \n",
    "    hyena_model.load_state_dict(checkpoint)\n",
    "    hyena_model.eval()\n",
    "\n",
    "    sequence_encoded=hyena_tokenizer(sequence,\n",
    "                        add_special_tokens= False,  # this is what controls adding eos\n",
    "                        padding=\"max_length\",\n",
    "                        max_length=max_length*batch_size,\n",
    "                        truncation=True,\n",
    "                    )\n",
    "\n",
    "\n",
    "    sequence_tensor=torch.tensor(sequence_encoded['input_ids']).to('cuda')\n",
    "    sequence_tensor=torch.reshape(sequence_tensor,(batch_size,max_length))\n",
    "    mat_hyena=hyena_model(sequence_tensor)\n",
    "\n",
    "    #calculate the \n",
    "    corr=genomic_structure_corr(mat_hyena,target)\n",
    "    print(corr)\n",
    "    mse=genomic_structure_mse(mat_hyena,target)\n",
    "    print(mse)\n",
    "from matplotlib import pyplot as plt\n",
    "def figshow(x, np=False):\n",
    "        if np:\n",
    "            plt.imshow(x.squeeze())\n",
    "        else:\n",
    "            plt.imshow(x.squeeze().cpu().detach().numpy())\n",
    "        #remove boudary\n",
    "        plt.axis('off')\n",
    "        plt.show()\n",
    "#plot the target and predicted genomic structure\n",
    "# figure,(axes1, axes2, axes3, axes4, axes5)=plt.subplots(1, 5,figsize=(20,5))\n",
    "# axes1.imshow(target[0].squeeze().cpu().detach().numpy())\n",
    "# axes1.axis('off')\n",
    "# axes1\n",
    "\n",
    "plt.imshow(mat_hyena[0].cpu().detach().numpy(),cmap = 'Spectral_r')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.353693892414759\n",
      "tensor(0.1047, device='cuda:0')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f0fcf51c910>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGdCAYAAAAPGjobAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1CUlEQVR4nO3df3xU1Z3/8fdMQhKUzMRAyJgSRLZoggh8DRJGaVEJBqVWlEXkwRaIecBju4lFo27BUqntdmO3VUGhoC1KXU2hWKULtXHTIFAl4UeQKhRYtWoiOAlIk0CAJGbm+wcydUoImXszM5eZ17OP86C5OWfu50blk8+5555r8/l8PgEAgKhhj3QAAACgZ5HcAQCIMiR3AACiDMkdAIAoQ3IHACDKkNwBAIgyJHcAAKIMyR0AgCgTH+4Ter1eHTp0SMnJybLZbOE+PQDABJ/Pp2PHjikjI0N2e+jqw1OnTqmtrc305yQkJCgpKakHIrqwhD25Hzp0SJmZmeE+LQCgB9XV1WnAgAEh+exTp07J1dupJplP7i6XSx9++GHMJfiwJ/fk5GRJ0p7pY5WcEPbTR1T1g9MjHULYpfVuj3QIEfH/aj+IdAiREYu7WX/WGOkIwqr5RJsum13m/7s8FNra2tSkNj1hv169TaSpk/pcJZ631NbWRnIPtTNT8ckJ8XLEWHK/KLl3pEMIuz69Y+uf8RmOPomRDiEyvDGY3E8mRDqCiAjHbdXeildvm4m/Q2LwX8czYvNvXgCA5dnjJLuJ3yHsPkneHgvngkJyBwBYks1uk93EDIHNF7uLtknuAABLiouT4kzk57gYnpbnOXcAAKIMlTsAwJLsJqfl7UzLAwBgLXZ7Dyyoi1FMywMAEGWo3AEAlmSPY1reKJI7AMCS7LbTU/OGx8foM+4S0/IAAEQdKncAgCXZ42yym1hRZ2ZK/0JHcgcAWFKc/XQzPL7nQrngMC0PAECUoXIHAFiS3W5yWl5MywMAYCm2OHOr5WP4ljvJHQBgTXa7TFbusSuWrx0AgKhEcgcAWNLpyt1cM2LZsmUaNGiQkpKSlJubq+3bt3fZf+3atcrKylJSUpKuvvpqvfbaa2f12bdvn775zW/K6XTq4osv1rXXXqva2lpjAXYDyR0AYEl2u+30s+5Gm4Ep/TVr1qikpESLFi3Srl27NGLECOXn56uhoaHT/lu3btX06dNVWFiot99+W5MnT9bkyZO1Z88ef58PPvhAY8eOVVZWljZt2qR33nlH3//+95WUlGT4Z3M+Np/PF9b35jQ3N8vpdOrjWTfIkRBbt/zf/N7MSIcQdum92yMdQkTkfPxepEOIDG8MvobryN8iHUFYNZ9o0yV3rVJTU5McDkdozvFFnlg/eIIujutl+HNaOtp1218rgoo1NzdX1157rZYuXSpJ8nq9yszM1L333qv58+ef1X/atGlqaWnRhg0b/MfGjBmjkSNHasWKFZKku+++W7169dJ///d/G76WYFG5AwAs6cwmNmaadPqXhS+31tbWTs/X1tammpoa5eXl+Y/Z7Xbl5eWpqqqq0zFVVVUB/SUpPz/f39/r9er3v/+9rrjiCuXn56t///7Kzc3VunXrzP+AukByBwBYkqkp+S+aJGVmZsrpdPpbaWlpp+c7cuSIOjo6lJ6eHnA8PT1dHo+n0zEej6fL/g0NDTp+/Lgee+wxTZw4Uf/7v/+rO+64Q3feeac2b95s9kd0TrE1Lw4AiDl1dXUB0/KJiYlhO7fXe/rVdLfffrvuv/9+SdLIkSO1detWrVixQuPGjQvJeQ1V7sGuJAQAIFg9tVre4XAEtHMl9379+ikuLk719fUBx+vr6+VyuTod43K5uuzfr18/xcfHa+jQoQF9srOzrbVaPtiVhAAAGHFm+1kzLRgJCQnKyclRZWWl/5jX61VlZaXcbnenY9xud0B/SaqoqPD3T0hI0LXXXqsDBw4E9Pm///s/XXbZZUHFF4ygk/sTTzyhOXPmqKCgQEOHDtWKFSt00UUX6bnnngtFfAAAhE1JSYl+8Ytf6Fe/+pX27dunb3/722ppaVFBQYEkaebMmVqwYIG//7x581ReXq7HH39c+/fv1w9+8APt3LlTxcXF/j4PPfSQ1qxZo1/84hd6//33tXTpUq1fv17/9m//FrLrCOqe+5mVhF++sPOtJGxtbQ1Ymdjc3GwwVABALLHFSXYT7201snHttGnTdPjwYT3yyCPyeDwaOXKkysvL/YvmamtrZf/S7jjXXXedysrKtHDhQj388MMaMmSI1q1bp2HDhvn73HHHHVqxYoVKS0v1ne98R1deeaV++9vfauzYscYv7jyCSu5drSTcv39/p2NKS0v16KOPGo8QABCTTL8VzmdsbHFxcUDl/WWbNm0669jUqVM1derULj/znnvu0T333GMoHiNC/ijcggUL1NTU5G91dXWhPiUAIArY48y3WBVU5W5kJWFiYmJYHzsAACDWBVW5G1lJCACAEXE2m+LsJloMv9A96E1sSkpKNGvWLI0aNUqjR4/W4sWLA1YSAgDQE8y82e3M+FgVdHI/30pCAAAQWYa2n+1qJSEAAD3hy/vDGxpvcLV8NGBveQCAJTEtb1wMXzoAANGJyh0AYEk2u082u8/U+FhFcgcAWJLNfrqZGR+rYvjSAQCITlTuAABLstl8stlMTMubGHuhI7kDACyJaXnjSO4AAEuy2X2ys6DOkBj+vQYAgOhE5Q4AsCSbzeS0fOxuUEdyBwBYE8+5G8e0PAAAUYbKHQBgSayWN47kDgCwJF4cY1wMXzoAANGJyh0AYEnsUGccyR0AYEncczcuhi8dAIDoROUOALAknnM3juQOALAkpuWNI7kDACzp9KNwxqtvHoUDAABRg8odAGBJTMsbR3IHAFiSTSafcxcL6sKu+sHpuii5d6ROHxFfvzQp0iGE3b13RTqCyFi16upIhxARvoPvRzqEsLMNviLSIYSV7dhJSasiHQbOg8odAGBJTMsbF8OXDgCwsjPPuZtpRixbtkyDBg1SUlKScnNztX379i77r127VllZWUpKStLVV1+t11577Zx9//Vf/1U2m02LFy82FFt3kdwBAPjCmjVrVFJSokWLFmnXrl0aMWKE8vPz1dDQ0Gn/rVu3avr06SosLNTbb7+tyZMna/LkydqzZ89ZfV999VVVV1crIyMj1JdBcgcAWJPN/vfXvhppRqbln3jiCc2ZM0cFBQUaOnSoVqxYoYsuukjPPfdcp/2XLFmiiRMn6qGHHlJ2drZ+9KMf6ZprrtHSpUsD+h08eFD33nuvXnrpJfXq1cvIjyMoJHcAgCX11LR8c3NzQGttbe30fG1tbaqpqVFeXp7/mN1uV15enqqqqjodU1VVFdBfkvLz8wP6e71efetb39JDDz2kq666yuyPpVtI7gCAqJaZmSmn0+lvpaWlnfY7cuSIOjo6lJ6eHnA8PT1dHo+n0zEej+e8/X/yk58oPj5e3/nOd0xeSfexWh4AYEk9tVq+rq5ODofDfzwxMdFkZN1XU1OjJUuWaNeuXbLZbGE7L5U7AMCSbDaf6SZJDocjoJ0ruffr109xcXGqr68POF5fXy+Xy9XpGJfL1WX/P/3pT2poaNDAgQMVHx+v+Ph4ffzxx3rggQc0aNAgkz+hcyO5AwAsycxiujMtGAkJCcrJyVFlZaX/mNfrVWVlpdxud6dj3G53QH9Jqqio8Pf/1re+pXfeeUe7d+/2t4yMDD300EN6/fXXgwswCEzLAwDwhZKSEs2aNUujRo3S6NGjtXjxYrW0tKigoECSNHPmTH3lK1/x37efN2+exo0bp8cff1yTJk3S6tWrtXPnTj377LOSpL59+6pv374B5+jVq5dcLpeuvPLKkF0HyR0AYElmNqI5Mz5Y06ZN0+HDh/XII4/I4/Fo5MiRKi8v9y+aq62tlf1LUwLXXXedysrKtHDhQj388MMaMmSI1q1bp2HDhhmOuyeQ3AEAlhSp7WeLi4tVXFzc6fc2bdp01rGpU6dq6tSp3f78jz76yFhgQeCeOwAAUYbKHQBgTXabFGfi8TF7+B49sxqSOwDAkmx2m2wmErSZsRc6puUBAIgyVO4AAGuKs59uZsbHKJI7AMCa7DZz982ZlgcAANGCyh0AYEm2OMlmYrW8La4Hg7nAkNwBANbEtLxhJHcAgDXFmXzO3czYCxz33AEAiDJU7gAAS7LZTG5iY4vdyp3kDgCwJp5zNyx2rxwAgCgVdHLfsmWLbrvtNmVkZMhms2ndunUhCAsAEOvO7C1vpsWqoJN7S0uLRowYoWXLloUiHgAATjuzWt5Mi1FB33O/5ZZbdMstt4QiFgAA0ANCvqCutbVVra2t/q+bm5tDfUoAQDTgOXfDQr6grrS0VE6n098yMzNDfUoAQBTgnrtxIU/uCxYsUFNTk7/V1dWF+pQAAMS0kE/LJyYmKjExMdSnAQBEG6blDWMTGwCANdnskt3EBLMtdrdyCTq5Hz9+XO+//77/6w8//FC7d+9WamqqBg4c2KPBAQBily3OZvKVr1Tu3bZz507deOON/q9LSkokSbNmzdKqVat6LDAAAGBM0Mn9hhtukM/nC0UsAAD8He9zN4x77gAAa2JBnWGxu9oAAIAoReUOALAksxvRxPImNiR3AIA18T53w2L3ygEAiFJU7gAAa4qTyQV1PRbJBYfkDgCwJJvN5D13W+zec2daHgCAL1m2bJkGDRqkpKQk5ebmavv27V32X7t2rbKyspSUlKSrr75ar732mv977e3t+u53v6urr75aF198sTIyMjRz5kwdOnQopNdAcgcAWNOZ59zNtCCtWbNGJSUlWrRokXbt2qURI0YoPz9fDQ0NnfbfunWrpk+frsLCQr399tuaPHmyJk+erD179kiSTpw4oV27dun73/++du3apVdeeUUHDhzQN7/5TVM/mvOx+cK83Vxzc7OcTqfWvPuMLkruHc5TR9zXL02KdAhhd+9dpyIdQkSsWhWbN/t8B98/f6coY+s/INIhhFXzsZNK+adiNTU1yeFwhOYcX+SJo6umy3FRgvHPOdGm1Nm/DirW3NxcXXvttVq6dKkkyev1KjMzU/fee6/mz59/Vv9p06appaVFGzZs8B8bM2aMRo4cqRUrVnR6jh07dmj06NH6+OOPQ/ZOFip3AIAlnXlxjJkWjLa2NtXU1CgvL89/zG63Ky8vT1VVVZ2OqaqqCugvSfn5+efsL0lNTU2y2WxKSUkJKr5gsKAOABDVmpubA75OTExUYmLiWf2OHDmijo4OpaenBxxPT0/X/v37O/1sj8fTaX+Px9Np/1OnTum73/2upk+fHrKZD4nKHQBgVXa7+SYpMzNTTqfT30pLSyNyOe3t7brrrrvk8/m0fPnykJ6Lyh0AYE099Fa4urq6gCq5s6pdkvr166e4uDjV19cHHK+vr5fL5ep0jMvl6lb/M4n9448/1saNG0NatUtU7gCAKOdwOALauZJ7QkKCcnJyVFlZ6T/m9XpVWVkpt9vd6Ri32x3QX5IqKioC+p9J7O+9957++Mc/qm/fvj1wVV2jcgcAWNOXptYNjw9SSUmJZs2apVGjRmn06NFavHixWlpaVFBQIEmaOXOmvvKVr/in9ufNm6dx48bp8ccf16RJk7R69Wrt3LlTzz77rKTTif2f//mftWvXLm3YsEEdHR3++/GpqalKSDD+NEBXSO4AAGuKQHKfNm2aDh8+rEceeUQej0cjR45UeXm5f9FcbW2t7F/63Ouuu05lZWVauHChHn74YQ0ZMkTr1q3TsGHDJEkHDx7U//zP/0iSRo4cGXCuN954QzfccIOxazsPkjsAAF9SXFys4uLiTr+3adOms45NnTpVU6dO7bT/oEGDFObtZCSR3AEAVmUzuaAuhveWJ7kDAKwpAtPy0SJ2rxwAgChF5Q4AsCYqd8NI7gAAa+qhTWxiUcSSe1rvdvXpHVu/W9x7V6QjCL+HnmuNdAgRUT1uXaRDiIjcJwZHOoTwS0mNdATh1XYyfOey20xW7rGb3GN3zgIAgCgVW6UzAODCwT13w0juAABr4p67YbH7aw0AAFGKyh0AYE1MyxtGcgcAWBPJ3bDYvXIAAKIUlTsAwJJsNptsJl7+YmbshY7kDgCwJpvJaXlb7E5Ox+6VAwAQpajcAQDWxII6w0juAABrYhMbw0juAABronI3LHavHACAKEXlDgCwJqblDSO5AwCsife5G8a0PAAAUYbKHQBgTSyoM4zkDgCwJu65Gxa7v9YAABClqNwBANbEtLxhJHcAgDXx4hjDSO4AAGvinrthsftrDQAAUYrKHQBgTdxzN4zkDgCwJpvd3H3zGL7nHrtXDgBAlAoquZeWluraa69VcnKy+vfvr8mTJ+vAgQOhig0AEMvOVO5mmgHLli3ToEGDlJSUpNzcXG3fvr3L/mvXrlVWVpaSkpJ09dVX67XXXgv4vs/n0yOPPKJLL71UvXv3Vl5ent577z1DsXVXUFe+efNmFRUVqbq6WhUVFWpvb9fNN9+slpaWUMUHAIhVZ+65m2lBWrNmjUpKSrRo0SLt2rVLI0aMUH5+vhoaGjrtv3XrVk2fPl2FhYV6++23NXnyZE2ePFl79uzx9/mv//ovPfXUU1qxYoW2bdumiy++WPn5+Tp16pThH835BHXl5eXlmj17tq666iqNGDFCq1atUm1trWpqakIVHwAAYfPEE09ozpw5Kigo0NChQ7VixQpddNFFeu655zrtv2TJEk2cOFEPPfSQsrOz9aMf/UjXXHONli5dKul01b548WItXLhQt99+u4YPH64XXnhBhw4d0rp160J2HabuuTc1NUmSUlNTz9mntbVVzc3NAQ0AgPOy2UxOy59+zv0fc1Bra2unp2tra1NNTY3y8vL8x+x2u/Ly8lRVVdXpmKqqqoD+kpSfn+/v/+GHH8rj8QT0cTqdys3NPedn9gTDyd3r9eq+++7T9ddfr2HDhp2zX2lpqZxOp79lZmYaPSUAIJb00D33zMzMgDxUWlra6emOHDmijo4OpaenBxxPT0+Xx+PpdIzH4+my/5k/g/nMnmD4UbiioiLt2bNHb775Zpf9FixYoJKSEv/Xzc3NJHgAQNjU1dXJ4XD4v05MTIxgNOFhKLkXFxdrw4YN2rJliwYMGNBl38TExJj4QQIAelgPPefucDgCkvu59OvXT3Fxcaqvrw84Xl9fL5fL1ekYl8vVZf8zf9bX1+vSSy8N6DNy5MhuX0qwgvqp+Xw+FRcX69VXX9XGjRt1+eWXhyouAECsC/Nq+YSEBOXk5KiystJ/zOv1qrKyUm63u9Mxbrc7oL8kVVRU+PtffvnlcrlcAX2am5u1bdu2c35mTwiqci8qKlJZWZl+97vfKTk52X+/wOl0qnfv3iEJEAAQoyKwQ11JSYlmzZqlUaNGafTo0Vq8eLFaWlpUUFAgSZo5c6a+8pWv+O/bz5s3T+PGjdPjjz+uSZMmafXq1dq5c6eeffbZ0yHYbLrvvvv0H//xHxoyZIguv/xyff/731dGRoYmT55s/NrOI6jkvnz5cknSDTfcEHD8+eef1+zZs3sqJgAAImLatGk6fPiwHnnkEXk8Ho0cOVLl5eX+BXG1tbWyf2lG4LrrrlNZWZkWLlyohx9+WEOGDNG6desCFpr/+7//u1paWjR37lw1NjZq7NixKi8vV1JSUsiuw+bz+Xwh+/RONDc3y+l0auP7S9UnObaq/aVze0U6hLB76LnOHzmJdsduXBfpECIi94nBkQ4h7GxZV0Q6hLBqPnZSKdn/rqampm7dxzZ0ji/yRONfl8phIk80HzuplMHFIY3VqnhxDADAmnhxjGGxe+UAAEQpKncAgEWZrNxjuH4luQMArMngy18Cxseo2L1yAACiFJU7AMCSbDa7bLY4U+NjFckdAGBNrJY3LHavHACAKEXlDgCwJip3w0juAABrIrkbRnIHAFgTj8IZFrtXDgBAlKJyBwBYE9PyhpHcAQDWRHI3LHavHACAKEXlDgCwJip3w0juAABrYrW8YbF75QAARCkqdwCANdlsJqflbT0XywWG5A4AsCbuuRsWu1cOAECUonIHAFgTlbthJHcAgDWR3A2LWHL/f7UfyNEnMVKnj4hVq66OdAhhVz1uXaRDiIiR3/RFOoSIsI25IdIhhN1ncccjHUJYHevVK3wns33RzIyPUbH7aw0AAFGKaXkAgCX5fD75fMZnwcyMvdCR3AEAluSTVz55TY2PVUzLAwAQZajcAQCW5Pvif2bGxyqSOwDAknw+r3w+E9PyJsZe6JiWBwAgylC5AwAsiWl546jcAQCWdPpROK+JFtrkfvToUc2YMUMOh0MpKSkqLCzU8eNdb2p06tQpFRUVqW/fvurTp4+mTJmi+vp6//f//Oc/a/r06crMzFTv3r2VnZ2tJUuWBB0byR0AAANmzJihvXv3qqKiQhs2bNCWLVs0d+7cLsfcf//9Wr9+vdauXavNmzfr0KFDuvPOO/3fr6mpUf/+/fXiiy9q7969+t73vqcFCxZo6dKlQcXGtDwAwJKsPC2/b98+lZeXa8eOHRo1apQk6emnn9att96qn/3sZ8rIyDhrTFNTk1auXKmysjLddNNNkqTnn39e2dnZqq6u1pgxY3TPPfcEjBk8eLCqqqr0yiuvqLi4uNvxUbkDACzpzCY2ZpokNTc3B7TW1lbTsVVVVSklJcWf2CUpLy9Pdrtd27Zt63RMTU2N2tvblZeX5z+WlZWlgQMHqqqq6pznampqUmpqalDxkdwBAJZk7n773x+jy8zMlNPp9LfS0lLTsXk8HvXv3z/gWHx8vFJTU+XxeM45JiEhQSkpKQHH09PTzzlm69atWrNmzXmn+/8R0/IAgKhWV1cnh8Ph/zox8dxvJJ0/f75+8pOfdPl5+/bt67HYurJnzx7dfvvtWrRokW6++eagxpLcAQCW1FP33B0OR0By78oDDzyg2bNnd9ln8ODBcrlcamhoCDj++eef6+jRo3K5XJ2Oc7lcamtrU2NjY0D1Xl9ff9aYv/zlLxo/frzmzp2rhQsXdiv2LyO5AwAsKRI71KWlpSktLe28/dxutxobG1VTU6OcnBxJ0saNG+X1epWbm9vpmJycHPXq1UuVlZWaMmWKJOnAgQOqra2V2+3299u7d69uuukmzZo1Sz/+8Y+DvgaJe+4AAAQtOztbEydO1Jw5c7R9+3a99dZbKi4u1t133+1fKX/w4EFlZWVp+/btkiSn06nCwkKVlJTojTfeUE1NjQoKCuR2uzVmzBhJp6fib7zxRt18880qKSmRx+ORx+PR4cOHg4qPyh0AYElWfhROkl566SUVFxdr/PjxstvtmjJlip566in/99vb23XgwAGdOHHCf+zJJ5/0921tbVV+fr5+/vOf+7//8ssv6/Dhw3rxxRf14osv+o9fdtll+uijj7odm80X5rfZNzc3y+l06m8b75ejz7kXNUQj25CrIx1C2FWPK4t0CBEx8puxue1l0vzCSIcQdp/Fdb0jWbQ51nxSg9P/VU1NTd2+jx2sM3nikyOr5HBcZOJzTmhAv9khjdWqmJYHACDKMC0PALCk03vLm5iWD+/EtKWQ3AEAlvTlXeaMjo9VTMsDABBlqNwBAJZk9dXyVkZyBwBYUiQ2sYkWJHcAgCX5ZK76jt26Pch77suXL9fw4cP9+/S63W794Q9/CFVsAADAgKAq9wEDBuixxx7TkCFD5PP59Ktf/Uq333673n77bV111VWhihEAEItMTsuLafnuue222wK+/vGPf6zly5erurqa5A4A6FEsqDPO8D33jo4OrV27Vi0tLQFvs/lHra2tam1t9X/d3Nxs9JQAAKAbgk7u7777rtxut06dOqU+ffro1Vdf1dChQ8/Zv7S0VI8++qipIAEAsYdNbIwLehObK6+8Urt379a2bdv07W9/W7NmzdJf/vKXc/ZfsGCBmpqa/K2urs5UwACA2HBm+1kzLVYFXbknJCToq1/9qqTTL57fsWOHlixZomeeeabT/omJiUpMjK23vwEAEEmmn3P3er0B99QBAOgJTMsbF1RyX7BggW655RYNHDhQx44dU1lZmTZt2qTXX389VPEBAGIUb4UzLqjk3tDQoJkzZ+rTTz+V0+nU8OHD9frrr2vChAmhig8AAAQpqOS+cuXKUMUBAEAApuWNY295AIAleX2nm5nxsYrkDgCwpA6fTR0+m6nxsSro59wBAIC1UbkDACyJaXnjSO4AAEvy+mzymphaNzP2Qse0PAAAUYbKHQBgSV6f1MG0vCEkdwCAJTEtbxzT8gAARBkqdwCAJbFa3jiSOwDAktjExjim5QEAiDJU7gAAS/J+0cyMj1UkdwCAJXllcrW8mJYHAMBSziyoM9NC6ejRo5oxY4YcDodSUlJUWFio48ePdznm1KlTKioqUt++fdWnTx9NmTJF9fX1nfb97LPPNGDAANlsNjU2NgYVG8kdAAADZsyYob1796qiokIbNmzQli1bNHfu3C7H3H///Vq/fr3Wrl2rzZs369ChQ7rzzjs77VtYWKjhw4cbio3kDgCwpDOr5c20UNm3b5/Ky8v1y1/+Urm5uRo7dqyefvpprV69WocOHep0TFNTk1auXKknnnhCN910k3JycvT8889r69atqq6uDui7fPlyNTY26sEHHzQUH8kdAGBJPpNT8r4vpuWbm5sDWmtrq+nYqqqqlJKSolGjRvmP5eXlyW63a9u2bZ2OqampUXt7u/Ly8vzHsrKyNHDgQFVVVfmP/eUvf9EPf/hDvfDCC7LbjaVpkjsAIKplZmbK6XT6W2lpqenP9Hg86t+/f8Cx+Ph4paamyuPxnHNMQkKCUlJSAo6np6f7x7S2tmr69On66U9/qoEDBxqOj9XyAABL6qm95evq6uRwOPzHExMTzzlm/vz5+slPftLl5+7bt89wTOezYMECZWdn61/+5V9MfQ7JHQBgSR0m3wp3ZqzD4QhI7l154IEHNHv27C77DB48WC6XSw0NDQHHP//8cx09elQul6vTcS6XS21tbWpsbAyo3uvr6/1jNm7cqHfffVcvv/yyJMn3xb2Ffv366Xvf+54effTRbl0HyR0AgC+kpaUpLS3tvP3cbrcaGxtVU1OjnJwcSacTs9frVW5ubqdjcnJy1KtXL1VWVmrKlCmSpAMHDqi2tlZut1uS9Nvf/lYnT570j9mxY4fuuece/elPf9I//dM/dfs6SO4AAEuy8itfs7OzNXHiRM2ZM0crVqxQe3u7iouLdffddysjI0OSdPDgQY0fP14vvPCCRo8eLafTqcLCQpWUlCg1NVUOh0P33nuv3G63xowZI0lnJfAjR474z/eP9+q7QnIHAFiS1d8K99JLL6m4uFjjx4+X3W7XlClT9NRTT/m/397ergMHDujEiRP+Y08++aS/b2trq/Lz8/Xzn/+8x2Oz+c5M6IdJc3OznE6n/lZ5nxwXn3tRQ1Tqc1GkIwg/z5FIRxARtjE3RDqEiKi+7rlIhxB2zk2TIx1CWB1vPqncQd9RU1NTt+9jB+tMntiwf7kuTu5t+HNajp3UN7K+HdJYrYrKHQBgST21oC4WkdwBAJZk5XvuVkdyBwBYktXvuVsZO9QBABBlqNwBAJbkNfnyF6blAQCwGKbljWNaHgCAKEPlDgCwJFbLG0dyBwBYEs+5G8e0PAAAUYbKHQBgSacX1JmZlu/BYC4wJHcAgCWxWt44puUBAIgyVO4AAEtiQZ1xJHcAgCV1yGRy77FILjwkdwCAJflM3nP3xXDlzj13AACiDJU7AMCSuOduHMkdAGBJJHfjmJYHACDKULkDACyJyt04kjsAwJLYoc44puUBAIgyVO4AAEtiWt44kjsAwJK8JpM70/IAACBqmErujz32mGw2m+67774eCgcAgNPOLKgz02KV4Wn5HTt26JlnntHw4cN7Mh4AACRxz90MQ5X78ePHNWPGDP3iF7/QJZdc0tMxAQDgT+5mWqwylNyLioo0adIk5eXlnbdva2urmpubAxoAAAidoKflV69erV27dmnHjh3d6l9aWqpHH3006MAAALGtw2dTh89manysCqpyr6ur07x58/TSSy8pKSmpW2MWLFigpqYmf6urqzMUKAAgtrCgzrigkntNTY0aGhp0zTXXKD4+XvHx8dq8ebOeeuopxcfHq6Oj46wxiYmJcjgcAQ0AgAvd0aNHNWPGDDkcDqWkpKiwsFDHjx/vcsypU6dUVFSkvn37qk+fPpoyZYrq6+vP6rdq1SoNHz5cSUlJ6t+/v4qKioKKLahp+fHjx+vdd98NOFZQUKCsrCx997vfVVxcXFAnBwDgXKy+Wn7GjBn69NNPVVFRofb2dhUUFGju3LkqKys755j7779fv//977V27Vo5nU4VFxfrzjvv1FtvveXv88QTT+jxxx/XT3/6U+Xm5qqlpUUfffRRULEFldyTk5M1bNiwgGMXX3yx+vbte9ZxAADMsPIOdfv27VN5ebl27NihUaNGSZKefvpp3XrrrfrZz36mjIyMs8Y0NTVp5cqVKisr00033SRJev7555Wdna3q6mqNGTNGf/vb37Rw4UKtX79e48eP948N9rFzdqgDACBIVVVVSklJ8Sd2ScrLy5Pdbte2bds6HVNTU6P29vaAJ82ysrI0cOBAVVVVSZIqKirk9Xp18OBBZWdna8CAAbrrrruCXq9mem/5TZs2mf0IAADO4vWebmbGSzrrEezExEQlJiaaiEzyeDzq379/wLH4+HilpqbK4/Gcc0xCQoJSUlICjqenp/vH/PWvf5XX69V//ud/asmSJXI6nVq4cKEmTJigd955RwkJCd2Kj8odAGBJPbWJTWZmppxOp7+Vlpae85zz58+XzWbrsu3fvz9k1+z1etXe3q6nnnpK+fn5GjNmjH7961/rvffe0xtvvNHtz+GtcACAqFZXVxfwpFZXVfsDDzyg2bNnd/l5gwcPlsvlUkNDQ8Dxzz//XEePHpXL5ep0nMvlUltbmxobGwOq9/r6ev+YSy+9VJI0dOhQ//fT0tLUr18/1dbWdhnXl5HcAQCW1CGTq+W/+DOYx7DT0tKUlpZ23n5ut1uNjY2qqalRTk6OJGnjxo3yer3Kzc3tdExOTo569eqlyspKTZkyRZJ04MAB1dbWyu12S5Kuv/56//EBAwZIOv3I3ZEjR3TZZZd16xokpuUBABZl5U1ssrOzNXHiRM2ZM0fbt2/XW2+9peLiYt19993+lfIHDx5UVlaWtm/fLklyOp0qLCxUSUmJ3njjDdXU1KigoEBut1tjxoyRJF1xxRW6/fbbNW/ePG3dulV79uzRrFmzlJWVpRtvvLHb8ZHcAQCWZPUXx7z00kvKysrS+PHjdeutt2rs2LF69tln/d9vb2/XgQMHdOLECf+xJ598Ut/4xjc0ZcoUff3rX5fL5dIrr7wS8LkvvPCCcnNzNWnSJI0bN069evVSeXm5evXq1e3YmJYHAMCA1NTULjesGTRokHy+wN8wkpKStGzZMi1btuyc4xwOh1auXKmVK1cajo3kDgCwJKvvUGdlJHcAgCVZeYc6q+OeOwAAUYbKHQBgSWZXvMdy5U5yBwBYEvfcjWNaHgCAKEPlDgCwJK/XJq/XZmp8rCK5AwAsydthk7fDRHI3MfZCx7Q8AABRhsodAGBJTMsbR3IHAFgSyd04kjsAwJJI7sZxzx0AgChD5Q4AsCSf19xqeV8MV+4kdwCAJTEtbxzT8gAARBkqdwCAJVG5Gxe55P5Zo3QyIWKnjwTb4CsiHUL4paRGOoKI+CzueKRDiAjnpsmRDiHstvX9daRDCKuTvs/Ddi6v93QzMz5WMS0PAECUYVoeAGBJ7C1vHMkdAGBJ3HM3jml5AACiDJU7AMCSqNyNI7kDACypo8Mmu4n75h3ccwcAwFq8PpOVuy92kzv33AEAiDJU7gAAS/KZvOfOi2MAALAYFtQZx7Q8AABRhsodAGBJ7FBnHMkdAGBJTMsbx7Q8AABRhsodAGBJp1/5aqZy78FgLjBU7gAASzozLW+mhdLRo0c1Y8YMORwOpaSkqLCwUMePH+9yzKlTp1RUVKS+ffuqT58+mjJliurr6wP67NixQ+PHj1dKSoouueQS5efn689//nNQsZHcAQAwYMaMGdq7d68qKiq0YcMGbdmyRXPnzu1yzP3336/169dr7dq12rx5sw4dOqQ777zT//3jx49r4sSJGjhwoLZt26Y333xTycnJys/PV3t7e7djY1oeAGBJVl4tv2/fPpWXl2vHjh0aNWqUJOnpp5/Wrbfeqp/97GfKyMg4a0xTU5NWrlypsrIy3XTTTZKk559/XtnZ2aqurtaYMWO0f/9+HT16VD/84Q+VmZkpSVq0aJGGDx+ujz/+WF/96le7FR+VOwDAks7sLW+4fbG3fHNzc0BrbW01HVtVVZVSUlL8iV2S8vLyZLfbtW3btk7H1NTUqL29XXl5ef5jWVlZGjhwoKqqqiRJV155pfr27auVK1eqra1NJ0+e1MqVK5Wdna1BgwZ1Oz6SOwDAknrqnntmZqacTqe/lZaWmo7N4/Gof//+Acfi4+OVmpoqj8dzzjEJCQlKSUkJOJ6enu4fk5ycrE2bNunFF19U79691adPH5WXl+sPf/iD4uO7P9lOcgcARLW6ujo1NTX524IFC87Zd/78+bLZbF22/fv3hyzWkydPqrCwUNdff72qq6v11ltvadiwYZo0aZJOnjzZ7c/hnjsAwJJ8Ju+5+74Y63A45HA4ujXmgQce0OzZs7vsM3jwYLlcLjU0NAQc//zzz3X06FG5XK5Ox7lcLrW1tamxsTGgeq+vr/ePKSsr00cffaSqqirZ7Xb/sUsuuUS/+93vdPfdd3frOkjuAABL8nptsoV5h7q0tDSlpaWdt5/b7VZjY6NqamqUk5MjSdq4caO8Xq9yc3M7HZOTk6NevXqpsrJSU6ZMkSQdOHBAtbW1crvdkqQTJ07IbrfLZvt77Ge+9gbx4D7T8gAABCk7O1sTJ07UnDlztH37dr311lsqLi7W3Xff7V8pf/DgQWVlZWn79u2SJKfTqcLCQpWUlOiNN95QTU2NCgoK5Ha7NWbMGEnShAkT9Le//U1FRUXat2+f9u7dq4KCAsXHx+vGG2/sdnwkdwCANXl95lsIvfTSS8rKytL48eN16623auzYsXr22Wf9329vb9eBAwd04sQJ/7Enn3xS3/jGNzRlyhR9/etfl8vl0iuvvOL/flZWltavX6933nlHbrdbX/va13To0CGVl5fr0ksv7XZsTMsDACzJ7vXJbiJB+0Kc3FNTU1VWVnbO7w8aNEg+X2AMSUlJWrZsmZYtW3bOcRMmTNCECRNMxRZU5f6DH/zgrFWDWVlZpgIAAAA9K+jK/aqrrtIf//jHv39AEM/dAQDQXbYOn2wdxqtvM2MvdEFn5vj4+HMu8wcAoKfYTE7Le0M8LW9lQS+oe++995SRkaHBgwdrxowZqq2t7bJ/a2vrWVv/AQCA0Akquefm5mrVqlUqLy/X8uXL9eGHH+prX/uajh07ds4xpaWlAdv+ndkIHwCArth9Pv+iOkPNR+XeLbfccoumTp2q4cOHKz8/X6+99poaGxv1m9/85pxjFixYELDtX11dnemgAQDRz+b1mW6xytRquJSUFF1xxRV6//33z9knMTFRiYmJZk4DAIhB9g7JbmJRnL2jB4O5wJjaxOb48eP64IMPgnqwHgAAhFZQyf3BBx/U5s2b9dFHH2nr1q264447FBcXp+nTp4cqPgBAjDJ1v93kSvsLXVDT8p988ommT5+uzz77TGlpaRo7dqyqq6u7tck+AADBMJugSe7dtHr16lDFAQAAegjbywEALMnsindWywMAYDFMyxvHK18BAIgyVO4AAEuyd/hMPuceu5U7yR0AYElntp81Mz5WMS0PAECUoXIHAFiT2f3hY3hBHckdAGBJrJY3juQOALAkFtQZxz13AACiDJU7AMCS2KHOOJI7AMCS7F6v7F6vqfGximl5AACiDJU7AMCSWC1vHMkdAGBJrJY3jml5AACiDJU7AMCSbD6Tq+VjeG95kjsAwJK4524c0/IAAEQZKncAgCVRuRtHcgcAWJK9wye7ndXyRpDcAQDW5DW5hWzsblDHPXcAAKJN2Ct33xePJjSfaAv3qSPOduxkpEMIv7YYvGZJx3r1inQIEXH8ZOz9d33S93mkQwirkzp9vb4wPGbW0XZCn5uo3Ds+j82/fyTJ5gvHP6Ev+eSTT5SZmRnOUwIAelhdXZ0GDBgQks8+deqULr/8cnk8HtOf5XK59OGHHyopKakHIrtwhD25e71eHTp0SMnJybLZbGE7b3NzszIzM1VXVyeHwxG280Ya1x071x2L1yzF5nVH8pp9Pp+OHTumjIwM2e2hu7N76tQptbWZnwlKSEiIucQuRWBa3m63h+y3ve5wOBwx8xfAl3HdsSMWr1mKzeuO1DU7nc6QnyMpKSkmk3JPYUEdAABRhuQOAECUiZnknpiYqEWLFikxMTHSoYQV1x071x2L1yzF5nXH4jUjOGFfUAcAAEIrZip3AABiBckdAIAoQ3IHACDKkNwBAIgyMZPcly1bpkGDBikpKUm5ubnavn17pEMKqS1btui2225TRkaGbDab1q1bF+mQQq60tFTXXnutkpOT1b9/f02ePFkHDhyIdFght3z5cg0fPty/oYnb7dYf/vCHSIcVVo899phsNpvuu+++SIcSUj/4wQ9ks9kCWlZWVqTDggXFRHJfs2aNSkpKtGjRIu3atUsjRoxQfn6+GhoaIh1ayLS0tGjEiBFatmxZpEMJm82bN6uoqEjV1dWqqKhQe3u7br75ZrW0tEQ6tJAaMGCAHnvsMdXU1Gjnzp266aabdPvtt2vv3r2RDi0sduzYoWeeeUbDhw+PdChhcdVVV+nTTz/1tzfffDPSIcGKfDFg9OjRvqKiIv/XHR0dvoyMDF9paWkEowofSb5XX3010mGEXUNDg0+Sb/PmzZEOJewuueQS3y9/+ctIhxFyx44d8w0ZMsRXUVHhGzdunG/evHmRDimkFi1a5BsxYkSkw8AFIOor97a2NtXU1CgvL89/zG63Ky8vT1VVVRGMDKHW1NQkSUpNTY1wJOHT0dGh1atXq6WlRW63O9LhhFxRUZEmTZoU8N93tHvvvfeUkZGhwYMHa8aMGaqtrY10SLCgsL84JtyOHDmijo4OpaenBxxPT0/X/v37IxQVQs3r9eq+++7T9ddfr2HDhkU6nJB799135Xa7derUKfXp00evvvqqhg4dGumwQmr16tXatWuXduzYEelQwiY3N1erVq3SlVdeqU8//VSPPvqovva1r2nPnj1KTk6OdHiwkKhP7ohNRUVF2rNnT8zcj7zyyiu1e/duNTU16eWXX9asWbO0efPmqE3wdXV1mjdvnioqKmLqzWG33HKL//8PHz5cubm5uuyyy/Sb3/xGhYWFEYwMVhP1yb1fv36Ki4tTfX19wPH6+nq5XK4IRYVQKi4u1oYNG7Rly5aIvl44nBISEvTVr35VkpSTk6MdO3ZoyZIleuaZZyIcWWjU1NSooaFB11xzjf9YR0eHtmzZoqVLl6q1tVVxcXERjDA8UlJSdMUVV+j999+PdCiwmKi/556QkKCcnBxVVlb6j3m9XlVWVsbEPclY4vP5VFxcrFdffVUbN27U5ZdfHumQIsbr9aq1tTXSYYTM+PHj9e6772r37t3+NmrUKM2YMUO7d++OicQuScePH9cHH3ygSy+9NNKhwGKivnKXpJKSEs2aNUujRo3S6NGjtXjxYrW0tKigoCDSoYXM8ePHA36b//DDD7V7926lpqZq4MCBEYwsdIqKilRWVqbf/e53Sk5OlsfjkSQ5nU717t07wtGFzoIFC3TLLbdo4MCBOnbsmMrKyrRp0ya9/vrrkQ4tZJKTk89aS3HxxRerb9++Ub3G4sEHH9Rtt92myy67TIcOHdKiRYsUFxen6dOnRzo0WExMJPdp06bp8OHDeuSRR+TxeDRy5EiVl5eftcgumuzcuVM33nij/+uSkhJJ0qxZs7Rq1aoIRRVay5cvlyTdcMMNAceff/55zZ49O/wBhUlDQ4NmzpypTz/9VE6nU8OHD9frr7+uCRMmRDo09LBPPvlE06dP12effaa0tDSNHTtW1dXVSktLi3RosBhe+QoAQJSJ+nvuAADEGpI7AABRhuQOAECUIbkDABBlSO4AAEQZkjsAAFGG5A4AQJQhuQMAEGVI7gAARBmSOwAAUYbkDgBAlCG5AwAQZf4/WOAzlnWI/QoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#import autotokenizer\n",
    "from transformers import AutoTokenizer,AutoModel\n",
    "import torch\n",
    "from torch import nn\n",
    "from model.mamba import Mamba\n",
    "\n",
    "with torch.no_grad():\n",
    "    Mamba_tokenizer=AutoTokenizer.from_pretrained(\"weight/mamba/caduceus-ph_seqlen-131k_d_model-256_n_layer-16\",trust_remote_code=True)\n",
    "    Mamba_model=Mamba(pretrained_weight='weight/mamba/caduceus-ph_seqlen-131k_d_model-256_n_layer-16',max_length=6000).to('cuda')\n",
    "    full_sequence=[]\n",
    "    checkpoint=torch.load(mamba_path)\n",
    "    # torch.nn.modules.utils.consume_prefix_in_state_dict_if_present(\n",
    "    #         checkpoint, \"model.backbone.\"\n",
    "    #     )\n",
    "    # torch.nn.modules.utils.consume_prefix_in_state_dict_if_present(\n",
    "    #         checkpoint, \"decoder.0.length_transform_2.\"\n",
    "    #     )\n",
    "\n",
    "    #edit key name in hyena_decoder\n",
    "    \n",
    "    Mamba_model.load_state_dict(checkpoint)\n",
    "    Mamba_model.eval()\n",
    "\n",
    "    sequence_encoded=Mamba_tokenizer(sequence,\n",
    "                        add_special_tokens= False,  # this is what controls adding eos\n",
    "                        padding=\"max_length\",\n",
    "                        max_length=max_length*batch_size,\n",
    "                        truncation=True,\n",
    "                    )\n",
    "\n",
    "\n",
    "    sequence_tensor=torch.tensor(sequence_encoded['input_ids']).to('cuda')\n",
    "    sequence_tensor=torch.reshape(sequence_tensor,(batch_size,max_length))\n",
    "    mat_mamba=Mamba_model(sequence_tensor)\n",
    "\n",
    "    #calculate the \n",
    "    corr=genomic_structure_corr(mat_mamba,target)\n",
    "    print(corr)\n",
    "    mse=genomic_structure_mse(mat_mamba,target)\n",
    "    print(mse)\n",
    "from matplotlib import pyplot as plt\n",
    "def figshow(x, np=False):\n",
    "        if np:\n",
    "            plt.imshow(x.squeeze())\n",
    "        else:\n",
    "            plt.imshow(x.squeeze().cpu().detach().numpy())\n",
    "        #remove boudary\n",
    "        plt.axis('off')\n",
    "        plt.show()\n",
    "#plot the target and predicted genomic structure\n",
    "# figure,(axes1, axes2, axes3, axes4, axes5)=plt.subplots(1, 5,figsize=(20,5))\n",
    "# axes1.imshow(target[0].squeeze().cpu().detach().numpy())\n",
    "# axes1.axis('off')\n",
    "# axes1\n",
    "\n",
    "plt.imshow(mat_hyena[0].cpu().detach().numpy(),cmap = 'Spectral_r')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAGnCAYAAABioYn3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYsUlEQVR4nO3dd3gUVdsG8Hs3ZTdtUyCF0EIChA5KidSARLqIKIqEKiC+r4gUUQSkqoAvIKAigjQVRFHpgnQp8tGRIr1IDSSE9L473x+YMUsmyZzsbBLi/buuua5k5pkzZze7m2fPOXOOTpIkCURERERkN/rirgARERFRaceEi4iIiMjOmHARERER2RkTLiIiIiI7Y8JFREREZGdMuIiIiIjsjAkXERERkZ0x4SIiIiKyMyZcRERERHbGhIuohOvfvz90Oh10Oh0mTZpU3NUhAgBcu3ZNfl3qdDqrYzn3X7t2zS7XL8nvi6+++gr169eHh4eHXMfdu3cXd7WomDHhohJn69at6Nu3L6pXrw4PDw8YjUZUrFgRjRs3xogRI7Bjx47irmKpcuvWLQwdOhTVq1eHi4sLXFxcUL58eTRu3BgDBw7Exo0breKvXbuGSZMmYdKkSZgzZ07xVLoAJ06ckOu4bNmy4q5OoQQFBVklLtmbh4cHnnjiCUyePBmJiYnFXU27mTNnjvw3tFfSZg8///wzBg8ejJMnTyIpKalY65LzNaSUlE6aNEk+HhQUlOe5eW05k8icCXBeW0lLjIuaY3FXgChbdHQ0evfuja1bt+Y6dvPmTdy8eRNHjhzBnDlzkJqaCqPRWAy1LHrjxo3DoEGDAACVKlXStOwrV67gqaeeQnR0tNX+27dv4/bt2zhy5AjS09PRpUsX+di1a9cwefJkAEDlypUxfPhwTeukhRMnTsh1DA8PR//+/Yu3QhpKSkrCiRMncOLECaxcuRL79+9H2bJli7taVvbu3Sv/XK5cuUKVMWfOHPz1118AgNatW+dKCOz5vrDFmjVr5J+7du2KESNGwNHREXXr1i3GWlFJwISLSoS0tDR07NgRR48elfe9+OKLeP7551GuXDkkJyfjzJkz2LhxI/bt21eMNS161apVQ7Vq1exS9pQpU+Rkq169enj77bdRoUIFJCYm4sSJE1izZk2u7iJbWCwWpKenw8XFRbMyS4KkpCS4u7vbrfwBAwbg1VdfRVpaGjZv3ozZs2cDAC5cuIAPP/wQn3zySYFlJCcnw83NzW51zKlFixZ2v4Y93xe2uHXrlvxzt27d0Lp1a7tcJzMzE5IkwdnZ2S7l55T9+ntUXklkgwYN8Omnn+baX5IS42IhEZUA06dPlwDI21dffZVn7PHjx6WsrCyrfffv35fGjx8v1atXT3Jzc5OMRqNUq1YtaeLEiVJiYqJV7NKlS+XrhIeHS6dPn5a6du0qmUwmydXVVerYsaN08eLFXNeNioqS3n77balmzZqSi4uLZDQapdDQUGn48OHSrVu38r3G//3f/0ktW7aUXFxcpPLly0sTJ06UsrKypDt37ki9e/eWvL29JTc3N6lTp07SpUuXrMrq16+fXNbEiROtjqWlpUlz586VmjdvLnl5eUlOTk5SuXLlpM6dO0u///57gc97zZo15bLXrVunGJOQkCD/XLlyZau/06Pb1atXJUmSrPb98ccf0rBhw6TAwEBJr9dLa9askXbt2iUfr1y5stX1Jk6cKB/r16+f1TGz2SwtW7ZMatu2rVSmTBnJyclJ8vPzk55++mlpw4YNua6ttEmSVKjrP/o3PXTokBQRESF5eHhIXl5eclxycrI0Y8YMqXHjxpKHh4fk7OwsVa1aVRoxYoR07969Av8mSs/1o3/3Vq1aycfq16+veM6vv/4qTZgwQapSpYrk4OAgffLJJ3LcwYMHpZ49e0oVKlSQnJycJC8vL6lt27Z5vgZ++OEHqX79+pLBYJAqVKggjR07Vjp79myu5zWb0msi26VLl6T//ve/UmhoqOTi4iK5urpKoaGh0uDBg6W0tDSr519py34u8ntfJCQkSFOmTJGeeOIJyd3dXXJ2dpaqVKkiDRo0SLpw4YJV7KOvhb/++kvq3bu35OPjIxmNRqlFixbS4cOHC/x75Xx95PW6kyRJSk9Pl+bMmSM99dRTkslkkpycnKQKFSpIr7zyinTkyBGrMq9evWpVxu3bt6V+/fpJvr6+kk6nk44fP55vnfJ7DUmS9Wv90fdBQec+KuffIzw8vMD4fyMmXFQi1KhRo9Bv1osXL0oVKlTI84OuTp060v379+X4nB+M5cqVk9zc3HKdU6tWLclsNsvn/Pnnn5Kfn1+e1yhbtqz0xx9/KF6jQoUKkqura65zhgwZIgUHBxd47bz+sdy/f1964okn8qxTzn+weXnqqafk+LCwMGn9+vVSbGxsnvGFSbiqVatm9XthE660tDSpXbt2eV77rbfeynVteyVc5cuXl1xcXOTfPT09JUmSpOjoaKlOnTp5Xrt8+fLSlStXCvy7PPpcP/oPr2vXrvKxmjVrKp7z6POe/Xr4/PPPJb1en2cd33vvPatrLV68WDHu0ddeTkqvCUmSpI0bNyq+F7K3Bw8e2Jxw3blzJ9djz7m5urpK27Ztk+NzvhZMJpPi+7xs2bJWXzyUqEm4kpKSrN5zj26Ojo7S8uXL5TIfTbgefVxMuB4vHDRPxS4lJQXnzp2Tf+/QoYPV8ePHj2Pfvn1WW/bYDgDo3bs3bt68CQBo06YN1qxZgw0bNiA8PBwAcPr06TzHGd25cwehoaH46aefMGfOHDg5OQEA/vzzT2zbts3qGvfu3QPwsCvju+++w+rVq1GrVi0AQExMDCIjI2GxWHJd4+bNm2jatCk2bNiA//znP/L+L7/8EnFxcViyZAm++eYbuZvt0WvnZejQoTh+/DgAwNnZGaNHj8amTZuwatUqDBw4EAaDocAyOnfuLP988OBBdO3aFT4+PggJCcGrr76a686qH3/8EfPmzZN/DwgIwN69e+VNabzO5cuXMWrUKGzevBnffPMNQkJCCqyXksmTJ8vj+3Q6HV577TWsX78eP/30E4YPHw4PDw8AD8cPjR07Vj6vQYMGVnXUwq1bt+Dj44NFixZh69atmDJlCgDgjTfewOnTp+Xrfvfdd9i8eTNeeOEF+bx+/foV+rrp6elYt24dtmzZIu974oknFGMvXryIAQMGYOPGjfjhhx/QsGFDnDlzBm+++SYsFgv0ej3GjRuHrVu34ssvv4S3tzcAYNq0adi5cycAIDExESNGjJDLbNGiBdatW4cvv/wSV69eFap7dHQ0evXqhZSUFABAcHAwvvzyS/z666/44osvEBYWBgB49dVXsXfvXgQEBMjnzps3T/77KXVt5fTf//4XFy9eBAD4+/tjyZIlWLt2rdzNmZKSgsjISCQnJ+c6NyEhAW5ubli5ciWWLl0KT09PAA/f3ytXrsz3up06dcLevXvRoEEDed/YsWOtXnfvv/8+/u///g8A4O7ujrlz52Ljxo3o1q0bACArKwuvvfYabty4oXiN69evY8qUKfj111+xcOFCobF7kydPzjWIPXucY2HO9fLyyjP+t99+Uxw0f+LECdX1LZWKO+MjunnzptW3tkWLFlkdr127dq5vgu+++64kSZJ06tQpeZ+Tk5P066+/Snv37pX27t0r/fjjj1bHsrsWc34TdXJykm7evClfq0OHDvKxefPmSZIkSX/88YfVtY8ePSrHnz592urYoUOHcl3DaDTKLWz37t2zip8/f75cVufOnXNdW5KUv8nHxcVJjo6O8v65c+cW6rlPT0+Xunfvnu838zfffNPqnPxah7LlPH/kyJG5jou2MFksFsnX11feP2LEiHwf16Pdf7Ze/9EydTqdVYumJEnSgwcPJAcHBzlm5cqV8mtx165dkpOTk3zs3Llz+dZfkgpuTQQgubm5SadOnVI8p3v37rnKHDVqlHw8IiJCrt/evXulV199VT7Ws2dPSZIkq/eQwWCw6hKdP39+rhacbDn3Z7dwffbZZ/I+d3d3q/ddQY9/165duY4rvS9iY2OtWu9++uknOT46OtqqVfKHH36QJMn6tZDzPSxJkvT666/n+zpWEh4eLp+zdOlSeb/FYpHKlCkjH5s1a5Z8LD09XQoMDJSPffzxx5Ik5W7hyvm5oIaa11Be74OCzs1u1c2W8++R11ZQi1xpx0HzVOwe/aZ0//591ef++eef8s+ZmZlo3769YlxmZibOnz+Phg0bWu2vUaMGypcvL/9epkwZ+efY2FgAsGp9c3FxwZNPPin/Xrt2bXh5eSEuLk6Obdy4ca5r+Pj45CofAJo2bSr/nPPbava183LhwgVkZWXJv3fv3j3f+Lw4Ozvjp59+wtGjR/Hzzz9j//79OHz4sNwKAQCffvopevXqhaeeeqpQ18hu3bFFTEyM1Z2UhX28WqlatSrq1atnte/ChQswm83y77169crz/NOnTyM0NLTQ19fpdGjdujVmzZqFOnXqKMYoPe853y/bt2/H9u3b86wfAFy6dEneFxISAl9fX/n35s2bC9U557XDwsKs3ndauXjxolUrc87B+2XLlkVoaKjcypLzfZ3Nw8PD6v2r9HlQWNHR0VafbTnr5uzsjCZNmmDt2rV51g2w7b2kNPB9yZIlWLp0aaHOdXTMO33Ia9B8SbzJoSgx4aJi5+bmhtDQUJw/fx4AsG3bNrz77rvy8ewP/9atW+O3334r9HWU5sTJToSy5fwQkSSp0NfKKbtbAgD0eute/Lya5bW6tloNGzaUk9H09HSsWbMG/fr1Q0ZGBoCH3Y2FTbiUuhlz3vmYM3EEkGuKCnuw9fqFneogm+j8TNn/8HQ6Hdzc3BAcHAyTyZTvObbUsbjnjyouRfF5YAtb/qaVKlXKdfdoXgm3mnPz4+npWSR3qj5uOIaLSoSc8yTt2LEDq1atUnVezZo15Z9dXFwQFxcH6eHNIFZbUlKSPKZLVI0aNeSfU1NT5XFTwMNv7dmtW4/G2lP16tXh4OAg/55z7p9sav5BbNmyxao1CwAMBgN69uyJ4OBgeV/OVoOcSaPSmLVHKU0rkT1eCHjYepWeng7gYfKTc3xStrJly1q1rhT0eAuqo+j1H6X0mB79m5w/fz7P16LoOK7sf3jNmzdHgwYNCky28qpjzvfLK6+8olg/SZLkLzk5x9tdvnwZMTEx8u+///670GPIHu8IPEzgb9++nW+86OsMeNiCkvO8/fv3yz/fv39f/lIHFN17NZuvr69Vi1nOumVmZuLw4cMF1k3LKVqo6LGFi0qEt956C6tWrcIff/wBAIiMjMTmzZvx7LPPomzZsoiJicH169dznVe3bl00btwYhw8fRmpqKp5++mkMGzYMFStWRHR0NK5evYqdO3fCYrGo/jb3qHr16uHJJ5/EsWPHADz8RzV58mQ4ODhYDTqtU6dOri5Le/H09ESPHj3kxHT06NG4desWwsPDkZSUhB07dqB+/fpWg/SVfPDBBzh9+jS6deuGli1bonLlysjKysKmTZusujVydn3m/Kdx+/ZtfP311wgODoaLi4vqxx8cHAxHR0dkZWUhPT0dPXr0QIcOHbB69WpcuXIlV7xOp8PAgQMxffp0AA8nxUxJSUHnzp2RlZWF/fv3w2g0YurUqbnqePLkSfz888/w8/ODl5cX6tSpI3x9Nby8vNC9e3esXr0awMNB1KNHj0bVqlURFxeHv/76C3v27MG5c+fy7DKyt/79+2POnDmwWCz47rvv4OHhgS5dusBgMODmzZv4888/sX79eowdOxb9+/dH+/bt4eHhgcTERKSnp+OFF17AqFGjcPfuXYwbN07o2i+99BLee+89JCYmyl+A3nnnHQQFBeHatWtYunQpfvnlF7nVt0yZMvLA/OXLl0Ov18PR0RH16tXLM+H09vbGc889Jyfkb7zxBuLj4+Hj44NZs2YhNTUVwMPkp1OnToV8FgtHp9Ohb9++8pxpEydOhJOTE4KDg7F48WJ5/q7sLzxUChXpiDGifERFRVkNOM1ve//99+XzLly4kO+0EHhk4HR+A6rzutX8zJkz+U4LUaZMmTynhXj0GjnPy3nLfF7Xzmt/TEyMVK9evTzrpGZaiObNmxf4XPfq1cvqnKysLMXnOyQkpMDHmNPgwYNzlaHT6axuksg5aD01NVVq27ZtnvXMnhZCkh4OnlaafqBt27aFvn5BA/El6eFNEflNCwGFwcl5Eb0t/9FzlAaaS9LDwev5TQsBWA/2XrRokWJMaGio1e855fX3X7dundXA9Ue3Bw8eyLHvvfeeYszevXslSbJtWoitW7fK8YWdEy4veQ2alyTbp4UQxWkhShZ2KVKJ4e/vj127dmHdunV46aWXEBQUBKPRCCcnJ/j6+qJp06YYPnw4tm/fbtWyVK1aNZw8eRITJkzAE088AXd3dxgMBlSqVAmtWrXChx9+iAULFthUt1q1auHkyZMYNWoUatSoAaPRCKPRiOrVq+Ott97CyZMncw2itrcyZcrg4MGDmD17Npo2bQpPT084OTmhXLly6NSpk3ybfX7mz5+P6dOno2PHjqhevTq8vLzg4OCAMmXKoE2bNli8eDG++eYbq3McHBywZs0atGrVCq6uroWu/+zZszF48GCUKVMGRqMRTz31FDZt2oQXX3xRMd5oNGLr1q1YvHgx2rRpAx8fHzg6OsLX1xdPP/00IiIi5Fhvb2/8/PPPaNSoUZ7TY4heXw1fX18cOnQIM2fOxFNPPSX/TQIDA/HUU09h3Lhx+OmnnwpdvhbeeOMN/N///R8iIyNRqVIlODs7w2QyITQ0FD169MDXX39tdVPCoEGDsGrVKtSrVw/Ozs4oV64c3nrrLfz444/C1+7atStOnDiB119/HdWqVYPRaISrqyuqV6+OQYMGWa1AMH78eAwZMgR+fn5CXWkBAQE4cuQIJk+ejAYNGsDV1RXOzs4ICgrCwIEDcfz4cTzzzDPCddeCm5sbfvvtN3zyyScICwuDh4cHHB0dERgYiJ49e+LAgQPo27dvsdSN7E8nSSVgJCARERFRKcYWLiIiIiI7Y8JFREREZGdMuIiIiIjsjAkXERERkZ0x4SIiIiKyMyZcRERERHbGhIuIiIjIzphwEREREdkZEy4iIiIiO2PCRURERGRnTLiIiIiI7IwJFxEREZGdMeEiIiIisjMmXERERER2xoSLiIiIyM6YcBERERHZGRMuIiIiIjtzLO4KEBERUemQlpaGjIwMTcpydnaG0WjUpKySgAkXERER2SwtLQ1VKrsj6p5Zk/ICAgJw9erVUpN0MeEiIiIim2VkZCDqnhl/HQ2CycO2EUsJiRZUbngNGRkZTLiIiIiIHuXuoYO7h86mMiyw7fySiAkXERERacYsWWCWbC+jtOFdikRERER2xhYuIiIi0owFEiywrYnL1vNLIiZcj7GZf7ZXHVvZOUao7DTJSXXsuugGQmU/SHdVHXsn3iRU9p/dJgnFZ/vgdBfVsU76LKGyPR1SVcfOO9tGqGyH3zzVBwsOiTg1e4TYCQCqrPxIdWwl/1ihsu8nq3/dNA64IVR2qln96/18rK9Q2cc7fygUDwBLLjRXHXs3U+A1AECvU99VcyA2WKjsLMlBKF7ExpafCp/z+tE+qmMPRlUSKtttsZfq2Ji6Yv9q0/zU3+XndkPsOT8zXfx9LcoCC2ztELS9hJKHXYpEREREdsYWLiIiItKMWZJglmzrErT1/JKICRcRERFphmO4lLFLkYiIiMjO2MJFREREmrFAgpktXLkw4SIiIiLNsEtRGRMuIiIi0gwHzSvjGC4iIiIiO2MLFxEREWnG8vdmaxmlDRMuIiIi0oxZg0Hztp5fEjHheoyJLNdT13BHqOzDaeqXuWjoeV2o7IMPglTHxjm7CJVdWCLL9TzncVKo7F+SaquO9TclCpWddttddWymm/1HEIgs1+PmlCFUdlSm+mWeqriKLWW1+XYt1bHpmeqXASoskeV6ahpvCZXd1S1FdezQ1LJCZV9L9lEde/62v1DZhSGyXI/jD2WEyq4w+pz64CQvobJvRauP110rms9Ish0TLiIiItKMWXq42VpGacOEi4iIiDTDMVzKeJciERERkZ2xhYuIiIg0Y4EOZuhsLqO0YcJFREREmrFIDzdbyyhtmHARERGRZswatHDZen5JxDFcRERERHbGFi4iIiLSDFu4lDHhIiIiIs1YJB0sko2D5m08vyRilyIRERE99j7//HMEBQXBaDQiLCwMhw4dUnXeqlWroNPp0K1bN7vWjwkXERERaSa7S9HWTcT333+PkSNHYuLEiTh27Bjq16+P9u3b4969e/med+3aNbz99tto2bKlLQ9ZFXYpPsbSJPVruomsjQgAQU7RqmPrGsTWcqvvon7txYQAo1DZwIeC8Q95OqSqjhVZGxEA3vS6ojr2Vrq3UNl/HlX/ncniKvpciruf7Ko6VmRtRADw91S/zmSm5CBUtqtTpurY2/Fif6PC0OvUz7MtsjYiALx7t4HqWNG5kJ7yuao6tpLbA6GyC8NtsZfqWKG1EQGUM8arjg1yvS9Utms59euMLtc/JVR2UTBDD7ON7TlmwfjZs2dj8ODBGDBgAABgwYIF2LRpE5YsWYIxY8YoX8NsRmRkJCZPnoy9e/ciLi7OpjoXhC1cRERE9NjKyMjA0aNHERERIe/T6/WIiIjAgQMH8jxvypQp8PPzw8CBA4uimmzhIiIiIu1IGgyal/4+PyEhwWq/wWCAwWCw2hcTEwOz2Qx/f3+r/f7+/jh3Trnlct++fVi8eDFOnDhhUz1FsIWLiIiINKPlGK6KFSvC09NT3qZNm2Zz/RITE9GnTx8sWrQIZcuWtbk8tdjCZUdXr17FrFmzsHXrVty8eRMAEBQUhDZt2mDIkCEYNmwYfvvttwLLmThxIiZNmmTn2hIREZUsN27cgMn0z3jPR1u3AKBs2bJwcHDA3bt3rfbfvXsXAQEBueIvX76Ma9eu4dlnn5X3WSwPx006Ojri/PnzCAkJ0eohyJhw2cnGjRvx8ssvw9HREZGRkahfvz70ej3OnTuHn3/+GV988QWWLl2KQYMGyeccPnwY8+bNw9ixY1GzZk15f7169YrjIRAREQkzS3qYJRsHzf+9lqLJZLJKuJQ4OzujYcOG2LFjhzy1g8ViwY4dOzB06NBc8TVq1MCpU6es9o0fPx6JiYmYO3cuKlasaFPd88KEyw4uX76Mnj17onLlytixYwfKlStndXzGjBmYP38+nn76aas/rNFoxLx58/DMM8+gdevWRVxrIiIi21mgg8XGEUsWiK1ePXLkSPTr1w+NGjVCkyZNMGfOHCQnJ8t3Lfbt2xfly5fHtGnTYDQaUadOHavzvby8ACDXfi0x4bKDjz/+GMnJyVi6dGmuZAt42GQ5bNiwYqgZERGRfRXH0j4vv/wyoqOjMWHCBERFRaFBgwbYsmWLPJD++vXr0OuLd9g6Ey472LhxI6pWrYqwsLDirgoREdG/wtChQxW7EAFg9+7d+Z67bNky7Sv0CCZcGktISMDt27cVlwiIi4tDVlaW/LubmxtcXFyKsHZERET2pc0YLrEuxccBp4XQWPacIe7u7rmOtW7dGr6+vvL2+eefF3X1iIiI7OrhGC7bt9KGLVwa8/DwAAAkJSXlOvbll18iMTERd+/eRe/evW2+1rroBqpjG3qqX04HEFuuJ1jwVZQpqV+iJSrLU6zwQpp3to3qWH+T+voDYsv1/LS/iVDZ0jj1C2BUH3hEqOzCaBxwQ3VsFdcYobJFlus5+kBsKaukDGfVsS4e6UJlF8aB2GDVsUNTxeYREvlH5uss9lp/wXRMdezOlFChsgsjpq7Ah1OSl1DZIsv1dPM6KlT297Hqh6OYM8SWsaLiw4RLY56enihXrhxOnz6d61j2mK5r164Vca2IiIiKhkWDtRRF71J8HLBL0Q46d+6MS5cu4dChQ8VdFSIioiKVPYbL1q20KX2PqAR455134OrqildffTXXzLcAIJXCwYBERESUN3Yp2kG1atWwcuVKvPLKKwgNDZVnmpckCVevXsXKlSuh1+tRoUKF4q4qERGRpizQF/nEp48DJlx28txzz+HUqVPyWopLliyBTqdD5cqV0blzZ7z++uuoX79+cVeTiIhIU2ZJB7Nk48SnNp5fEjHhsqOQkBDMnz9fdfyLL77I7kYiIqJSiAkXERERacaswV2KZnYpEhEREeXNIulhsfEuQ0sp7O1hwkVERESaYQuXMk4LQURERGRnbOEiIiIizVhg+12GFm2qUqIw4XqMPUh3VR178EGQUNn1XdSvvSiyNiIA+DpkqI71ckgRKruwHH5Tv2Zj2u3cC5Pn58+j6huSRdZGBAC/35xUxyZvUb8+X2GlmtXXZ/PtWkJluzplqo4VWRsRAO7FmFTHtq9xVqjswsgSWDfyWrKPUNlP+VxVHSuyNiIAbEyqqzr2SZdrQmUXRpqf+vfTrWgvobJdy6n/HBNZGxEAqrnknjA7L86u6utRVLSZh6v0dcCVvkdEREREVMKwhYuIiIg0o8VaiKVxLUUmXERERKQZC3SwwNYxXKVvpvnSl0ISERERlTBs4SIiIiLNsEtRGRMuIiIi0ow2E5+WvoSr9D0iIiIiohKGLVxERESkGYukg8XWiU9tPL8kYsJFREREmrFo0KVYGic+ZcJFREREmrFIelhsHPRu6/klEROux9idePXLkcQ5uwiVnRBgVB0blaV+WRxAbLkeX4cEobILTaD1OtNN7IPA4qr+uaw+8IhQ2SLL9bxS6bBQ2YVxPtZXdWx6pvplgADgdry36lgXj3ShskWW6/nt5yeFykZDsXBR52/7C8VXcnugOnZnSqhQ2SLL9ZxLDxQq+xmh6IfcbqhfIkl3Tewzcrn+KdWx5gz19QDElutpUvEvobKp+DDhIiIiIs2YoYPZxolLbT2/JGLCRURERJphl6Ky0veIiIiIiEoYtnARERGRZsywvUvQrE1VShQmXERERKQZdikqK32PiIiIiKiEYQsXERERaYaLVytjwkVERESakaCDxcYxXFIpnBai9KWQRERERCUMW7iIiIhIM+xSVFb6HhEREREVG4uk02QT9fnnnyMoKAhGoxFhYWE4dOhQnrGLFi1Cy5Yt4e3tDW9vb0REROQbrwWdJEmSXa9AREREpV5CQgI8PT0xfH9XGNzF1kp9VHpSJuY0X4/4+HiYTAWvG/z999+jb9++WLBgAcLCwjBnzhysXr0a58+fh5+fX674yMhING/eHM2aNYPRaMSMGTOwZs0anDlzBuXLl7ep7nlhwkVEREQ2K86EKywsDI0bN8Znn30GALBYLKhYsSLefPNNjBkzpsDzzWYzvL298dlnn6Fv37421T0vHMNFREREmilsl+CjZQAPk7icDAYDDAaD1b6MjAwcPXoU7733nrxPr9cjIiICBw4cUHW9lJQUZGZmwsfHx6Z654djuIiIiEgzFug12QCgYsWK8PT0lLdp06blul5MTAzMZjP8/f2t9vv7+yMqKkpVnd99910EBgYiIiLC9icgD2zhIiIiohLpxo0bVl2Kj7ZuaWH69OlYtWoVdu/eDaPRqHn52ZhwERERkWbMkg5mG7sUs883mUwFjuEqW7YsHBwccPfuXav9d+/eRUBAQL7nzpw5E9OnT8f27dtRr149m+pcEHYpEhERkWaKeloIZ2dnNGzYEDt27PinDhYLduzYgaZNm+Z53scff4ypU6diy5YtaNSokU2PWQ22cBEREdFjbeTIkejXrx8aNWqEJk2aYM6cOUhOTsaAAQMAAH379kX58uXlMWAzZszAhAkTsHLlSgQFBcljvdzd3eHu7m6XOjLhIiIiIs1Ikh4WG2eKlwTPf/nllxEdHY0JEyYgKioKDRo0wJYtW+SB9NevX4de/0+ZX3zxBTIyMvDiiy9alTNx4kRMmjTJprrnhfNwERERkc2y5+Ea+NtLcLZxHq6MpEwsDv9B9TxcjwOO4SIiIiKyM3YpEhERkWYsEjSY+FSjypQgTLiIiIhIMxYNxnDZen5JxISLiIiINGOBDhbY2MJl4/klUelLIYmIiIhKGLZwERERkWa0nGm+NGHCRURERJrhGC5lpe8REREREZUwbOEiIiIizVggthZiXmWUNky4iIiISDOSBncpSqUw4WKXIhEREZGdsYWLiIiINGORNOhS5F2KRERERHnjXYrKSt8jIiIiIiph2MJFREREmmGXojImXERERKQZrqWojAkXERERaYYtXMqYcD3Gqv4wVXWsn1eSUNkpGU6qY71dU4XKNlvUDx100FuEyt719CyheCIioqLAhIsIQJe9b9qt7JvxnqpjA00JQmXfT3VVHRuf7CJU9rnuE4XiAaDyoo9Vx+rdsoTKlizqv/G6uKULla3XS6pjQ3zuC5W9vuVnQvEA0GLbO6pjkwW+HAGAm3Om6tj4VKNQ2Y4OZtWxTg5iX6YOdZgmFA8AjTaPVR374HRZobI9asaqjs0S+JIJAJJA646zo9j76HjnD4XiC4MtXMqYcBEREZFmmHAp47QQRERERHbGFi4iIiLSDFu4lDHhIiIiIs1IsH1aB/WjKh8f7FIkIiIisjO2cBEREZFm2KWojAkXERERaYYJlzJ2KRIRERHZGVu4iIiISDNs4VLGhOsxJrJcz91Yk1DZOp36e0QMTmIzHbs6qZ/p2kmvfubqoiIyuzsAJCSqjxddJklEZob93+4is8c7OIr9bS033FTHpgeK3eNkjlc/W/tVgfdGYYnMHq8X/L8U5ntNdeyR+5WEys4wO4hVxs5EZo93qCK2/NlT5f5SHXs5QWwW+wdp6leF0BfB61EUEy5lTLiIiIhIM5KkE1qeKK8yShuO4SIiIiKyM7ZwERERkWYs0Nk88amt55dETLiIiIhIMxzDpYxdikRERER2xhYuIiIi0gwHzStjwkVERESaYZeiMnYpEhER0WPv888/R1BQEIxGI8LCwnDo0KF841evXo0aNWrAaDSibt26+OWXX+xaPyZcREREpJnsLkVbNxHff/89Ro4ciYkTJ+LYsWOoX78+2rdvj3v37inG//7773jllVcwcOBAHD9+HN26dUO3bt1w+vRpLZ4CRUy4iIiISDPS312KtmyiCdfs2bMxePBgDBgwALVq1cKCBQvg6uqKJUuWKMbPnTsXHTp0wOjRo1GzZk1MnToVTz75JD777DMtngJFTLiIiIioREpISLDa0tPTc8VkZGTg6NGjiIiIkPfp9XpERETgwIEDiuUeOHDAKh4A2rdvn2e8Fjho/jGWIrDemsjaiIDYgMVybglCZadkOQvFF4Wb8Z6qY0XWRgSAkIBo1bF1vG4Llf1/0UGqY82Z9v9+JVnUv25E1kYEAH2FFNWxLi4ZQmVLbmmqY+253mU2N2f1642KrI0IAOsu1FUdm5mq/jMGALx8klXH6vUWobILw6NmrOpYkbURAWDPjRDVsakpYp95zgb1a5I6Ca5lWxQkAJKNSzxmn16xYkWr/RMnTsSkSZOs9sXExMBsNsPf399qv7+/P86dO6dYflRUlGJ8VFSUTfXODxMuIiIi0owFOug0mmn+xo0bMJlM8n6DwWBTucWJCRcRERFpRst5uEwmk1XCpaRs2bJwcHDA3bt3rfbfvXsXAQEBiucEBAQIxWuBY7iIiIjoseXs7IyGDRtix44d8j6LxYIdO3agadOmiuc0bdrUKh4Atm3blme8Fphw5ePMmTPo3bs3ypcvD4PBgMDAQERGRuLMmTPFXTUiIqISydY7FAszcerIkSOxaNEiLF++HGfPnsV//vMfJCcnY8CAAQCAvn374r333pPj33rrLWzZsgWzZs3CuXPnMGnSJBw5cgRDhw7V9LnIiV2Kefj555/xyiuvwMfHBwMHDkSVKlVw7do1LF68GD/++CNWrVqF559/vrirSUREVKJIkgaD5gXPf/nllxEdHY0JEyYgKioKDRo0wJYtW+SB8devX4de/08bU7NmzbBy5UqMHz8eY8eORbVq1bB27VrUqVPHtorngwmXgsuXL6NPnz4IDg7Gnj174OvrKx9766230LJlS/Tp0wcnT55EcHCwYhnJyclwcxO7C4uIiIgKZ+jQoXm2UO3evTvXvh49eqBHjx52rtU/2KWo4H//+x9SUlKwcOFCq2QLeDg478svv0RycjI+/vhjAMCkSZOg0+nw559/olevXvD29kaLFi3kc7799ls0adIErq6u8Pb2RqtWrbB161b5+Lp169C5c2cEBgbCYDAgJCQEU6dOhdlsLpoHTEREpJHimGn+ccAWLgUbNmxAUFAQWrZsqXi8VatWCAoKwqZNm6z29+jRA9WqVcNHH30E6e/20MmTJ2PSpElo1qwZpkyZAmdnZxw8eBA7d+5Eu3btAADLli2Du7s7Ro4cCXd3d+zcuRMTJkxAQkIC/ve//9n3wRIREWlIy7sUSxMmXI+Ij4/H7du38dxzz+UbV69ePaxfvx6JiYnyvvr162PlypXy75cuXcKUKVPw/PPP48cff7TqP5ZydFCvXLkSLi4u8u+vv/46Xn/9dcyfPx8ffPDBYz3vCBEREbFLMZfsBMrDwyPfuOzjCQn/zLL++uuvW8WsXbsWFosFEyZMsEq2AECn+yd7z5lsJSYmIiYmBi1btkRKSkqes+QSERGVRMVxl+LjgC1cj8hOpHK2XClRSsyqVKliFXP58mXo9XrUqlUr37LOnDmD8ePHY+fOnVYJHPCwxS0vIkuMGASXfxBZrud2kvplcQAg06w+zzcZc6+bZQ+BJvWPV3RpF5Hleg4KLNUDAHfvq3/uHe7av6XUxU393ys9UOw2JJHlehKi3YXKdnRTv5SOxUX9MkCFFZ9qVB175H4lobJFlutxMop9bjg6qF+ux91g//d2lkX9Z83lhLJCZYss1+PuLvaaSbiv/oYrY1n1r92iUhx3KT4OmHA9wtPTE+XKlcPJkyfzjTt58iTKly9vNQNuzpYqteLi4hAeHg6TyYQpU6YgJCQERqMRx44dw7vvvguLxf7rjREREZF9MeFS0KVLFyxatAj79u2zutsw2969e3Ht2jUMGTIk33JCQkJgsVjw559/okGDBooxu3fvxv379/Hzzz+jVatW8v6rV6/a9BiIiIiKw8MWLlsHzWtUmRKEY7gUjB49Gi4uLhgyZAju379vdSw2Nhavv/46XF1dMXr06HzL6datG/R6PaZMmZKrpSp70LyDg4PV7wCQkZGB+fPna/FQiIiIihSnhVDGFi4F1apVw/LlyxEZGYm6devmmmk+JiYG3333HUJCQvItp2rVqhg3bhymTp2Kli1bonv37jAYDDh8+DACAwMxbdo0NGvWDN7e3ujXrx+GDRsGnU6Hb775xioBIyIielxIf2+2llHaMOHKQ48ePVCjRg1MmzZNTrLKlCmDNm3aYOzYsaqn/58yZQqqVKmCTz/9FOPGjYOrqyvq1auHPn36AADKlCmDjRs3YtSoURg/fjy8vb3Ru3dvtG3bFu3bt7fnQyQiIqIiwoQrH3Xr1rWaVysvkyZNwqRJk/I8PmDAAHkBTSXNmjXDgQMHcu1nKxcRET1uOPGpMiZcREREpB32KSrioHkiIiIiO2MLFxEREWlHi7sM2aVIRERElDfONK+MXYpEREREdsYWrseYWWCdMFcnsfW2UrLUrxMmsjYiAPi7J6mOvRUvtk5jYd1PdbVb2f8nsD6iyNqIAGBOUL8uHrzF1sUrDL1e/ddSc7xA3QFIburXoxNZGxEAdNfUL8vVODT/Zb+04OhgVh2bYXYQKtvLJ1mgHmJLiyWnqf/cEHmMhSXSrfUgTWxpNmeD+veTyNqIABBW44rq2CtxZYTKLgq8S1EZEy4iIiLSjqSzfQxWKUy42KVIREREZGds4SIiIiLNcNC8MiZcREREpB1OfKqICRcRERFphoPmlXEMFxEREZGdsYWLiIiItFUKuwRtxYSLiIiINMMuRWXsUiQiIiKyM7ZwERERkXZ4l6IiJlyPMQe9+mU3nPT2W0bDZEwXihdZrsddsOzCik9Wv6xHZobY28acqb4h2eGuQahskeV6vPwSxcouhBCf+6pjr+rEPlG9XVNVx1pc1C8DBIgt17N+e5hQ2Z88IRQOAHASXFJHhF7gc8PdIPb+E1muJy1DbGmnwnB2VP/+0Au+Hp2c1JdtLCu21JTIcj2i9S4aur83W8soXdilSERERGRnbOEiIiIi7bBLURETLiIiItIOEy5F7FIkIiIisjO2cBEREZF2JN3DzdYyShkmXERERKQZSXq42VpGacOEi4iIiLTDMVyKOIaLiIiIyM6YcBEREZF2ssdw2brZSWxsLCIjI2EymeDl5YWBAwciKSkp3/g333wToaGhcHFxQaVKlTBs2DDEx8cLXZcJFxEREWlGJ2mz2UtkZCTOnDmDbdu2YePGjdizZw9ee+21PONv376N27dvY+bMmTh9+jSWLVuGLVu2YODAgULX5RguIiIi+lc4e/YstmzZgsOHD6NRo0YAgE8//RSdOnXCzJkzERgYmOucOnXq4KeffpJ/DwkJwYcffojevXsjKysLjo7qUikmXI+xXU/PKu4qlBrnuk8s7iqUCutbflbcVbC7wqyNKOpQh2n2v8i/wPHOHxZ3Ff6dNBw0n5CQYLXbYDDAYBBcczaHAwcOwMvLS062ACAiIgJ6vR4HDx7E888/r6qc+Ph4mEwm1ckWwC5FIiIi0pKGY7gqVqwIT09PeZs2zbYvI1FRUfDz87Pa5+joCB8fH0RFRakqIyYmBlOnTs23G1IJW7iIiIioRLpx4wZMJpP8e16tW2PGjMGMGTPyLevs2bM21ychIQGdO3dGrVq1MGnSJKFzmXARERGRdjTsUjSZTFYJV15GjRqF/v375xsTHByMgIAA3Lt3z2p/VlYWYmNjERAQkO/5iYmJ6NChAzw8PLBmzRo4OTkVWK+cmHARERGRdoph4lNfX1/4+voWGNe0aVPExcXh6NGjaNiwIQBg586dsFgsCAsLy/O8hIQEtG/fHgaDAevXr4fRaBSrIDiGi4iIiP4latasiQ4dOmDw4ME4dOgQ9u/fj6FDh6Jnz57yHYq3bt1CjRo1cOjQIQAPk6127dohOTkZixcvRkJCAqKiohAVFQWz2az62mzhIiIiIu2U8KV9VqxYgaFDh6Jt27bQ6/V44YUXMG/ePPl4ZmYmzp8/j5SUFADAsWPHcPDgQQBA1apVrcq6evUqgoKCVF2XCRcRERFpR4uZ4u0407yPjw9WrlyZ5/GgoCBIOVbPbt26tdXvhcWEi4iIiDSjxUzx9pxpvrhwDBcRERGRnbGFi4iIiLRTwsdwFRe2cBERERHZGRMuIiIiIjtjlyIRERFpRgcNBs1rUpOShQkXERERaaeETwtRXNilSERERGRnbOEiIiIi7fAuRUVMuIiIiEg7TLgUsUuRiIiIyM7YwkVERESa4dI+yphwERERkXbYpaiICRcRERFphwmXIo7hIiIiIrIztnARERGRZjiGSxkTLiIiItIOZ5pXxC5FIiIiIjtjCxcRERFph4PmFTHhIiIiIs1wDJcyJlyPsbBfx6iOvR/nLlS2syFLdayLc6ZQ2Q56i+pYi2A//pGOHwnFZ2u29V3VsY4C9QeAxHRn1bEehgyhskWenwyzg1DZhzpME4oHgMrLp6uOdTKqf40BYo/V25QiVLZO4NM90D1BqOz1LT8TigeAZ3aNUB2bJYmNDHFxVP9+jUtzESrbLPA3yhJ8PR7t9KFQPADU3/C+6tikc95CZbuExgnWRj0nB7Pdyj7eWfx5JG0w4SIiIiLtsEtRERMuIiIi0o4GXYqlMeHiXYpEREREdsYWLiIiItIOuxQVMeEiIiIi7TDhUsSEi4iIiDTDaSGUcQwXERERkZ0x4SIiIiKyM3YpEhERkXY4hksRW7iIiIiI7IwtXI8xkeV6PNxThcrOFFh2wyL4TUQnsPyHvohGToos1yO67EZahpPqWGdHsbIlgefSbLH/9yuR5XqyMsSWdtFHq18iKV5w+aWsdPUfhUXxmhRZric9S+xjvLbXHdWxJ9IrCJUt8tyIvHYLS2S5HqeQRKGyO1Y+qzp2/91gobJFlMSGIA6aV8aEi4iIiLRVChMmW7FLkYiIiMjO2MJFRERE2uGgeUVs4SIiIiLNZI/hsnWzl9jYWERGRsJkMsHLywsDBw5EUlKSqnMlSULHjh2h0+mwdu1aoesy4SIiIqJ/jcjISJw5cwbbtm3Dxo0bsWfPHrz22muqzp0zZw50usLd8MEuRSIiItJOCe5SPHv2LLZs2YLDhw+jUaNGAIBPP/0UnTp1wsyZMxEYGJjnuSdOnMCsWbNw5MgRlCtXTvjabOEiIiIizZTkLsUDBw7Ay8tLTrYAICIiAnq9HgcPHszzvJSUFPTq1Quff/45AgICCnVttnARERGRdjRs4UpISLDabTAYYDAYCl1sVFQU/Pz8rPY5OjrCx8cHUVFReZ43YsQINGvWDM8991yhr80WLiIiIiqRKlasCE9PT3mbNm2aYtyYMWOg0+ny3c6dO1eoOqxfvx47d+7EnDlzbHgkbOEiIiIiLWnYwnXjxg2YTCZ5d16tW6NGjUL//v3zLTI4OBgBAQG4d++e1f6srCzExsbm2VW4c+dOXL58GV5eXlb7X3jhBbRs2RK7d+/O97rZmHARERGRZrRc2sdkMlklXHnx9fWFr69vgXFNmzZFXFwcjh49ioYNGwJ4mFBZLBaEhYUpnjNmzBgMGjTIal/dunXxySef4Nlnny3wmtmYcD3GnA3q160TWRsRAFydM1XHOgiuW1cSJaarX6dPZG1EAPAzqZvfBQD8XcXWc7ud5Kk6NkOo5MKxiKyTKbA2IgBIAemqYz3d04TKNruqr7ebs/2fSRdH9e8/kbURAeBwdCXVsVG31a9FCAA+fgkFB/3NUgRre7qExqmOFVkbEQA2Xa2tOjYlwShUttFd/WvdxVAU7+zSo2bNmujQoQMGDx6MBQsWIDMzE0OHDkXPnj3lOxRv3bqFtm3b4uuvv0aTJk0QEBCg2PpVqVIlVKlSRfW1OYaLiIiItCNptNnJihUrUKNGDbRt2xadOnVCixYtsHDhQvl4ZmYmzp8/j5SUFE2vyxYuIiIi0k4JnocLAHx8fLBy5co8jwcFBUGS8q9AQceVlNoWrqCgoAIH0BEREREVhWJJuC5fvowhQ4YgODgYRqMRJpMJzZs3x9y5c5GamlocVSIiIiINlOSJT4tTkXcpbtq0CT169IDBYEDfvn1Rp04dZGRkYN++fRg9ejTOnDlj1ZdKREREj5ES3qVYXIo04bp69Sp69uyJypUrY+fOnVZrEb3xxhu4dOkSNm3aVJRVIiIiIrK7Iu1S/Pjjj5GUlITFixcrLvxYtWpVvPXWWwCApUuX4umnn4afnx8MBgNq1aqFL774Itc5kiThgw8+QIUKFeDq6oo2bdrgzJkzueImTZqkuML3smXLoNPpcO3aNav9mzdvRnh4ODw8PGAymdC4ceNcg+wOHjyIDh06wNPTE66urggPD8f+/futYvr374+goCBV9dm2bRtatGgBLy8vuLu7IzQ0FGPHjs11LhERUUnFLkVlRdrCtWHDBgQHB6NZs2YFxn7xxReoXbs2unbtCkdHR2zYsAH//e9/YbFY8MYbb8hxEyZMwAcffIBOnTqhU6dOOHbsGNq1a4eMjMLPTbJs2TK8+uqrqF27Nt577z14eXnh+PHj2LJlC3r16gXg4URpHTt2RMOGDTFx4kTo9Xo5Sdy7dy+aNGkidM0zZ86gS5cuqFevHqZMmQKDwYBLly7lSuCIiIhKNHYpKiqyhCshIQG3bt1SvfDjb7/9BhcXF/n3oUOHokOHDpg9e7accEVHR+Pjjz9G586dsWHDBrnFaNy4cfjoo48KVc/4+HgMGzYMTZo0we7du2E0/jNhXfZtoJIk4fXXX0ebNm2wefNm+bpDhgxB7dq1MX78eGzdulXoutu2bUNGRgY2b96MsmXLFqruRERExY4Jl6Ii61LMXvHbw8NDVXzOZCs+Ph4xMTEIDw/HlStXEB8fDwDYvn07MjIy8Oabb1p1zw0fPrzQ9dy2bRsSExMxZswYq2QLgHyNEydO4OLFi+jVqxfu37+PmJgYxMTEIDk5GW3btsWePXtgsYjNvp69RtO6deuEzyUiIqKSrchauLLXQkpMVLd0yf79+zFx4kQcOHAg12yv8fHx8PT0xF9//QUAqFatmtVxX19feHuLLUmR7fLlywCAOnXq5Blz8eJFAEC/fv3yjImPjxeqw8svv4yvvvoKgwYNwpgxY9C2bVt0794dL774IvR65bzYRWD5HYvgtwWR5XqcHcxCZadmii2NUxQ8BJbHcHYUe7wiy/XEpbsUHJSDyHOZlmn/t7u3Sf3MzPGCS0KJLNcTc6/gtdesCAwYcS+CpVTi0tS/Dk6kVxAqW2S5nrL+6pfqAQAXJ/WfSWaBZaCKwv67wULxIsv1lCkrtmTX/Vh31bHuLuqXASoqur83W8sobYo04QoMDMTp06cLjL18+TLatm2LGjVqYPbs2ahYsSKcnZ3xyy+/4JNPPilUC5DSgHkAMJvF/nkCkK//v//9Dw0aNFCMcXd3F7qui4sL9uzZg127dmHTpk3YsmULvv/+ezz99NPYunUrHBzE1kIkIiIqFuxSVFSkg+a7dOmChQsX4sCBA2jatGmecRs2bEB6ejrWr1+PSpX+WWh1165dVnGVK1cG8LDFKTj4n28n0dHRePDggVVsdmtTXFyc3H0HQG4lyxYSEgIAOH36NKpWrapYv+wYk8mEiIiIPB9H9nXj4uJy7X/0ugCg1+vRtm1btG3bFrNnz8ZHH32EcePGYdeuXQVeh4iIiEquIp0W4p133oGbmxsGDRqEu3fv5jp++fJlzJ07V27NyblWUXx8PJYuXWoVHxERAScnJ3z66adWsXPmzMlVdnaStGfPHnlfcnIyli9fbhXXrl07eHh4YNq0aUhLs+7CyL5Gw4YNERISgpkzZyIpKSnXtaKjo62uGx8fj5MnT8r77ty5gzVr1lidExsbm6uc7Naz9PSS12RMRESkhNNCKCvSFq6QkBCsXLkSL7/8MmrWrGk10/zvv/+O1atXo3///hg5ciScnZ3x7LPPYsiQIUhKSsKiRYvg5+eHO3fuyOX5+vri7bffxrRp09ClSxd06tQJx48fV7zTr127dqhUqRIGDhyI0aNHw8HBAUuWLIGvry+uX78ux5lMJnzyyScYNGgQGjdujF69esHb2xt//PEHUlJSsHz5cuj1enz11Vfo2LEjateujQEDBqB8+fK4desWdu3aBZPJhA0bNgAAevbsiXfffRfPP/88hg0bhpSUFHzxxReoXr06jh07Jl93ypQp2LNnDzp37ozKlSvj3r17mD9/PipUqIAWLVrY+S9DRESkEXYpKirypX26du2KkydP4n//+x/WrVuHL774AgaDAfXq1cOsWbMwePBgGAwG/Pjjjxg/fjzefvttBAQE4D//+Q98fX3x6quvWpX3wQcfwGg0YsGCBdi1axfCwsKwdetWdO7c2SrOyckJa9aswX//+1+8//77CAgIwPDhw+Ht7Y0BAwZYxQ4cOBB+fn6YPn06pk6dCicnJ9SoUQMjRoyQY1q3bo0DBw5g6tSp+Oyzz5CUlISAgACEhYVhyJAhclyZMmWwZs0ajBw5Eu+88w6qVKmCadOm4eLFi1YJV9euXXHt2jUsWbIEMTExKFu2LMLDwzF58mR4enpq+ScgIiKiIqaTcvbF0WOl4S/jVMeK3qVodMpSHVuS7lI81GFaoc5rtX206th0s9gNDJVNDwoO+pvoXYr3U9xUx4repXjmuclC8YDYazI+Sf1dXkDJuUuxYmDu7v/87HvmY7G6AGi29V3VsXrBvpcbN8uoji1JdykeaDdDKB4A6qyfoDrWZBQbunE7Sv3dnva8S7FsGbGyj3Qs3ByVaiQkJMDT0xO1h3wEB2ex9/ejzBlpOPPlWMTHx8uzHDzuiryFi4iIiEovLcZglcYxXEU6aJ6IiIjo34gtXERERKQdDppXxISLiIiINMMuRWVMuIiIiEg7bOFSxITrMSay3qHOjuuWid51WNY1WXWsyJpytrAIPD+S4HN5O0n9tB6iz6XIa0C03oWhE/hampUu9vFjdhWov+DXY12S+rpU84wuOMhGInfwid6l6OOn/s5DkbsOAbHXr8hrt7CcBO+gFmF0V39Xo8hdhwBQudx91bGORfA8kjaYcBEREZFm2KWojAkXERERaYddioo4LQQRERGRnbGFi4iIiLTDFi5FTLiIiIhIMxzDpYxdikRERER2xhYuIiIi0g67FBUx4SIiIiLN6CQJOsm2jMnW80sidikSERER2RlbuIiIiEg77FJUxITrMSayHI3o8h/2JLJcj49Lih1r8o8Ms4PqWLNFrGE4QyA2LVPsLSmyXE9RLKUS6K5+2RjR16Sbs/pn0t0g8qyLLdfz++Z6QmWjiVg4AGQJvB5Fl2yyCLx+RZYYAsReY5nmktXBIvoJ6SLwGnN3Ub8MECC2XI+T3n7LFxUW71JUxoSLiIiItMMWLkUl6ysGERERUSnEFi4iIiLSDLsUlTHhIiIiIu2wS1ERuxSJiIjoXyM2NhaRkZEwmUzw8vLCwIEDkZSUVOB5Bw4cwNNPPw03NzeYTCa0atUKqampqq/LhIuIiIg0k92laOtmL5GRkThz5gy2bduGjRs3Ys+ePXjttdfyPefAgQPo0KED2rVrh0OHDuHw4cMYOnQo9Hr1aRS7FImIiEg7JbhL8ezZs9iyZQsOHz6MRo0aAQA+/fRTdOrUCTNnzkRgYKDieSNGjMCwYcMwZswYeV9oaKjQtdnCRURERCVSQkKC1ZaeLjan2aMOHDgALy8vOdkCgIiICOj1ehw8eFDxnHv37uHgwYPw8/NDs2bN4O/vj/DwcOzbt0/o2ky4iIiISFNadSdWrFgRnp6e8jZt2jSb6hUVFQU/Pz+rfY6OjvDx8UFUVJTiOVeuXAEATJo0CYMHD8aWLVvw5JNPom3btrh48aLqa7NLkYiIiLQjSQ83W8sAcOPGDZhMJnm3wWBQDB8zZgxmzJiRb5Fnz54tVFUslocz/w8ZMgQDBgwAADzxxBPYsWMHlixZojoJZMJFREREmtFyHi6TyWSVcOVl1KhR6N+/f74xwcHBCAgIwL1796z2Z2VlITY2FgEBAYrnlStXDgBQq1Ytq/01a9bE9evXC6xbNiZcj7EjHT8q7iqUGoc62NZMTQ+tb/lZcVfB/gqxNqKoo50+tP9F/gWOd+bz+G/h6+sLX1/fAuOaNm2KuLg4HD16FA0bNgQA7Ny5ExaLBWFhYYrnBAUFITAwEOfPn7faf+HCBXTs2FF1HTmGi4iIiLQjabTZQc2aNdGhQwcMHjwYhw4dwv79+zF06FD07NlTvkPx1q1bqFGjBg4dOgQA0Ol0GD16NObNm4cff/wRly5dwvvvv49z585h4MCBqq/NFi4iIiLSjM7ycLO1DHtZsWIFhg4dirZt20Kv1+OFF17AvHnz5OOZmZk4f/48UlJS5H3Dhw9HWloaRowYgdjYWNSvXx/btm1DSEiI6uvqJMnWkW1ERET0b5eQkABPT080fv4DODoZbSorKzMNh9eMR3x8vKoxXI8DtnARERGRdkrwxKfFiQkXERERaUbLuxRLEw6aJyIiIrIztnARERGRdjSc+LQ0YcJFREREmmGXojJ2KRIRERHZGVu4iIiISDu8S1EREy4iIiLSDLsUlTHhIiIiIu1w0LwijuEiIiIisjO2cBEREZFm2KWojAkXERERaYeD5hWxS5GIiIjIztjCRURERJphl6IyJlxERESkHYv0cLO1jFKGXYpEREREdsYWLiIiItIOB80rYsJFREREmtFBgzFcmtSkZGGXIhEREZGdsYWLiIiItMOlfRQx4SIiIiLNcFoIZUy4iIiISDscNK+IY7iIiIiI7IwtXI+x8B1vq46NSXITKlsn0J7rZsgQKluEXrBd+UC7GXaqCRERqaGTJOhsHINl6/klERMuIgAdfntLdWxsqqtQ2akZTqpjPV3ShMoWkWkRa9A+2H668DWCZ89SHetQMUWo7MxU9c+j0T1dqGyDU5bqWD/3JKGyt7X5RCgeALrv/4/q2HspHkJluzhmqo6NSRF7rRsFnkeLJHbjf2Fej622j1Yde/2Cv1DZzn7qX78uBvXPOQBkZDmojjU6i5V9vPOHQvGFYvl7s7WMUoZdikRERER2xhYuIiIi0gy7FJUx4SIiIiLt8C5FRexSJCIiIrIztnARERGRdjjTvCImXERERKQZzjSvjF2KRERERHbGFi4iIiLSDrsUFTHhIiIiIs3oLA83W8sobZhwPcZElutx1Iu9epOSjapjDY5mobIdBOoiOjt6YYnMHi86Q3ZqqrPqWNGZ5kWeH3MRPJcis8dnJKl/XgDA+bb6mebTKol9O07NVP96F/37F4bI7PEBbglCZXs7q/8bpZvF/kVkFdH7VS2R2eN1nmJLlAWVjVUdK7o6hQhnwc/fIsEWLkUl691BREREVAox4SIiIiLtSBptdhIbG4vIyEiYTCZ4eXlh4MCBSErKf43UqKgo9OnTBwEBAXBzc8OTTz6Jn376Sei6TLiIiIhIM9lL+9i62UtkZCTOnDmDbdu2YePGjdizZw9ee+21fM/p27cvzp8/j/Xr1+PUqVPo3r07XnrpJRw/flz1dZlwERER0b/C2bNnsWXLFnz11VcICwtDixYt8Omnn2LVqlW4fft2nuf9/vvvePPNN9GkSRMEBwdj/Pjx8PLywtGjR1VfmwkXERERaSd70LytG4CEhASrLT093aaqHThwAF5eXmjUqJG8LyIiAnq9HgcPHszzvGbNmuH7779HbGwsLBYLVq1ahbS0NLRu3Vr1tZlwERERkXYkABYbt797FCtWrAhPT095mzZtmk1Vi4qKgp+fn9U+R0dH+Pj4ICoqKs/zfvjhB2RmZqJMmTIwGAwYMmQI1qxZg6pVq6q+NqeFICIiohLpxo0bMJlM8u8Gg0ExbsyYMZgxY0a+ZZ09e7bQ9Xj//fcRFxeH7du3o2zZsli7di1eeukl7N27F3Xr1lVVBhMuIiIi0owWg96zzzeZTFYJV15GjRqF/v375xsTHByMgIAA3Lt3z2p/VlYWYmNjERAQoHje5cuX8dlnn+H06dOoXbs2AKB+/frYu3cvPv/8cyxYsEDFI2LCRURERFqSoMHEp2Lhvr6+8PX1LTCuadOmiIuLw9GjR9GwYUMAwM6dO2GxWBAWFqZ4TkrKw8mC9XrrUVgODg6wWNRP5M0xXERERPSvULNmTXTo0AGDBw/GoUOHsH//fgwdOhQ9e/ZEYGAgAODWrVuoUaMGDh06BACoUaMGqlatiiFDhuDQoUO4fPkyZs2ahW3btqFbt26qr82Ei4iIiLSj4V2K9rBixQrUqFEDbdu2RadOndCiRQssXLhQPp6ZmYnz58/LLVtOTk745Zdf4Ovri2effRb16tXD119/jeXLl6NTp06qr8suxceYTqf+BSmyNiIAuLiqv/XWx0X92mwAkJSpfg09i9lBqOzCSs1Qv06fyNqIAFC93L2Cg/7mY0gWKvtOiqfq2MR05cGmWspMVf88iqyNCAAZAVmqY908xG4d1wus7+ninClUdmG4OKq/hsjaiACw93qI6tjUBLHPjYByD1THpmfa/9+Ps5/650ZkbUQAiE5Wv5btg/vq18YEAIOb+nUdRV67RcYCwNYlR+34sHx8fLBy5co8jwcFBUF6JOGrVq2a8Mzyj2LCRURERJrRctB8acIuRSIiIiI7YwsXERERaUeLMVhs4QKWLVsGnU6HI0eOKB5v3bo16tSpY3PFitvu3buh0+nkzWAwwN/fH61bt8ZHH32E6OjoXOdkPzdGoxG3bt3KdTy/58ZsNiMwMBA6nQ6bN2/W/PEQEREViRI+aL64sEuxAMOGDcM333yDhQsXYvTo0fDx8cHEiRNRs2ZN7Ny5U/Gc9PR0TJ8+Xeg6O3fuxJ07dxAUFIQVK1ZoUXUiIiIqIdilWICWLVvixRdftNr3xx9/oF27dnjhhRfw559/oly5clbHGzRogEWLFuG9996T5/UoyLfffosnn3wS/fr1w9ixY5GcnAw3N/V3wRAREZUI7FJUZNcWrvDwcNSvX1/xWGhoKNq3by//brFYMGfOHNSuXRtGoxH+/v4YMmQIHjywvs04KCgIXbp0wb59+9CkSRMYjUYEBwfj66+/toqLjY3F22+/jbp168Ld3R0mkwkdO3bEH3/8YfPjql+/PubMmYO4uDh89tlnuY6PHTsWZrNZdStXamoq1qxZg549e+Kll15Camoq1q1bZ3M9iYiIipytC1dnb6VMoROu+Ph4xMTE5NoyM/+ZP6ZPnz44efIkTp8+bXXu4cOHceHCBfTu3VveN2TIEIwePRrNmzfH3LlzMWDAAKxYsQLt27e3KhMALl26hBdffBHPPPMMZs2aBW9vb/Tv3x9nzpyRY65cuYK1a9eiS5cumD17NkaPHo1Tp04hPDwct2/fLuzDlr344otwcXHB1q1bcx2rUqUK+vbti0WLFqm61vr165GUlISePXsiICAArVu3ZrciERFRKVLohCsiIkJeuyjn9vvvv8sxPXr0gNFoxLfffmt17rfffgs3Nzd0794dALBv3z589dVXWL58ORYuXIghQ4Zg+vTp+Omnn3D48GGsXr3a6vzz589j9erV+PDDD/HGG29gy5YtcHZ2xtKlS+WYunXr4sKFC5g2bRpee+01vP/++9i3bx/S0tKwePHiwj5smZOTE6pXr47Lly8rHh83bhyysrIKXL0cePh8NGvWDBUrVgQA9OzZE1u3blUcmE9ERFSSZc/DZetW2hQ64fr888+xbdu2XFu9evXkGE9PTzz33HP47rvv5FlbzWYzvv/+e3Tr1k0eo7R69Wp4enrimWeesWota9iwIdzd3bFr1y6ra9eqVQstW7aUf/f19UVoaCiuXLki7zMYDPJCk2azGffv34e7uztCQ0Nx7Nixwj5sK+7u7khMTFQ8FhwcjD59+mDhwoW4c+dOnmXcv38fv/76K1555RV53wsvvACdTocffvhBk3oSEREVGd6lqKjQg+abNGmCRo0a5drv7e2NmJgY+fe+ffvi+++/x969e9GqVSts374dd+/eRZ8+feSYixcvIj4+Hn5+forXunfPemmUSpUqKV4353gvi8WCuXPnYv78+bh69SrMZrN8rEyZMvLPUVFRVuV4enrCxcUlr4dtJSkpCR4eeS/ZMH78eHzzzTeYPn065s6dqxjz/fffIzMzE0888QQuXbok7w8LC8OKFSvwxhtv5Fm+m0H98g8GR3PBQTmILNeTbhZ7GWVkqY/PMhfNjbSeLml2iQXElus5F+svVLabs/rXQIrA8kWFZXRXv6ROWiWxD1SR5XqS79nvhpPMMql2KztbTIqr6ljR95/Icj0Vy98XKttBYJkZg4PYZ1JhuBjUL5EUm6r+OQfElutx8xR7zSTdV18XD1exzyMqPna/S7F9+/bw9/fHt99+i1atWuHbb79FQEAAIiIi5BiLxQI/P788xy35+vpa/e7goLy+Xs61jz766CO8//77ePXVVzF16lT4+PhAr9dj+PDhsFj++VB49A7DpUuXon///gU+rszMTFy4cCHfOceCg4PRu3dvLFy4EGPGjFGMyX7MzZs3Vzx+5coVBAcHF1gfIiKiEsEiAQJr/eZZRilj94TLwcEBvXr1wrJlyzBjxgysXbsWgwcPtkqaQkJCsH37djRv3lx161JBfvzxR7Rp0ybXeK24uDiULVtW/n3btm1Wx2vXrq26/NTUVKs7LZWMHz8e3377reJYrqtXr+L333/H0KFDER4ebnXMYrGgT58+WLlyJcaPH6+qTkRERMWO00IoKpJ5uPr06YNPPvkEQ4YMQVJSktXdiQDw0ksvYf78+Zg6dSo++ugjq2NZWVlISkqCl5eX0DUdHBxyrfa9evVq3Lp1C1WrVpX35WxpU+uPP/7A8OHD4e3tnW+XH/Awmezduze+/PJLVK5cGY6O/zzl2a1b77zzjjxgPqevvvoKK1asYMJFRESPES3GYDHhKpQnnngCderUwerVq1GzZk08+eSTVsfDw8MxZMgQTJs2DSdOnEC7du3g5OSEixcvYvXq1Zg7d26uyUcL0qVLF0yZMgUDBgxAs2bNcOrUKaxYsUK4e27v3r1IS0uTB97v378f69evh6enJ9asWYOAgIACyxg3bhy++eYbnD9/3qoFbcWKFWjQoIFisgUAXbt2xZtvvoljx47les6IiIjo8VFkS/v07dsXAKwGy+e0YMECLFy4EPfu3cPYsWPx3nvvYefOnejdu3ee45vyM3bsWIwaNQq//vor3nrrLRw7dgybNm3KM7nJy7x589CnTx8MHjwY06dPR0xMDCZPnoyzZ8/m6gbMS9WqVXO16h07dgznzp3Ds88+m+d52ccenVaDiIioxOJdioqEW7j69++f76Dy3bt3K+53dnaGTqdDZGRknucOHjwYgwcPzvf6165dU3Vdg8GAmTNnYubMmarq96jWrVvn6pIsSH7PzbJly7Bs2TL59yeffLLA8itXrixcByIiomJlkWBzl2ApHDRfJC1ckiRh8eLFCA8PV5zSgYiIiKg0s+sYruTkZKxfvx67du3CqVOnuD4gERFRaSdZHm62llHK2DXhio6ORq9eveDl5YWxY8eia9eu9rwcERERFTdOC6HIrglXUFAQxyARERHRv16RTAtBRERE/xIcNK+ICde/hMgaZwCQlOmsOlZkbURAbJ3GO4nq1ysrKpkWsXtN7qR4qo4VWRsRAJIz1P+dEu+6C5VdGAanLNWxqZnq1/QDAL3ga1iEc4zycmFKqoSKrS9YGEaB5zFL8PUYUO5BwUF/E/3ciEvRZqUQrWRkqf+7ijK4qX+viqyNCADBQfcKDvqbo64EjnVil6KiIpuHi4iIiOjfii1cREREpB0JGrRwaVKTEoUJFxEREWmHXYqKmHARERGRdiwWADaOLbOUwLFpNuIYLiIiIiI7YwsXERERaYddioqYcBEREZF2mHApYpciERERkZ2xhYuIiIi0w5nmFTHhIiIiIs1IkgWSZNtdhraeXxIx4XqM6XXqvwGILkdjMatfEiPLLLjUjcByPRU844XKLiyR58cs+FwmphtUx6ZkOImVLbBcT9ny9n8u/dyTVMdaJJ1Q2S7OmapjM8ukCpUtslzP5f2VhcpGuFg4IP7ciEjPVP+xb3Aw260ejg72/4dqFHjNODuKPVaRpaY8XNOEyhZZrifdzH/jjwv+pYiIiEg7kmR7l2ApHDTPhIuIiIi0I2kwhqsUJly8S5GIiIj+NT788EM0a9YMrq6u8PLyUnWOJEmYMGECypUrBxcXF0RERODixYtC12XCRURERNqxWLTZ7CQjIwM9evTAf/7zH9XnfPzxx5g3bx4WLFiAgwcPws3NDe3bt0damvrxeexSJCIiIu2U8C7FyZMnAwCWLVumsioS5syZg/Hjx+O5554DAHz99dfw9/fH2rVr0bNnT1XlsIWLiIiINCNZLJpsAJCQkGC1paenF/njuXr1KqKiohARESHv8/T0RFhYGA4cOKC6HCZcREREVCJVrFgRnp6e8jZt2rQir0NUVBQAwN/f32q/v7+/fEwNdikSERGRdjTsUrxx4wZMJpO822BQntdwzJgxmDFjRr5Fnj17FjVq1LCtXjZgwkVERETasUiAwMTciv5OuEwmk1XClZdRo0ahf//++cYEBwcXqioBAQEAgLt376JcuXLy/rt376JBgwaqy2HCRURERI81X19f+Pr62qXsKlWqICAgADt27JATrISEBBw8eFDoTkeO4SIiIiLtSBIgWWzc7HeX4vXr13HixAlcv34dZrMZJ06cwIkTJ5CU9M/SZDVq1MCaNWsAADqdDsOHD8cHH3yA9evX49SpU+jbty8CAwPRrVs31dfVSVIpnM6ViIiIilRCQgI8PT3RxvFFOOrE1oV9VJaUiV1ZPyI+Pl5Vl6KI/v37Y/ny5bn279q1C61btwbwMMlaunSp3E0pSRImTpyIhQsXIi4uDi1atMD8+fNRvXp11ddlwkVEREQ2e1wSruLCMVxERESkHckCwMaZ4iX7zTRfXJhwERERkWYkiwTJxrsUS2PnGwfNExEREdkZW7iIiIhIM1lSus1dglnI1Kg2JQcTLiIiIrKZs7MzAgICsC/qF03KCwgIgLOzsyZllQS8S5GIiIg0kZaWhoyMDE3KcnZ2htFo1KSskoAJFxEREZGdcdA8ERERkZ0x4SIiIiKyMyZcRERERHbGhIuIiIjIzphwEREREdkZEy4iIiIiO2PCRURERGRn/w9Ww9N5wZMofQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 25 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot nxn image\n",
    "rows = ['GT','Orca','Caduceus','Hyena-DNA']\n",
    "out_collection=[]\n",
    "normmat_bydist = np.exp(\n",
    "                np.load(\"orca/resources/resources/4DNFI9GMP2J8.rebinned.mcool.expected.res1000.npy\")\n",
    "            )[:6]\n",
    "normmat = normmat_bydist[np.abs(np.arange(6)[:, None] - np.arange(6)[None, :])]\n",
    "normmat_r = torch.from_numpy(normmat).float().cuda()\n",
    "eps = torch.min(normmat_r)\n",
    "target_r=torch.log(((target + eps) / (normmat_r + eps)))\n",
    "out_collection.append(target_r.cpu().detach().numpy())\n",
    "out_collection.append(prediction.cpu().detach().numpy())\n",
    "out_collection.append(mat_mamba.cpu().detach().numpy())\n",
    "out_collection.append(mat_hyena.cpu().detach().numpy())\n",
    "_min, _max = np.amin(out_collection), np.amax(out_collection)\n",
    "fig,axes = plt.subplots(4,5)\n",
    "for ax, row in zip(axes[:,0], rows):\n",
    "    ax2 = ax.twinx()\n",
    "    ax2.yaxis.set_label_position('left')\n",
    "    ax2.spines['left'].set_visible(False)\n",
    "    ax2.spines['right'].set_visible(False)\n",
    "    ax2.spines['top'].set_visible(False)\n",
    "    ax2.spines['bottom'].set_visible(False)\n",
    "    #set ticks invisible\n",
    "    ax2.tick_params(axis='both', which='both', length=0)\n",
    "    ax2.set_yticks([])\n",
    "    ax2.set_ylabel(row, rotation=0, size='large',\n",
    "                   ha='right', va='center')\n",
    "for i in range(4):\n",
    "    for j in range(5):\n",
    "        ax_n=axes[i,j].imshow(out_collection[i][j], vmin = _min, vmax = _max)\n",
    "        axes[i,j].axis('off')\n",
    "        axes[i,j].autoscale(False)\n",
    "fig.colorbar(ax_n,ax=axes,orientation='vertical')\n",
    "\n",
    "plt.suptitle('Genomic Structure Prediction for HFF',fontsize=13,y=0.95,fontweight='bold')\n",
    "plt.show()\n",
    "        #show value bar\n",
    "fig.savefig('assets/genomic_structure_hff.pdf',dpi=300,bbox_inches='tight')\n",
    "fig.savefig('assets/genomic_structure_hff.png',dpi=300,bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}