StochasticFactorizatoin.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6bad5e5e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count= 2257\n",
      "98.82 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAANmUlEQVR4nO3df6zdd13H8efLlhEVkB+9kNkftpiK9g/AcR0YFTFEaOcflYQ/NgzDhaVZQg3+YbIaEiXhH5BgDGHQVGwAY+g/TClQmIaoxMBknRnbytJx2ZBdurBODBJJnIW3f5zv7OFw7z3ntqfr7rvPR3Jyz/f7/dxzP+eTs+fOvufHUlVIkja+n7jcE5AkzYdBl6QmDLokNWHQJakJgy5JTWy+XH94y5YttXPnzsv15yVpQ7r77rsfr6qFlY5dtqDv3LmTkydPXq4/L0kbUpJ/X+2Yp1wkqQmDLklNGHRJasKgS1ITBl2SmjDoktTE1KAnOZrksST3r3I8Sd6fZCnJvUmumf80JUnTzPIM/SPA3jWO7wN2D5cDwIcuflqSpPWaGvSq+gLwnTWG7Ac+ViN3As9NcvW8JihJms08Pim6FXhkbHt52Pfo5MAkBxg9i2fHjh0X/Ad3HvrMBf/u09k33v077Dz0mafVzyfN4/a0tgtZp6fDY+Sp+Pl0dDFzG/9na57m8aJoVti34v8GqaqOVNViVS0uLKz4VQSSpAs0j6AvA9vHtrcBZ+Zwu5KkdZhH0I8DNw7vdnkV8N2q+rHTLZKkS2vqOfQkHwdeA2xJsgz8KfAMgKo6DJwArgOWgO8DN12qyUqSVjc16FV1w5TjBbxtbjOSJF0QPykqSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJamJmYKeZG+S00mWkhxa4fjPJPlUkq8kOZXkpvlPVZK0lqlBT7IJuA3YB+wBbkiyZ2LY24CvVtXLgNcA70ty1ZznKklawyzP0K8Flqrqoap6AjgG7J8YU8CzkwR4FvAd4NxcZypJWtMsQd8KPDK2vTzsG/cB4JeAM8B9wNur6oeTN5TkQJKTSU6ePXv2AqcsSVrJLEHPCvtqYvv1wD3AzwIvBz6Q5Dk/9ktVR6pqsaoWFxYW1jlVSdJaZgn6MrB9bHsbo2fi424Cbq+RJeBh4BfnM0VJ0ixmCfpdwO4ku4YXOq8Hjk+M+SbwWoAkLwJeAjw0z4lKkta2edqAqjqX5CBwB7AJOFpVp5LcMhw/DLwL+EiS+xidorm1qh6/hPOWJE2YGnSAqjoBnJjYd3js+hngdfOdmiRpPfykqCQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJmYKepK9SU4nWUpyaJUxr0lyT5JTSf55vtOUJE2zedqAJJuA24DfBpaBu5Icr6qvjo15LvBBYG9VfTPJCy/RfCVJq5jlGfq1wFJVPVRVTwDHgP0TY94E3F5V3wSoqsfmO01J0jSzBH0r8MjY9vKwb9wvAM9L8k9J7k5y47wmKEmazdRTLkBW2Fcr3M4rgNcCPwl8KcmdVfXgj9xQcgA4ALBjx471z1aStKpZnqEvA9vHtrcBZ1YY87mq+u+qehz4AvCyyRuqqiNVtVhViwsLCxc6Z0nSCmYJ+l3A7iS7klwFXA8cnxjzSeA3kmxO8lPAK4EH5jtVSdJapp5yqapzSQ4CdwCbgKNVdSrJLcPxw1X1QJLPAfcCPwQ+XFX3X8qJS5J+1Czn0KmqE8CJiX2HJ7bfC7x3flOTJK2HnxSVpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJmYKepK9SU4nWUpyaI1xv5LkB0neOL8pSpJmMTXoSTYBtwH7gD3ADUn2rDLuPcAd856kJGm6WZ6hXwssVdVDVfUEcAzYv8K4PwA+ATw2x/lJkmY0S9C3Ao+MbS8P+/5fkq3AG4DDa91QkgNJTiY5efbs2fXOVZK0hlmCnhX21cT2XwC3VtUP1rqhqjpSVYtVtbiwsDDjFCVJs9g8w5hlYPvY9jbgzMSYReBYEoAtwHVJzlXV381jkpKk6WYJ+l3A7iS7gG8B1wNvGh9QVbuevJ7kI8CnjbkkPbWmBr2qziU5yOjdK5uAo1V1Ksktw/E1z5tLkp4aszxDp6pOACcm9q0Y8qr6/YufliRpvfykqCQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJmYKepK9SU4nWUpyaIXjv5fk3uHyxSQvm/9UJUlrmRr0JJuA24B9wB7ghiR7JoY9DPxmVb0UeBdwZN4TlSStbZZn6NcCS1X1UFU9ARwD9o8PqKovVtV/Dpt3AtvmO01J0jSzBH0r8MjY9vKwbzVvBT670oEkB5KcTHLy7Nmzs89SkjTVLEHPCvtqxYHJbzEK+q0rHa+qI1W1WFWLCwsLs89SkjTV5hnGLAPbx7a3AWcmByV5KfBhYF9V/cd8pidJmtUsz9DvAnYn2ZXkKuB64Pj4gCQ7gNuBN1fVg/OfpiRpmqnP0KvqXJKDwB3AJuBoVZ1Kcstw/DDwJ8ALgA8mAThXVYuXbtqSpEmznHKhqk4AJyb2HR67fjNw83ynJklaDz8pKklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU3MFPQke5OcTrKU5NAKx5Pk/cPxe5NcM/+pSpLWMjXoSTYBtwH7gD3ADUn2TAzbB+weLgeAD815npKkKWZ5hn4tsFRVD1XVE8AxYP/EmP3Ax2rkTuC5Sa6e81wlSWtIVa09IHkjsLeqbh623wy8sqoOjo35NPDuqvqXYfvzwK1VdXLitg4wegYP8BLg9AXOewvw+AX+bjeuxXmuxYjrcF7Htfi5qlpY6cDmGX45K+yb/LfALGOoqiPAkRn+5toTSk5W1eLF3k4HrsV5rsWI63DelbYWs5xyWQa2j21vA85cwBhJ0iU0S9DvAnYn2ZXkKuB64PjEmOPAjcO7XV4FfLeqHp3zXCVJa5h6yqWqziU5CNwBbAKOVtWpJLcMxw8DJ4DrgCXg+8BNl27KwBxO2zTiWpznWoy4DuddUWsx9UVRSdLG4CdFJakJgy5JTWy4oE/7GoJuknwjyX1J7klyctj3/CT/kORrw8/njY3/42FtTid5/eWb+cVLcjTJY0nuH9u37vue5BXDGi4NX1Gx0ttsn9ZWWYt3JvnW8Ni4J8l1Y8darkWS7Un+MckDSU4lefuw/4p8XPyYqtowF0Yvyn4deDFwFfAVYM/lntclvs/fALZM7Psz4NBw/RDwnuH6nmFNngnsGtZq0+W+Dxdx318NXAPcfzH3Hfgy8KuMPi/xWWDf5b5vc1qLdwJ/tMLYtmsBXA1cM1x/NvDgcH+vyMfF5GWjPUOf5WsIrgT7gY8O1z8K/O7Y/mNV9T9V9TCjdx1d+9RPbz6q6gvAdyZ2r+u+D19B8Zyq+lKN/in+2NjvbBirrMVq2q5FVT1aVf82XP8e8ACwlSv0cTFpowV9K/DI2PbysK+zAv4+yd3DVycAvKiG9/kPP1847L8S1me9933rcH1yfxcHh284PTp2muGKWIskO4FfBv4VHxfAxgv6TF8x0MyvVdU1jL7R8m1JXr3G2CtxfZ602n3vvCYfAn4eeDnwKPC+YX/7tUjyLOATwB9W1X+tNXSFfa3WYtxGC/oV9xUDVXVm+PkY8LeMTqF8+8lvsxx+PjYMvxLWZ733fXm4Prl/w6uqb1fVD6rqh8Bfcv70Wuu1SPIMRjH/m6q6fdjt44KNF/RZvoagjSQ/neTZT14HXgfcz+g+v2UY9hbgk8P148D1SZ6ZZBej76f/8lM760tuXfd9+M/v7yV51fAuhhvHfmdDm/iK6jcwemxA47UY5v1XwANV9edjh3xcwMZ6l8votQuuY/TK9teBd1zu+Vzi+/piRq/QfwU49eT9BV4AfB742vDz+WO/845hbU6zwV+1Bz7O6FTC/zJ6RvXWC7nvwCKj2H0d+ADDJ6Q30mWVtfhr4D7gXkbhurr7WgC/zujUyL3APcPluiv1cTF58aP/ktTERjvlIklahUGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1IT/wdZCXDVi2th1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD5CAYAAADGMZVsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtv0lEQVR4nO3deXhU1fnA8e9LICCIbIZ9CSIFEWUxoigqSFFELLZVC60WVyzV7laxFrV1KXWr6w/EorivuFVQVlEQBAKy76uELQGEBAKEJO/vj5mEyeTOPsls7+d5eDJz7rn3vjfAO2fOPfccUVWMMcYktxqxDsAYY0zVs2RvjDEpwJK9McakAEv2xhiTAizZG2NMCrBkb4wxKaBmoAoi0gZ4DWgOlALjVfUZEWkMvAtkAluB61T1B4f9BwLPAGnAf1V1TKBznnrqqZqZmRn8VRhjTIpbvHjxXlXN8LVdAo2zF5EWQAtVXSIi9YHFwNXAjcB+VR0jIqOARqp6j9e+acB6YACQAywChqnqan/nzMrK0uzs7EDXZowxxk1EFqtqlq/tAbtxVHWXqi5xvy4A1gCtgCHAq+5qr+L6APDWC9ioqptVtQh4x72fMcaYahRSn72IZAI9gAVAM1XdBa4PBKCpwy6tgO0e73PcZcYYY6pR0MleRE4GJgF/VNX8YHdzKHPsNxKRESKSLSLZeXl5wYZljDEmCEElexGphSvRv6mqH7qL97j788v69XMdds0B2ni8bw3sdDqHqo5X1SxVzcrI8HmPwRhjTBgCJnsREWACsEZVn/LY9Ckw3P16OPCJw+6LgI4i0l5E0oGh7v2MMcZUo2Ba9hcCNwCXishS959BwBhggIhswDXaZgyAiLQUkSkAqloM3AlMxXVj9z1VXVUF12GMMcaPgOPsVXUuzn3vAP0d6u8EBnm8nwJMCTdAY4wxkbMnaI0xpoqoKpMW53CkqCTWoViyN8aYaCkpVYqKS8vfz9+8j7+8v4yHJvt9jrRaWLI3xpgoeHXeVjr8bQo/+vvn5WWHjhYDkJt/LFZhlbNkb4wxUfDAp/E99sSSvTHGpABL9sYYU8VmrNnDe4u2B65YhSzZG2NMNbh70vKYnj/gOHtjjDG+Ldi8jyYnp8c6jIAs2RtjUs6x4hKKS5R6tcNPgUXFpRSXlvKL8d+Gtf9/pq9nWc4BJt7UK+wYQmHJ3hiTcq5+YR5rduWzdcyV5WWqyua9h+mQcXKl+qrKo1PW0KXlKfzp3WW8d3tvHpm8mmU5B8OO4ZmZG8LeNxzWZ2+MSTlrdrlmaV+4ZX952UtzNtP/ya9Y4ZDAc344wktztvCnd5cB8OW63IgSvbcHPllJ/ydnR+14TizZG2NS1uJtJ5bNXrLtAAC/+m/lbhnv1VvHzt7k97i3vroopDhenb+NTXmHQ9onVNaNY4xJWUXFpdwwYQH3DOxcXpbvfuq1zBNT1/G/5Y7LcPg0Y00umaMmRyXGaLFkb4xJGUeKSjhcdCKZr951kDkb9pJ/dCUtTqlTqf7ibft5/suN1RlilbFuHGNMwnn9221kjppcoRsmGNe9OJ+sh2c4bluWc6BS2c/Hzg8nvLhkyd4Yk3Ae+p9rFsmfj53ns87R4yXk5h+tULZiR8Wbqsc8ZqjcdbBi3WRjyd4Yk3SmrdpN59Ff0OvRmeVl6n2XFZi9Ls9x/8xRk5m6aneVxRcLluyNMXFp7oa9fLEyvIT7gsNomfb3hrZg3jsLvw/r3PEq4A1aEXkZGAzkqmpXd9m7QCd3lYbAAVXt7rDvVqAAKAGKVTUrKlEbY5Le9RMWAFR48KlMUUlppTLjXzAt+4nAQM8CVf2FqnZ3J/hJwId+9u/nrmuJ3hhTbt3uAk67dzLb9xdGfKy9h47x4ZKcKESVvAIme1X9GtjvtE1EBLgOeDvKcRljkty7i7ZTqjBpSQ5frXfuOw/Wba9l8+f3lrGn7IasR//8pU/M5uCR4xEdPxlE2md/EbBHVX1N8qDANBFZLCIjIjyXMSYBHSw8ztrd+ZXKRVw/n56xgeEvL2TngSPl21bvrFzfn7Jl/447dO9s3nuYN77dFtLxklGkD1UNw3+r/kJV3SkiTYHpIrLW/U2hEveHwQiAtm3bRhiWMSZe/HTsN2zOO1yp71286r00ZzODzmpBeloNhrzwTUTnzCsIbc1X71iSUdgtexGpCfwMeNdXHVXd6f6ZC3wE+JzLU1XHq2qWqmZlZGSEG5YxJobu+WA5T89YX6Fss485X8Qrw77yzVauHTefBVv2RRzHziQfMx+OSLpxfgysVVXHuyIiUk9E6pe9Bi4DVkZwPmNMnFFV/v3FWrbsdSX0d7O38/SM4KbuFe9s7/bJ0tDmoakYT9i7VrLvcFH0DhYHAiZ7EXkbmA90EpEcEbnFvWkoXl04ItJSRMoGszYD5orIMmAhMFlVv4he6MaYWNu+/whjZ2/ilomhzfIIlVv2wdicd4iNuYd8bt+QW1A+fXGklkdxCmN/Skuj+AnlR8A+e1Ud5qP8RoeyncAg9+vNQLcI4zPGxDHFlaiKvRLW5rxDrN9TwMCuLXzu+/2+0IdcXvrkV3633zwxO+RjxtqkahoyarNeGmOiriwpOz0QtSLnINt/KOTzMJ+OrQrV07Z2Vl3DQi3ZG2PKqSq5Bcdo5jDdbzSs2ZXPVc/PrZJjG/9sbhxjTLk3FnzPeY/ODHmcuwbZNs4NYkhkNG+ymhMs2RuTwr5an1ehG2H+pr0A5aNrwNXaf3bmhgplZSRORqiXBviECHQzuKQ0+efasW4cY1LQku9/4Gf/55oL/oIOTXjrtvN91l24ZT9PTV/PU9PX8+MzmvHV+lw2PDKIp2esp8C9hN/2/Uf4dnPl8fGfLgt9GOXqEEfTlJZqxHPRHykqiWj/RGAte2NSxPf7Cnl6xnpUlfezT4wA2ZTneygjVBxpM2PNHo6XuN4/PWMDE+ZuKd82dHzlhbp///Z3kYbNS19v5uhx38n4iWnrIj5HdYrVPD2W7I1JETdNXMjTMzaQ88MRn3Xisb/8kSlreGam7we1/s9h7npv8XJd63YX0O0f03gve3u1n9uSvTFJJq/gmOOqTGVL8H29IY/lDuutegrUx/3d96Gt/VrGKa5gHDpaHLhSAli/pwAg4lk+w2HJ3pgksm3fYc59ZAbjv97ss859H61kVQijbZzy/lsLkmsVp1Rgyd6YJFLWRRNpy3Hf4SK27XOewAzg/cW2UEiisWRvTIpTdXUvXPCvmew75Jr8a/THK7nk8dmu7TGMzUSPDb00xvDiV5vZefBo3E4NHM6kaaHs72sGzmRiLXtj4lhhUTG7Dx7l1lcX0f2f04Ler7gk+PZ4bsEx/rfc93j4aKbBcL8lxMtomqiJwfVYy96YODVlxS5+++YSx223vrqIVTvzmX9vf8ftC7eeWDY6r+AYjeul+02YRcXOT5DujmJLf+EWx6WsU0osv0BYy96YOPX3j32v9TNjTS67Dh4lc9Rkdnis3eqdS/YdOsa5j8zgsS/WhhXD7a9Hb8rg616cH7VjmdBZsjcmxrK37ue9RcE9ZPP9vkIe+KTih0DZ+6PHSyo87aqq/FDouuE6Y82esGJblnMwqmPCw23YVnWLONzx/4nEunGMiYHiktLyxHzNOFeLt1ubhnRqXt/vfpc+ObvSQiEz1uSSOWpypbovfLmRy89sHnGsL/oZs58sNvlYJzeZWMvemGo2ZcUuznpwGp1Hf1FhzpfCosBPiXonen+emLae/KOueVhSIZklgmpagdBRMGvQviwiuSKy0qPsQRHZISJL3X8G+dh3oIisE5GNIjIqmoEbk0gOFh7n7YWup05/++YSjviY2OtAYRFnjP4iajczj3ncePXs24+V5O8s8e/1+VuB4Of/j6ZgunEmAs8Dr3mV/0dVn/C1k4ikAS8AA4AcYJGIfKqqq8OM1ZiEs31/IRc99mX5+7NbN6iw3bOr+I1vv6dH24YcOV7CdS/Op0m99MgDiLPseszP7JWpIJjFW6pKwJa9qn4NhNPM6AVsVNXNqloEvAMMCeM4xiSUklItv+HnmegBrnzW95J8k5bkUFwS3UU0fvnfBVE9XqQW2PDLmImkz/5OEVnu7uZp5LC9FeA5xCDHXWZM0jpQWESHv03xOxFZsPYdLopCRAbg4+92xDqEmAs32Y8FOgDdgV3Akw51nAZL+fxSKSIjRCRbRLLz8qp/+k9joqHsa/oHQU4UFme9LHEr0pGX6/f4X6AlFYSV7FV1j6qWqGop8BKuLhtvOUAbj/etAZ/PZKvqeFXNUtWsjIyMcMIyJuZSYLh2RML9/divNXJhJXsRaeHx9qeA06N+i4COItJeRNKBocCn4ZzPmESzJz+8aQZCGVppEs+2fYUA7PCzWlhVCWbo5dvAfKCTiOSIyC3AYyKyQkSWA/2AP7nrthSRKQCqWgzcCUwF1gDvqeqqKroOY6rU4WPFPD9rAyUBknHZkLr8o8XkhTHy4rlZG8OKzySWWNyPCTj0UlWHORRP8FF3JzDI4/0UYErY0RkTY8u2HyCthvB+9nZenb+NNo3rMqR7xXEGH3+3g9yCo4y4uEOF8rKpCvzxfkw/VotRx7vkn4C46tl0Ccb4MeSFbwC45pzWQMWHlMr88d2lAIy4uEPIfdKp1mmTAtPGxy1L9saEKOeHQuqm12RP/lFW5BwsL9936BhfrNwdw8hMoojFh54le2NC1OffX1IrTTjutUDIOQ/PCPlYqTb+20YrxY5NhGZMKNzJyjvRO1YNIrE9PHlNhAGZRCQxuAthyd4Ytz35R5m9Ltdxm3U1m0Rnyd4Yt6tf+IYbX1lUvljHml35MY7ImOixZG+S3s4DR3h6xvqAqxHtcq+3OvzlhUBsZyg0yS0WN2gt2ZukVFqq7Droekpx5JtLeHrGhvL5UX750reOKzuZ+CU2ZjNiluxNUhr39SZ6/2sWG3MLOFrkmkO97OnWeZv2xTI0E4ZUWCO2qlmyN0npsS/WAfCP/1VcK2f7/sJYhGNMBbH4nmLj7E1S23HgCLVquNo0eQXH2HrMkr2JvVh8T7GWvUlYq3YeZOzsTQHrHXev/nTDhIVVHZKpItZnHzlL9iZhXfnsXP79xVq/dTbnHaY0zP7eOesrL6ITi4Wik8nEeVtjHUJc8Pzoemr6+mo5pyV7k1JCaSD+d+6WsPYrM2Fu5EsTGpdk/pAoLKqeRdgt2ZuU8vu3v6u2c72XHdzShMZUB0v2Ju5ljprMXe8vC3t/z/5epymKAbK37vd7DBv5ZxKdJXsTl+54awnjvjpx8zXYBbzHzt4U1gNT14ybX+H9gSAWHjHVq+whuWRgo3GMcZu8fBdjPvd/89XJ2NnRWdbvyWkVb5rZYJDYsxlCIxPMGrQvi0iuiKz0KHtcRNaKyHIR+UhEGvrYd6t7rdqlIpIdxbhNEjp0rJjdB/0v1F1w9Dhb9x4O6bhbQqwPUGL9NnFn8vJdsQ4hamLRdgimZT8RGOhVNh3oqqpnA+uBe/3s309Vu6tqVnghmmR19HgJK3ecWOnpqufmcv6/Zvrd5+dj59H3idk+E36oKXrVzoPc/UH49wOMSRTBLDj+tYhkepVN83j7LXBNlOMySW72ulx+9/Z3FBwtZtF9Pyajfu2gWuBlk5mt31NQXpY5ajL//XUWGfVrU3C0OOgYtuw9zJXPznXcZg17U5Vi8ZBYNKZLuBl418c2BaaJiAIvqup4XwcRkRHACIC2bdtGISwTz258ZVH568KiYqB2SPvX8PrPcutrofcS7jsU+hTG9iFgElVEN2hF5D6gGHjTR5ULVbUncAVwh4hc7OtYqjpeVbNUNSsjIyOSsEyS8hyNEY2GUWkIiTsWy8gZE01ht+xFZDgwGOivPuYfVdWd7p+5IvIR0Av4OtxzmuQ08o0lXNTx1ID1ev9rVvnraCT7//h9TN2a8KbquL7NVq+wkr2IDATuAS5RVcdpBEWkHlBDVQvcry8D/hl2pCahHCsuQRXq1EoLWHf1rnxW+1gCMOcH51kqb54Y+eCu+Zt9z2uf80PFMd02J46Jpj351b8KWjBDL98G5gOdRCRHRG4BngfqA9PdwyrHueu2FJEp7l2bAXNFZBmwEJisql9UyVWYuNPv8dl0Hv0Fe8PoF/fU599fRimi0MzZsDcm5zWmqgQzGmeYQ/EEH3V3AoPcrzcD3SKKziSUo8dL6Dz6C/7zi27sdI+Xz3p4BlvHXFleJ3PUZP42qHOsQgyb9dmbRGdP0Jqo+Z17krE/ves8bn3/YdcUBI9OCf3JWGNMZCzZm6iZvnqP3+13vrWkmiIxxnizZG/CMnXVbjJHTWbwc3OC3icZFvp+6LPVgSsZE4dsDVoTkkPHiunz71kcKDwOwModzqNoPKlq0iwrd7iaFpowJtqsZW9CcudbS8oTfbCSYbbCJPmsMinMkr0Jym2vZTN11e4KE5cF69V5W/lybW4VRFV9bJoEk+gs2Ru/Plm6g4Kjx5m+eg+3v77YsU6gG7MAN01cFLCOMabqWLI3Pq3Zlc8f3lnK3R8s91vvtteyKTgaWteOMaZ6WbI3PpWtev/5yt3lZceOO6/hWupcnDTezd4e6xCMiYgle1PJ9NV7+PN7Sx23FRyr/gmcjDGRs2RvANcMkAu37Adc3TIfLtkR0v42UZgx8c2SvQHgmZkbuO7F+RXKbLihMcnDkn2KOlh4vMJiIE6KioPviN+2z3kqYmNMfLBkn6IueeLLCouBOBk6/tugjzfkhW8iDckYU4Us2aeoUJ+CNcYkNkv2KWL/4SIufXI2m/IOxToUY0wMWLJPARO/2ULPh6azOe8wL3292W/dNxdsq6aojDHVyZJ9Cnjwf8FPy3vfRyujfn4b1WNM7AWzBu3LIpIrIis9yhqLyHQR2eD+2cjHvgNFZJ2IbBSRUdEM3ITnyPHKU/T2fGh6lZ7TJhEzJvaCadlPBAZ6lY0CZqpqR2Cm+30FIpIGvABcAXQBholIl4iiNRH7ZOlOpqzYVaGsbLlAY0zyCpjsVfVrYL9X8RDgVffrV4GrHXbtBWxU1c2qWgS8497PVKPF236oVDZ7XWJPN2yMCV24ffbNVHUXgPtnU4c6rQDP2aNy3GWORGSEiGSLSHZeXl6YYRlvPx87r1LZroNHy18/OiX+FhZR6/cxJuqq8gat0205n/+LVXW8qmapalZGRkYVhpV8/vzuUl6dtzXo+nM27C1/PT7A6JxYWL0r8FKHxpjQhJvs94hICwD3T6d+gRygjcf71sDOMM9n/Pjwux088OmqWIfhU6ijccbO3lQ1gRiTwsJN9p8Cw92vhwOfONRZBHQUkfYikg4Mde9nouDLtbm88OXGWIdRJT5bvitwJWNMSIIZevk2MB/oJCI5InILMAYYICIbgAHu94hISxGZAqCqxcCdwFRgDfCeqsZv8zNBHCsu4bpx87lp4iIen7quwrbXv91Gzg+uCcm+2biX3PyjToeodtYFb0zs1QxUQVWH+djU36HuTmCQx/spwJSwozOV3P/xKhZu9R4c5TL6Y9ejEI9dczZ3f7Cclg3qVGdoxpg4Zk/QJphlOQcC1vl0qevWyM6D8dGyN8bEniX7BHegMP4fiCoutX4cY2LNkn0cKyoupdcjM8j20W0D0P2fVTvVgTEmOViyj2MD/vMVuQXHuGbcieUCdwdx09UmHjPGeAt4g9ZUv72HjpH18IwKZaWlylPT19uiI8aYsFiyj0PeiR7gnknLeX9xTgyiMcYkA+vGSRCW6I0xkbBkHydumLCAp6atC1zRGGPCYMk+xlSVRyavZs6GvTw7KzmnPzDGxJ4l+xibOG8rL83ZEtVjes5qaYwxYMk+5lbtrDidb26BPfVqjIk+S/YxcOhYMVe/8A3fbNxbaZKwXo/MjE1QxpikZsm+mhWXlNL1gaks3X6AmycuinU4xpgUYePsq9H+w0X8b5mt32KMqX6W7KvRba9lV1gA3KY1MMZUF0v21WD7/kIueuxLx23qe1leY4yJGuuzrwa+Er0xxlSXsJO9iHQSkaUef/JF5I9edfqKyEGPOvdHHLExxpiQhd2No6rrgO4AIpIG7AA+cqg6R1UHh3ueRJY5ajKXdm4a6zCMMSZq3Tj9gU2qui1Kx0t4maMmAzBrba7felv3Hq6OcIwxKS5ayX4o8LaPbb1FZJmIfC4iZ0bpfElBEJZ8fyDWYRhjUkDEyV5E0oGfAO87bF4CtFPVbsBzwMd+jjNCRLJFJDsvLy/SsGLmyWnrWJFzMNZhGGNMBdFo2V8BLFHVPd4bVDVfVQ+5X08BaonIqU4HUdXxqpqlqlkZGRlRCKv6qSrPzdrIVc/PjXUoxhhTQTSS/TB8dOGISHMR16NDItLLfb59UThnXPKe5yYQe6jKGFNdInqoSkTqAgOA2z3KfgOgquOAa4CRIlIMHAGGqoaaEuPfyh0HadO4LifXtmfUjDHxKaKWvaoWqmoTVT3oUTbOnehR1edV9UxV7aaq56vqvEgDjkeDn5vLL1/6liT8HDOmguG928U6BBMma4pGIOeHQh6f6lpKcNXOfA4XlQCu7hnL+yYZtWlcN9YhmDDZdAkReOCTVXyy9MQsls/N3ABYok9257VvHOsQombmXy4JqX7XVg2qKBJT1SzZR8A7px8tLolJHKZ6VcWN9Z5tG1bLPgBtPVrnp55cO+j9Rg/uwvmnNQnrnCb2LNlHUagt+sIi+3CIV3XT06r1fC0anhTyPuFOxeH5YdXgpFpB79eqYR33T9+xnpvZKKyYTNWzZB8mVa00FYL13iSP3/fv6HObEB9jZgd0ae5YvvahgVE7x097tKpU5u+DMNwuzL9feYbPbZ2a1Q/voKYCS/Zh8l4oHOBg4fEYRGLixVle/dkPXtUl5GM4JVdv57RztZ47NXdOgrVr+v5v/fOerf1+VN0/uGLM9WqncfmZzQLGFInnf9mDWy86zed2ex4lOizZh+mYQ//8N5v2xiASE2092jb020L1lXwu61IxKd54Yfugz1l2yH6dmzLu+p68/5vePHBVF9689TzG/qpnhbrvjDi/Uuvds2UsfrLjzX0y/cbRq31jPv/DRTzg/qBS9d9aPy2jnt/jBaPsPkDDusF3KZnQWbIP0g0TFpA5ajJTV+2OdSjGgVNrtlcYo2aW3X8Zb992vt8VxDxz6bQ/XQzApJEXRKUFqqoM7NqCczMbc9OF7bnw9FO54qwWbB1zZXmdWmk1qFPL1ZXyyo3nMu76nrRoEFyf/5ktG5R/GIy7vqdjnTNanELNGoEvZtz15zD47JYV4/dRz5+yMz0ztEd5mfcHp4mcJfsgzdngarXf/vpiwLm1U1pqvfaxcm1W60plngnj23v7B3WcBnVrUadWms/W7L1XdC7vs2/XpC4/alafrWOuLO9aAejY9GTG/OwswDW0cf69l3LaqfXK388bdWml4/prjZdpfkqdSmX9OjdlYNcWlT6cPr3zQp/94H8bdAYN69aib6fwbvA+cW03Lv5RBv3PaFp+07aM04OFA7s631so43TtWSl6o3fB5qqbTcaSfQA5PxRy4ysLg6prqb5qTBp5gWP5wDNPJJHrstr4PUbzBnW4qKPjHHyOerRpCMCE4Vm0bnSi1Xz5mc3LW/D/HNLVcd/LzmzG0F5tAeiQcXKlVndLh9Estdwt6Rp+kv6suy5h2f2XBRX/2a0b+uwHH9ClGUvvv6z824EjxzhcZd3aNOS1m3tRK61GwN97KDy/TJT9/vzpEebQ03j2xoLvq+zYluwDeHzqOmavS9wpl5NBz7YNme7uLhnSvSUnuZPUmS1PKa9zduuGlfbr45XcfSW3Ru6+4vS0E/8dLjj9VJaMHkD/M3x3J0Rreoxx15/D6MFduKVPe7+t4LrpNWkQQb/2h791/tAMpKxPvV2Tyk/PerfKw/mNlP3+L+jg/GHsK6nbXFShsWTvx/ivN1V4QjYga9pHTVk3CLgSSkd3d8kzQ3uw5qGBfPa7PtzR73Sf+3/91350bn5KhTJfbea0Gq7/Bj/pXrH/uXG99Arv+3duStvGdcsTnPdf9+XubxpXntUSb/dc0ZmTa9esNEZ965grGdi1OY3qpTN6cBdqpVXdf8mebQN3jTgl9JsuzGTeqEs5o8UpDntUdEodV+KeeNO5Aev2Of1Usto1Kv99ptUQ6jsk8Ad/UnnNo7JuMRM8+2j04fMVu3h0ylqf253yuuX66Bnaq63fr/Kej+079TjUSa+cNLu1aci01ZWWXSCtBiy9f0DAluKDPzmTGjWkvCXqfVO47APJyeVnNufyf/jvu46lTs3qM9X97QkgPc31S02vWQMRcex68nbxjzJ46rpufLhkB5f8KPCaFG/cel5QsdWuWfkb2ay7+nLDhAVB7Z9IqnIyRWvZ+zDyzSWxDiGhPXVdNwBu7dPe7wMzAFee1aL8dddWp1Toiw/kpV9nMfuuvpXKm9avfDNz5CUdyl97x9Swbjo1fbSqy75BlE0t8NDVXbl/cBd6RzB1wDNDu4c1Wihcb/lJrEtGD+DjOy6sUPbTHq0Z2bcDf7msU9DneO3mXjQ5uTa3XXwaIsKwXm2ZMDwr7JgB2jQO/cli48xa9mE64jDVQapOcfz1X/tReLyYgU/PAVyP0/+sZ2t+1tM1QubNBa516If1akNmk3r863PXN6YbL8gkvWYNbr/4NCav2AXAZ7+7KKRzD3AYouerdV3D4w7grRedxmkZ9bh5YnbAJ2KH9WrLMI9vGafUqcXNfYIfQ+9kSPdWDOke+AGqaLngdN83p727q8DVor9nYOeIzvkvj664MrPv6kvfJ2b73Gf0VV24/5OV1K2VxsL7+lM3PfQUNbx3O16dvy3k/ZKdJXsvG3MLGP7yooD1xn21qVJZqo28HPurnpxavzZt3f28zwztzh/eWVqpXtnvpYYIt1/Sgdsv6cC3m/fRs20j0t1dIR0y6kU8smPSyAvYtu+w3zo92jakaX1XC71fp6bcfvFpXHNO5WGbpmpkBuhrvy6rTfm/A6dvZ8Hw9Q0t1Vmy96Cq/N/sTew4cMRvvRsmLGDt7oJK5ak06+UzQ7tzhUf3C1BhrLmnsucPPIcVes+eOPMvfSOO6Zx2jXzGUOaj357orhAR7h3kv4spUfj6Uvm/O/tQcCz203hMGtmbgqPF1XKuId1bMmHuFsA1FcOdb31XLeeNd/YR6OEv7y3jwyU7Atabs2Ev+w8XVSpPhl6c3196YoTLM0O7l3dfeD5t+bOeoXVBlGpZso9SkMantl6Li5zVuoHPIY3RsuLBy1j+oP/x/+e0a1z+ENfV3SuPVgrkq7/2Dare1jFXVhiG6/2EbyqLdA3arUABUAIUq2qW13YBngEGAYXAjaoal3c+848e58PvAif6ZPb2befTu0MTnp21kau7t2RI91ZceVYLBnZtXmF0xU+6Of8HKuv79e7P/mmPVny+Yje3e9wgNdFV9qXJezK26lC/Tmhj/5+4tpvjcEp/2jUJPNTyr5cHfzM5FUWjG6efqvqaAewKoKP7z3nAWPfPuLNs+4FYhxBzvTu4ulY8b3DWTKtRaRidr8fs66bXdLw52rBuOu/9pncUIw3eKzedy7yNyT9BXSJ9q6yZVoOGdSvfFI7Ehac38fvcRaKoyr/Gqu6zHwK8pq5hKt+KSEMRaaGqu6r4vCFZuzufGyYENyVCsvlu9ADunrScTbmHYh1KlejXqSn9wpwDJiGlYFfZE9d2Y/DZLQJXTHGR9tkrME1EFovICIftrYDtHu9z3GVxpWzIYDKYe08/v9ufG9aDh64+MadLo3rpvPTrLGY5jFU3JhFk1K/tf54fL2kpevMo0mR/oar2xNVdc4eIXOy13em36vhNRURGiEi2iGTn5VXfXDTb9xdW27mqQ+tGdSuMj35umGva2MeuOZtb+rRn0Fkt6BLEY+/GJLLHrjmbD2LUdRivIurGUdWd7p+5IvIR0Av42qNKDuA5eLo14DjZjKqOB8YDZGVlVVsP5JHjiTdc8uGru/L3j1fSqG4tfnBYHWtk3w6M7HviZuhVXjdUyxaqfvoX3asyTFONOrtXrUqFLiunJ6a9RXM2zmQRdrIXkXpADVUtcL++DPinV7VPgTtF5B1cN2YPxlt/fSLd2Cpz/fntuP78dgC8+NUmpq7azZLvDwS9v4j4fMrUJKaOzeqz8h+Xp8RMkGUPZl3yowzmbNhbabhpIttbcKzKjh3Jv4xmwEfuGetqAm+p6hci8hsAVR0HTME17HIjrqGXN0UWbvT5W5EoHnVrXXFoXdkTqQ99tjrg04mRevGGc2jmsICGiQ+pkOg93dKnPVf3aFU+Z1EyWLBlf5UdO+x/Haq6GejmUD7O47UCd4R7jlT2j5+cyQOfrip//+rNvWjXuC6tGjlPDDV6cOiLW4fq8hAmKDOmqolIWIk+NW/PpvB0CR3vm1Jpgqt40a11A4ZfkMnwCzI5VlxCDZEqnefcmGT275+fRdvG9Rj20rexDiWmUjLZf7J0B8dLlNfmb+MX58bXjZz5915aYa1Rp7m8jTHB+8W58degi4WUbC5+vmJ3+etY3KD1dXP0jn4daNHgpKAWnzbGRC6VVrxKyWQfC043z3plnli8Yt3DA7krhIUijDHh8WxLzbqrL3cPTI3/dynZjeOpulr2V3Vrwe6DR/nnENfTq9/e258GJ9XiZ2PnsWZXPmki1qI3xo/ubRpyaefoP0dQq0ZqtHkt2VfB0Mvsv/+YJdt+YMTri8vL+nZqWmE0S/MGrn75N289j9U7823BBWMC8F46MdpaNTwp4FoWicySfZRz/TXntObUk2szoEszHr/mbC7v2pzc/KOc3rS+Y/3G9dLp07Fq5xs3xgR2Rdfm5BYc49Nljg/5J7yUb04eKy6N+Bh39DsxNcET17oePRARrs1qwyl1avlM9MaY+JLMPakp07LPzT/K/sIiOjevOAnYdS/Oj+i4w3u346+Xd2brvkJqpehsesYkEtcC8ye+0ifaU/ThSplk33vMLEpKla1jrozKX+6pJ9dm0sje5SvovPDLngH2MMbEUo+2DfnOzxxSIsn9dG1KJPvDx4opKY08wd90YSbfbNzLLX3a0/+MZkk1J4cxye7NW8/jQOFx+j4+O9ahxERKJPszH5ha4X3+kdBXuW/V8CQeuCq0dTONMfGjbnpN6qbX5KT0NIqOVL5Xp0pSD39OyRu08zfvC6l+11an8KWt5GRMUpg08oIK7yWpO29OSImWvac5G0JbBevac1rz+LWVJvc0xiSo05ue7HNbMqf9pE/2477aVOF9KAuLv3d7b872mj/eGGMSUdIn+zGfrw1rP1vJyRiTTJI+2Yfqrst+FOsQjDGxksT9OJbs3cqmOPjNJR1snhpjTNKJZMHxNsBrQHOgFBivqs941ekLfAJscRd9qKrei5LH3MNXdy1fwNsYk7qSeWROJC37YuAvqrpEROoDi0Vkuqqu9qo3R1UHR3CeKmV988aYVBB2f4Wq7lLVJe7XBcAaoFW0AqsOmx8dFOsQjDGmWkSlc1pEMoEewAKHzb1FZJmIfC4i1foI6sHC4z63TRp5ATVs4jJjjIdknhQt4hu0InIyMAn4o6rme21eArRT1UMiMgj4GOjo4zgjgBEAbdtGZ4HgA0eKHMtn/uUSOmT4frDCGJM6ymZISPbGX0QtexGphSvRv6mqH3pvV9V8VT3kfj0FqCUijit1qOp4Vc1S1ayMjIxIwipX7DD52ejBXSzRG2PK/eq8dvzqvLbceenpdoPWibhmDJoArFHVp3zUaQ7sUVUVkV64PlxCm5gmDNv2HeaSx2dXWq/yxRvOqbA0oDHGnJSexiM/PSvWYVS5SLpxLgRuAFaIyFJ32d+AtgCqOg64BhgpIsXAEWCoatUv8T1vk+vzZNba3ArlluiNMf4k8aSX4Sd7VZ1LgOfNVPV54PlwzxGq/87ZzMEjx3lu1sZK25bdf1l1hWGMSVBJnOuT6wnahyevcSy/Z2BnGtStVc3RGGMSTTK37FNiXoCRfTsErmSMMUks6ZP9ktEDYh2CMcbEXNIn+8b10mMdgjHGxFxSJ/tXbjw31iEYY0xcSOpk389rnL0xxviTzA9VJW2yb39qvViHYIwxcSNpk/2nd14Y6xCMMQnGhl4mmD/+uCP169i4emOMKZOUyT6ziXXhGGNCZy37BNO2Sd1Yh2CMSUC1a6bFOoQqk5TJvmvLBrEOwRgTx8b+qqfj0Oy7B3YC4NlhPSrUbVIvnWeGdueNW85j9l19AWhXBY3Kzs3rR/2YZaQaJqEMWVZWlmZnZ8c6DGOMSRgislhVs3xtT8qWvTHGmIos2RtjTAqwZG+MMSnAkr0xxqQAS/bGGJMCLNkbY0wKsGRvjDEpwJK9McakgLh8qEpE8oBtYe5+KrA3iuEkArvm5Jdq1wt2zaFqp6oZvjbGZbKPhIhk+3uKLBnZNSe/VLtesGuONuvGMcaYFGDJ3hhjUkAyJvvxsQ4gBuyak1+qXS/YNUdV0vXZG2OMqSwZW/bGGGO8JE2yF5GBIrJORDaKyKhYxxMJEWkjIl+KyBoRWSUif3CXNxaR6SKywf2zkcc+97qvfZ2IXO5Rfo6IrHBve1YkfhdeE5E0EflORD5zv0/2620oIh+IyFr333XvFLjmP7n/Ta8UkbdFpE6yXbOIvCwiuSKy0qMsatcoIrVF5F13+QIRyQwqMFVN+D9AGrAJOA1IB5YBXWIdVwTX0wLo6X5dH1gPdAEeA0a5y0cB/3a/7uK+5tpAe/fvIs29bSHQGxDgc+CKWF+fn+v+M/AW8Jn7fbJf76vAre7X6UDDZL5moBWwBTjJ/f494MZku2bgYqAnsNKjLGrXCPwWGOd+PRR4N6i4Yv2LidIvtzcw1eP9vcC9sY4ritf3CTAAWAe0cJe1ANY5XS8w1f07aQGs9SgfBrwY6+vxcY2tgZnApZxI9sl8vae4E594lSfzNbcCtgONgZrAZ8BlyXjNQKZXso/aNZbVcb+uieshLAkUU7J045T9IyqT4y5LeO6vaD2ABUAzVd0F4P7Z1F3N1/W3cr/2Lo9HTwN3A6UeZcl8vacBecAr7q6r/4pIPZL4mlV1B/AE8D2wCzioqtNI4mv2EM1rLN9HVYuBg0CTQAEkS7J36q9L+GFGInIyMAn4o6rm+6vqUKZ+yuOKiAwGclV1cbC7OJQlzPW61cT1VX+sqvYADuP6eu9Lwl+zu596CK7uipZAPRG53t8uDmUJdc1BCOcaw7r+ZEn2OUAbj/etgZ0xiiUqRKQWrkT/pqp+6C7eIyIt3NtbALnucl/Xn+N+7V0eby4EfiIiW4F3gEtF5A2S93rBFWuOqi5wv/8AV/JP5mv+MbBFVfNU9TjwIXAByX3NZaJ5jeX7iEhNoAGwP1AAyZLsFwEdRaS9iKTjumnxaYxjCpv7rvsEYI2qPuWx6VNguPv1cFx9+WXlQ9136dsDHYGF7q+LBSJyvvuYv/bYJ26o6r2q2lpVM3H93c1S1etJ0usFUNXdwHYR6eQu6g+sJomvGVf3zfkiUtcda39gDcl9zWWieY2ex7oG1/+XwN9sYn0jI4o3RAbhGrWyCbgv1vFEeC19cH0tWw4sdf8ZhKtfbiawwf2zscc+97mvfR0eIxOALGCle9vzBHEjJ8bX3pcTN2iT+nqB7kC2++/5Y6BRClzzP4C17nhfxzUKJamuGXgb1z2J47ha4bdE8xqBOsD7wEZcI3ZOCyYue4LWGGNSQLJ04xhjjPHDkr0xxqQAS/bGGJMCLNkbY0wKsGRvjDEpwJK9McakAEv2xhiTAizZG2NMCvh/rq4O5lo1yDIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count= 1225\n",
      "100.0 %\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAANmUlEQVR4nO3df6zdd13H8efLlhEVkB+9kNkftpiK9g/AcR0YFTFEaOcflYQ/NgzDhaVZQg3+YbIaEiXhH5BgDGHQVGwAY+g/TClQmIaoxMBknRnbytJx2ZBdurBODBJJnIW3f5zv7OFw7z3ntqfr7rvPR3Jyz/f7/dxzP+eTs+fOvufHUlVIkja+n7jcE5AkzYdBl6QmDLokNWHQJakJgy5JTWy+XH94y5YttXPnzsv15yVpQ7r77rsfr6qFlY5dtqDv3LmTkydPXq4/L0kbUpJ/X+2Yp1wkqQmDLklNGHRJasKgS1ITBl2SmjDoktTE1KAnOZrksST3r3I8Sd6fZCnJvUmumf80JUnTzPIM/SPA3jWO7wN2D5cDwIcuflqSpPWaGvSq+gLwnTWG7Ac+ViN3As9NcvW8JihJms08Pim6FXhkbHt52Pfo5MAkBxg9i2fHjh0X/Ad3HvrMBf/u09k33v077Dz0mafVzyfN4/a0tgtZp6fDY+Sp+Pl0dDFzG/9na57m8aJoVti34v8GqaqOVNViVS0uLKz4VQSSpAs0j6AvA9vHtrcBZ+Zwu5KkdZhH0I8DNw7vdnkV8N2q+rHTLZKkS2vqOfQkHwdeA2xJsgz8KfAMgKo6DJwArgOWgO8DN12qyUqSVjc16FV1w5TjBbxtbjOSJF0QPykqSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJamJmYKeZG+S00mWkhxa4fjPJPlUkq8kOZXkpvlPVZK0lqlBT7IJuA3YB+wBbkiyZ2LY24CvVtXLgNcA70ty1ZznKklawyzP0K8Flqrqoap6AjgG7J8YU8CzkwR4FvAd4NxcZypJWtMsQd8KPDK2vTzsG/cB4JeAM8B9wNur6oeTN5TkQJKTSU6ePXv2AqcsSVrJLEHPCvtqYvv1wD3AzwIvBz6Q5Dk/9ktVR6pqsaoWFxYW1jlVSdJaZgn6MrB9bHsbo2fi424Cbq+RJeBh4BfnM0VJ0ixmCfpdwO4ku4YXOq8Hjk+M+SbwWoAkLwJeAjw0z4lKkta2edqAqjqX5CBwB7AJOFpVp5LcMhw/DLwL+EiS+xidorm1qh6/hPOWJE2YGnSAqjoBnJjYd3js+hngdfOdmiRpPfykqCQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJmYKepK9SU4nWUpyaJUxr0lyT5JTSf55vtOUJE2zedqAJJuA24DfBpaBu5Icr6qvjo15LvBBYG9VfTPJCy/RfCVJq5jlGfq1wFJVPVRVTwDHgP0TY94E3F5V3wSoqsfmO01J0jSzBH0r8MjY9vKwb9wvAM9L8k9J7k5y47wmKEmazdRTLkBW2Fcr3M4rgNcCPwl8KcmdVfXgj9xQcgA4ALBjx471z1aStKpZnqEvA9vHtrcBZ1YY87mq+u+qehz4AvCyyRuqqiNVtVhViwsLCxc6Z0nSCmYJ+l3A7iS7klwFXA8cnxjzSeA3kmxO8lPAK4EH5jtVSdJapp5yqapzSQ4CdwCbgKNVdSrJLcPxw1X1QJLPAfcCPwQ+XFX3X8qJS5J+1Czn0KmqE8CJiX2HJ7bfC7x3flOTJK2HnxSVpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJmYKepK9SU4nWUpyaI1xv5LkB0neOL8pSpJmMTXoSTYBtwH7gD3ADUn2rDLuPcAd856kJGm6WZ6hXwssVdVDVfUEcAzYv8K4PwA+ATw2x/lJkmY0S9C3Ao+MbS8P+/5fkq3AG4DDa91QkgNJTiY5efbs2fXOVZK0hlmCnhX21cT2XwC3VtUP1rqhqjpSVYtVtbiwsDDjFCVJs9g8w5hlYPvY9jbgzMSYReBYEoAtwHVJzlXV381jkpKk6WYJ+l3A7iS7gG8B1wNvGh9QVbuevJ7kI8CnjbkkPbWmBr2qziU5yOjdK5uAo1V1Ksktw/E1z5tLkp4aszxDp6pOACcm9q0Y8qr6/YufliRpvfykqCQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJmYKepK9SU4nWUpyaIXjv5fk3uHyxSQvm/9UJUlrmRr0JJuA24B9wB7ghiR7JoY9DPxmVb0UeBdwZN4TlSStbZZn6NcCS1X1UFU9ARwD9o8PqKovVtV/Dpt3AtvmO01J0jSzBH0r8MjY9vKwbzVvBT670oEkB5KcTHLy7Nmzs89SkjTVLEHPCvtqxYHJbzEK+q0rHa+qI1W1WFWLCwsLs89SkjTV5hnGLAPbx7a3AWcmByV5KfBhYF9V/cd8pidJmtUsz9DvAnYn2ZXkKuB64Pj4gCQ7gNuBN1fVg/OfpiRpmqnP0KvqXJKDwB3AJuBoVZ1Kcstw/DDwJ8ALgA8mAThXVYuXbtqSpEmznHKhqk4AJyb2HR67fjNw83ynJklaDz8pKklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU3MFPQke5OcTrKU5NAKx5Pk/cPxe5NcM/+pSpLWMjXoSTYBtwH7gD3ADUn2TAzbB+weLgeAD815npKkKWZ5hn4tsFRVD1XVE8AxYP/EmP3Ax2rkTuC5Sa6e81wlSWtIVa09IHkjsLeqbh623wy8sqoOjo35NPDuqvqXYfvzwK1VdXLitg4wegYP8BLg9AXOewvw+AX+bjeuxXmuxYjrcF7Htfi5qlpY6cDmGX45K+yb/LfALGOoqiPAkRn+5toTSk5W1eLF3k4HrsV5rsWI63DelbYWs5xyWQa2j21vA85cwBhJ0iU0S9DvAnYn2ZXkKuB64PjEmOPAjcO7XV4FfLeqHp3zXCVJa5h6yqWqziU5CNwBbAKOVtWpJLcMxw8DJ4DrgCXg+8BNl27KwBxO2zTiWpznWoy4DuddUWsx9UVRSdLG4CdFJakJgy5JTWy4oE/7GoJuknwjyX1J7klyctj3/CT/kORrw8/njY3/42FtTid5/eWb+cVLcjTJY0nuH9u37vue5BXDGi4NX1Gx0ttsn9ZWWYt3JvnW8Ni4J8l1Y8darkWS7Un+MckDSU4lefuw/4p8XPyYqtowF0Yvyn4deDFwFfAVYM/lntclvs/fALZM7Psz4NBw/RDwnuH6nmFNngnsGtZq0+W+Dxdx318NXAPcfzH3Hfgy8KuMPi/xWWDf5b5vc1qLdwJ/tMLYtmsBXA1cM1x/NvDgcH+vyMfF5GWjPUOf5WsIrgT7gY8O1z8K/O7Y/mNV9T9V9TCjdx1d+9RPbz6q6gvAdyZ2r+u+D19B8Zyq+lKN/in+2NjvbBirrMVq2q5FVT1aVf82XP8e8ACwlSv0cTFpowV9K/DI2PbysK+zAv4+yd3DVycAvKiG9/kPP1847L8S1me9933rcH1yfxcHh284PTp2muGKWIskO4FfBv4VHxfAxgv6TF8x0MyvVdU1jL7R8m1JXr3G2CtxfZ602n3vvCYfAn4eeDnwKPC+YX/7tUjyLOATwB9W1X+tNXSFfa3WYtxGC/oV9xUDVXVm+PkY8LeMTqF8+8lvsxx+PjYMvxLWZ733fXm4Prl/w6uqb1fVD6rqh8Bfcv70Wuu1SPIMRjH/m6q6fdjt44KNF/RZvoagjSQ/neTZT14HXgfcz+g+v2UY9hbgk8P148D1SZ6ZZBej76f/8lM760tuXfd9+M/v7yV51fAuhhvHfmdDm/iK6jcwemxA47UY5v1XwANV9edjh3xcwMZ6l8votQuuY/TK9teBd1zu+Vzi+/piRq/QfwU49eT9BV4AfB742vDz+WO/845hbU6zwV+1Bz7O6FTC/zJ6RvXWC7nvwCKj2H0d+ADDJ6Q30mWVtfhr4D7gXkbhurr7WgC/zujUyL3APcPluiv1cTF58aP/ktTERjvlIklahUGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1IT/wdZCXDVi2th1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fin\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy\n",
    "import math\n",
    "from random import randint\n",
    "\n",
    "def Exp(N,expN):\n",
    "    if N>10:\n",
    "        pi=2*N**0.5/math.log(N)\n",
    "        r=randint(5,N)\n",
    "        if r==N:\n",
    "            r=N-2\n",
    "        s=N-r\n",
    "        if r<=10:\n",
    "            return expN[r],expN\n",
    "        else:\n",
    "            expN[N]=max(pi-expN[r]-expN[s],expN[N])\n",
    "        if s<=10:\n",
    "            return expN[s],expN\n",
    "        else:\n",
    "            expN[N]=max(pi-expN[r]-expN[s],expN[N])\n",
    "        return expN[N],expN\n",
    "    else:\n",
    "        return expN[N],expN\n",
    "\n",
    "n_iter=100\n",
    "Nmax=10000\n",
    "primes=list(sympy.primerange(2,Nmax+1))\n",
    "prime_train=[2,3,5,7]\n",
    "expN=[1]*(Nmax+2)\n",
    "expN[2]=1\n",
    "expN[3]=1\n",
    "expN[4]=1\n",
    "expN[5]=1\n",
    "expN[6]=2\n",
    "expN[7]=1\n",
    "expN[8]=1\n",
    "expN[9]=1\n",
    "expN[10]=2\n",
    "\n",
    "PrimeTrue=[]\n",
    "PrimeCheck=[]\n",
    "countskip=0\n",
    "for n in range(10,Nmax+1):\n",
    "    if (n%2==0 or n%3==0) or (n%5==0 or n%7==0):\n",
    "        continue\n",
    "    else:\n",
    "        countskip+=1\n",
    "    BOOL=True\n",
    "    if n in primes:\n",
    "        PrimeTrue.append(1)\n",
    "    else:\n",
    "        PrimeTrue.append(0)\n",
    "    pi=2*n**0.5/math.log(n)\n",
    "    \n",
    "    BOOL=True\n",
    "    for m in range(n_iter):\n",
    "        r=randint(5,n)\n",
    "        if r==n:\n",
    "            r=n-2\n",
    "        s=n-r\n",
    "        g=math.gcd(s,r)\n",
    "\n",
    "        if BOOL:\n",
    "            if g!=1:\n",
    "                PrimeCheck.append(0)\n",
    "                prob=0\n",
    "                BOOL=False\n",
    "                continue\n",
    "            if (n%2==0 or n%3==0) or (n%5==0 or n%7==0):\n",
    "                PrimeCheck.append(0)\n",
    "                prob=0\n",
    "                BOOL=False\n",
    "                continue\n",
    "        if m==0:\n",
    "            ermax,expN=Exp(r,expN)\n",
    "            esmax,expN=Exp(s,expN)\n",
    "        else:\n",
    "            er,expN=Exp(r,expN)\n",
    "            es,expN=Exp(s,expN)\n",
    "            ermax=max(er,ermax)\n",
    "            esmax=max(es,esmax)\n",
    "        prob=(ermax+esmax)/pi\n",
    "    if BOOL:\n",
    "        if prob >0.5:\n",
    "            PrimeCheck.append(1)\n",
    "        else:\n",
    "            PrimeCheck.append(0)\n",
    "\n",
    "x=[]\n",
    "Hit=[]\n",
    "count=0\n",
    "for n in range(len(PrimeCheck)):\n",
    "    C=PrimeCheck[n]\n",
    "    T=PrimeTrue[n]\n",
    "    if C==T:\n",
    "        count+=1\n",
    "        x.append(n)\n",
    "        Hit.append(1)\n",
    "    else:\n",
    "        x.append(n)\n",
    "        Hit.append(0)\n",
    "\n",
    "print('count=',count)\n",
    "#print(round(count/Nmax*100,2),'%')\n",
    "print(round(count/countskip*100,2),'%')\n",
    "\n",
    "plt.bar(x,Hit)\n",
    "plt.show()\n",
    "plt.plot(expN)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "count=0\n",
    "for n in range(len(PrimeCheck)):\n",
    "    C=PrimeCheck[n]\n",
    "    T=PrimeTrue[n]\n",
    "    if C==T and T==1:\n",
    "        count+=1\n",
    "        x.append(n)\n",
    "        Hit.append(1)\n",
    "    else:\n",
    "        x.append(n)\n",
    "        Hit.append(0)\n",
    "\n",
    "print('count=',count)\n",
    "print(round(count/(len(primes)-4)*100,2),'%')\n",
    "\n",
    "plt.bar(x,Hit)\n",
    "plt.show()\n",
    "\n",
    "print('Fin')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}