NDV.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_info=np.load(\"learned_ndv_estimator/model_training/training_data/rfs_F_infos.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_s=np.load(\"learned_ndv_estimator/model_training/training_data/rfs_f_s.npy\",allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.polynomial import Chebyshev\n",
    "from numpy.polynomial import Hermite\n",
    "from scipy import optimize\n",
    "from scipy.special import comb\n",
    "from estndv import ndvEstimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import NullFormatter\n",
    "from matplotlib import style\n",
    "#plt.rcParams['font.sans-serif']=['simhei']\n",
    "matplotlib.rcParams['text.usetex'] = True\n",
    "from collections import Counter\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimator\n",
    "def Chao(fi):\n",
    "    if fi.shape[0]==1:\n",
    "        return sum(fi[:,1])\n",
    "    if fi[0,0]==1 and fi[1,0]==2:\n",
    "        return sum(fi[:,1])+fi[0,1]*(fi[0,1]-1)/2/(fi[1,1]+1)\n",
    "    else:\n",
    "        return sum(fi[:,1])\n",
    "def GEE(fi,q):\n",
    "    import math\n",
    "    if fi[0,0]==1:\n",
    "        return sum(fi[:,1])+math.sqrt(1/q)*fi[0,1]\n",
    "    else:\n",
    "        return sum(fi[:,1])\n",
    "def Cheby4f0(fi,pmin=10**(-5)):\n",
    "    L=int(0.45*np.log(1/pmin))\n",
    "    ratio=0.5*np.log(1/pmin)\n",
    "    num = sum(fi[:,0]*fi[:,1])\n",
    "    def shift_polynomial(coeffs, shift):\n",
    "        length = len(coeffs)\n",
    "        coeffs_shift = np.zeros(length)\n",
    "        for i in range(length):\n",
    "            for j in range(i+1):\n",
    "                coeffs_shift[j] += coeffs[i] * (shift**(i-j)) * scipy.special.binom(i, j)\n",
    "        return coeffs_shift\n",
    "    if pmin<ratio/num:\n",
    "        cheb = np.polynomial.chebyshev.Chebyshev.basis(L)\n",
    "        cheb_coeffs = np.polynomial.chebyshev.cheb2poly(cheb.coef)\n",
    "        shift = (ratio+num*pmin)/(ratio-num*pmin)\n",
    "        a_coeffs = shift_polynomial(cheb_coeffs, -shift)\n",
    "        g_coeffs = -a_coeffs / a_coeffs[0]\n",
    "        g_coeffs[0] = 0\n",
    "        scaling = 2. / (ratio-num*pmin)\n",
    "        for j in range(1, L+1):\n",
    "            for i in range(1, j+1):\n",
    "                g_coeffs[j] *= (i*scaling)\n",
    "            g_coeffs[j]+=1\n",
    "        esti=0\n",
    "        for [i,j] in fi:\n",
    "            if i > L:\n",
    "                esti+=j\n",
    "            else:\n",
    "                esti+=g_coeffs[i]*j\n",
    "        if esti<sum(fi[:,1]):\n",
    "            esti=sum(fi[:,1])\n",
    "    else:\n",
    "        esti=sum(fi[:,1])\n",
    "    return esti\n",
    "def GT(fi):\n",
    "    t=0.1\n",
    "    return sum(fi[:,1])+ sum((-1)**(fi[:,0]+1)*t**fi[:,0]*fi[:,1])\n",
    "    \n",
    "def Shlosser(fi,q):\n",
    "    Ef=0\n",
    "    Exf=0\n",
    "    if fi[0,0]==1:\n",
    "        for item in fi:\n",
    "            Ef+=(1-q)**item[0]*item[1]\n",
    "            Exf+=q*item[0]*(1-q)**(item[0]-1)*item[1]\n",
    "        if Exf==0 or fi[0,1]==0:\n",
    "            return sum(fi[:,1])\n",
    "        else:\n",
    "            return sum(fi[:,1])+Ef/Exf*fi[0,1]\n",
    "    else:\n",
    "        return sum(fi[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jiajun_li/anaconda3/envs/py37/lib/python3.7/site-packages/ipykernel_launcher.py:84: UserWarning: Creating a tensor from a list of numpy.ndarrays is extremely slow. Please consider converting the list to a single numpy.ndarray with numpy.array() before converting to a tensor. (Triggered internally at ../torch/csrc/utils/tensor_new.cpp:230.)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prepare time 416.5573604106903\n",
      "epoch 0\n",
      "epoch 1\n",
      "epoch 2\n",
      "epoch 3\n",
      "epoch 4\n",
      "epoch 5\n",
      "epoch 6\n",
      "epoch 7\n",
      "epoch 8\n",
      "epoch 9\n",
      "epoch 10\n",
      "epoch 11\n",
      "epoch 12\n",
      "epoch 13\n",
      "epoch 14\n",
      "epoch 15\n",
      "epoch 16\n",
      "epoch 17\n",
      "epoch 18\n",
      "epoch 19\n",
      "epoch 20\n",
      "epoch 21\n",
      "epoch 22\n",
      "epoch 23\n",
      "epoch 24\n",
      "epoch 25\n",
      "epoch 26\n",
      "epoch 27\n",
      "epoch 28\n",
      "epoch 29\n",
      "epoch 30\n",
      "epoch 31\n",
      "epoch 32\n",
      "epoch 33\n",
      "epoch 34\n",
      "epoch 35\n",
      "epoch 36\n",
      "epoch 37\n",
      "epoch 38\n",
      "epoch 39\n",
      "epoch 40\n",
      "epoch 41\n",
      "epoch 42\n",
      "epoch 43\n",
      "epoch 44\n",
      "epoch 45\n",
      "epoch 46\n",
      "epoch 47\n",
      "epoch 48\n",
      "epoch 49\n",
      "epoch 50\n",
      "epoch 51\n",
      "epoch 52\n",
      "epoch 53\n",
      "epoch 54\n",
      "epoch 55\n",
      "epoch 56\n",
      "epoch 57\n",
      "epoch 58\n",
      "epoch 59\n",
      "epoch 60\n",
      "epoch 61\n",
      "epoch 62\n",
      "epoch 63\n",
      "epoch 64\n",
      "epoch 65\n",
      "epoch 66\n",
      "epoch 67\n",
      "epoch 68\n",
      "epoch 69\n",
      "epoch 70\n",
      "epoch 71\n",
      "epoch 72\n",
      "epoch 73\n",
      "epoch 74\n",
      "epoch 75\n",
      "epoch 76\n",
      "epoch 77\n",
      "epoch 78\n",
      "epoch 79\n",
      "epoch 80\n",
      "epoch 81\n",
      "epoch 82\n",
      "epoch 83\n",
      "epoch 84\n",
      "epoch 85\n",
      "epoch 86\n",
      "epoch 87\n",
      "epoch 88\n",
      "epoch 89\n",
      "epoch 90\n",
      "epoch 91\n",
      "epoch 92\n",
      "epoch 93\n",
      "epoch 94\n",
      "epoch 95\n",
      "epoch 96\n",
      "epoch 97\n",
      "epoch 98\n",
      "epoch 99\n",
      "epoch 100\n",
      "epoch 101\n",
      "epoch 102\n",
      "epoch 103\n",
      "epoch 104\n",
      "epoch 105\n",
      "epoch 106\n",
      "epoch 107\n",
      "epoch 108\n",
      "epoch 109\n",
      "epoch 110\n",
      "epoch 111\n",
      "epoch 112\n",
      "epoch 113\n",
      "epoch 114\n",
      "epoch 115\n",
      "epoch 116\n",
      "epoch 117\n",
      "epoch 118\n",
      "epoch 119\n",
      "epoch 120\n",
      "epoch 121\n",
      "epoch 122\n",
      "epoch 123\n",
      "epoch 124\n",
      "epoch 125\n",
      "epoch 126\n",
      "epoch 127\n",
      "epoch 128\n",
      "epoch 129\n",
      "epoch 130\n",
      "epoch 131\n",
      "epoch 132\n",
      "epoch 133\n",
      "epoch 134\n",
      "epoch 135\n",
      "epoch 136\n",
      "epoch 137\n",
      "epoch 138\n",
      "epoch 139\n",
      "epoch 140\n",
      "epoch 141\n",
      "epoch 142\n",
      "epoch 143\n",
      "epoch 144\n",
      "epoch 145\n",
      "epoch 146\n",
      "epoch 147\n",
      "epoch 148\n",
      "epoch 149\n",
      "epoch 150\n",
      "epoch 151\n",
      "epoch 152\n",
      "epoch 153\n",
      "epoch 154\n",
      "epoch 155\n",
      "epoch 156\n",
      "epoch 157\n",
      "epoch 158\n",
      "epoch 159\n",
      "epoch 160\n",
      "epoch 161\n",
      "epoch 162\n",
      "epoch 163\n",
      "epoch 164\n",
      "epoch 165\n",
      "epoch 166\n",
      "epoch 167\n",
      "epoch 168\n",
      "epoch 169\n",
      "epoch 170\n",
      "epoch 171\n",
      "epoch 172\n",
      "epoch 173\n",
      "epoch 174\n",
      "epoch 175\n",
      "epoch 176\n",
      "epoch 177\n",
      "epoch 178\n",
      "epoch 179\n",
      "epoch 180\n",
      "epoch 181\n",
      "epoch 182\n",
      "epoch 183\n",
      "epoch 184\n",
      "epoch 185\n",
      "epoch 186\n",
      "epoch 187\n",
      "epoch 188\n",
      "epoch 189\n",
      "epoch 190\n",
      "epoch 191\n",
      "epoch 192\n",
      "epoch 193\n",
      "epoch 194\n",
      "epoch 195\n",
      "epoch 196\n",
      "epoch 197\n",
      "epoch 198\n",
      "epoch 199\n",
      "epoch 200\n",
      "epoch 201\n",
      "epoch 202\n",
      "epoch 203\n",
      "epoch 204\n",
      "epoch 205\n",
      "epoch 206\n",
      "epoch 207\n",
      "epoch 208\n",
      "epoch 209\n",
      "epoch 210\n",
      "epoch 211\n",
      "epoch 212\n",
      "epoch 213\n",
      "epoch 214\n",
      "epoch 215\n",
      "epoch 216\n",
      "epoch 217\n",
      "epoch 218\n",
      "epoch 219\n",
      "epoch 220\n",
      "epoch 221\n",
      "epoch 222\n",
      "epoch 223\n",
      "epoch 224\n",
      "epoch 225\n",
      "epoch 226\n",
      "epoch 227\n",
      "epoch 228\n",
      "epoch 229\n",
      "epoch 230\n",
      "epoch 231\n",
      "epoch 232\n",
      "epoch 233\n",
      "epoch 234\n",
      "epoch 235\n",
      "epoch 236\n",
      "epoch 237\n",
      "epoch 238\n",
      "epoch 239\n",
      "epoch 240\n",
      "epoch 241\n",
      "epoch 242\n",
      "epoch 243\n",
      "epoch 244\n",
      "epoch 245\n",
      "epoch 246\n",
      "epoch 247\n",
      "epoch 248\n",
      "epoch 249\n",
      "epoch 250\n",
      "epoch 251\n",
      "epoch 252\n",
      "epoch 253\n",
      "epoch 254\n",
      "epoch 255\n",
      "epoch 256\n",
      "epoch 257\n",
      "epoch 258\n",
      "epoch 259\n",
      "epoch 260\n",
      "epoch 261\n",
      "epoch 262\n",
      "epoch 263\n",
      "epoch 264\n",
      "epoch 265\n",
      "epoch 266\n",
      "epoch 267\n",
      "epoch 268\n",
      "epoch 269\n",
      "epoch 270\n",
      "epoch 271\n",
      "epoch 272\n",
      "epoch 273\n",
      "epoch 274\n",
      "epoch 275\n",
      "epoch 276\n",
      "epoch 277\n",
      "epoch 278\n",
      "epoch 279\n",
      "epoch 280\n",
      "epoch 281\n",
      "epoch 282\n",
      "epoch 283\n",
      "epoch 284\n",
      "epoch 285\n",
      "epoch 286\n",
      "epoch 287\n",
      "epoch 288\n",
      "epoch 289\n",
      "epoch 290\n",
      "epoch 291\n",
      "epoch 292\n",
      "epoch 293\n",
      "epoch 294\n",
      "epoch 295\n",
      "epoch 296\n",
      "epoch 297\n",
      "epoch 298\n",
      "epoch 299\n",
      "epoch 300\n",
      "epoch 301\n",
      "epoch 302\n",
      "epoch 303\n",
      "epoch 304\n",
      "epoch 305\n",
      "epoch 306\n",
      "epoch 307\n",
      "epoch 308\n",
      "epoch 309\n",
      "epoch 310\n",
      "epoch 311\n",
      "epoch 312\n",
      "epoch 313\n",
      "epoch 314\n",
      "epoch 315\n",
      "epoch 316\n",
      "epoch 317\n",
      "epoch 318\n",
      "epoch 319\n",
      "epoch 320\n",
      "epoch 321\n",
      "epoch 322\n",
      "epoch 323\n",
      "epoch 324\n",
      "epoch 325\n",
      "epoch 326\n",
      "epoch 327\n",
      "epoch 328\n",
      "epoch 329\n",
      "epoch 330\n",
      "epoch 331\n",
      "epoch 332\n",
      "epoch 333\n",
      "epoch 334\n",
      "epoch 335\n",
      "epoch 336\n",
      "epoch 337\n",
      "epoch 338\n",
      "epoch 339\n",
      "epoch 340\n",
      "epoch 341\n",
      "epoch 342\n",
      "epoch 343\n",
      "epoch 344\n",
      "epoch 345\n",
      "epoch 346\n",
      "epoch 347\n",
      "epoch 348\n",
      "epoch 349\n",
      "epoch 350\n",
      "epoch 351\n",
      "epoch 352\n",
      "epoch 353\n",
      "epoch 354\n",
      "epoch 355\n",
      "epoch 356\n",
      "epoch 357\n",
      "epoch 358\n",
      "epoch 359\n",
      "epoch 360\n",
      "epoch 361\n",
      "epoch 362\n",
      "epoch 363\n",
      "epoch 364\n",
      "epoch 365\n",
      "epoch 366\n",
      "epoch 367\n",
      "epoch 368\n",
      "epoch 369\n",
      "epoch 370\n",
      "epoch 371\n",
      "epoch 372\n",
      "epoch 373\n",
      "epoch 374\n",
      "epoch 375\n",
      "epoch 376\n",
      "epoch 377\n",
      "epoch 378\n",
      "epoch 379\n",
      "epoch 380\n",
      "epoch 381\n",
      "epoch 382\n",
      "epoch 383\n",
      "epoch 384\n",
      "epoch 385\n",
      "epoch 386\n",
      "epoch 387\n",
      "epoch 388\n",
      "epoch 389\n",
      "epoch 390\n",
      "epoch 391\n",
      "epoch 392\n",
      "epoch 393\n",
      "epoch 394\n",
      "epoch 395\n",
      "epoch 396\n",
      "epoch 397\n",
      "epoch 398\n",
      "epoch 399\n",
      "epoch 400\n",
      "epoch 401\n",
      "epoch 402\n",
      "epoch 403\n",
      "epoch 404\n",
      "epoch 405\n",
      "epoch 406\n",
      "epoch 407\n",
      "epoch 408\n",
      "epoch 409\n",
      "epoch 410\n",
      "epoch 411\n",
      "epoch 412\n",
      "epoch 413\n",
      "epoch 414\n",
      "epoch 415\n",
      "epoch 416\n",
      "epoch 417\n",
      "epoch 418\n",
      "epoch 419\n",
      "epoch 420\n",
      "epoch 421\n",
      "epoch 422\n",
      "epoch 423\n",
      "epoch 424\n",
      "epoch 425\n",
      "epoch 426\n",
      "epoch 427\n",
      "epoch 428\n",
      "epoch 429\n",
      "epoch 430\n",
      "epoch 431\n",
      "epoch 432\n",
      "epoch 433\n",
      "epoch 434\n",
      "epoch 435\n",
      "epoch 436\n",
      "epoch 437\n",
      "epoch 438\n",
      "epoch 439\n",
      "epoch 440\n",
      "epoch 441\n",
      "epoch 442\n",
      "epoch 443\n",
      "epoch 444\n",
      "epoch 445\n",
      "epoch 446\n",
      "epoch 447\n",
      "epoch 448\n",
      "epoch 449\n",
      "epoch 450\n",
      "epoch 451\n",
      "epoch 452\n",
      "epoch 453\n",
      "epoch 454\n",
      "epoch 455\n",
      "epoch 456\n",
      "epoch 457\n",
      "epoch 458\n",
      "epoch 459\n",
      "epoch 460\n",
      "epoch 461\n",
      "epoch 462\n",
      "epoch 463\n",
      "epoch 464\n",
      "epoch 465\n",
      "epoch 466\n",
      "epoch 467\n",
      "epoch 468\n",
      "epoch 469\n",
      "epoch 470\n",
      "epoch 471\n",
      "epoch 472\n",
      "epoch 473\n",
      "epoch 474\n",
      "epoch 475\n",
      "epoch 476\n",
      "epoch 477\n",
      "epoch 478\n",
      "epoch 479\n",
      "epoch 480\n",
      "epoch 481\n",
      "epoch 482\n",
      "epoch 483\n",
      "epoch 484\n",
      "epoch 485\n",
      "epoch 486\n",
      "epoch 487\n",
      "epoch 488\n",
      "epoch 489\n",
      "epoch 490\n",
      "epoch 491\n",
      "epoch 492\n",
      "epoch 493\n",
      "epoch 494\n",
      "epoch 495\n",
      "epoch 496\n",
      "epoch 497\n",
      "epoch 498\n",
      "epoch 499\n",
      "epoch 500\n",
      "epoch 501\n",
      "epoch 502\n",
      "epoch 503\n",
      "epoch 504\n",
      "epoch 505\n",
      "epoch 506\n",
      "epoch 507\n",
      "epoch 508\n",
      "epoch 509\n",
      "epoch 510\n",
      "epoch 511\n",
      "epoch 512\n",
      "epoch 513\n",
      "epoch 514\n",
      "epoch 515\n",
      "epoch 516\n",
      "epoch 517\n",
      "epoch 518\n",
      "epoch 519\n",
      "epoch 520\n",
      "epoch 521\n",
      "epoch 522\n",
      "epoch 523\n",
      "epoch 524\n",
      "epoch 525\n",
      "epoch 526\n",
      "epoch 527\n",
      "epoch 528\n",
      "epoch 529\n",
      "epoch 530\n",
      "epoch 531\n",
      "epoch 532\n",
      "epoch 533\n",
      "epoch 534\n",
      "epoch 535\n",
      "epoch 536\n",
      "epoch 537\n",
      "epoch 538\n",
      "epoch 539\n",
      "epoch 540\n",
      "epoch 541\n",
      "epoch 542\n",
      "epoch 543\n",
      "epoch 544\n",
      "epoch 545\n",
      "epoch 546\n",
      "epoch 547\n",
      "epoch 548\n",
      "epoch 549\n",
      "epoch 550\n",
      "epoch 551\n",
      "epoch 552\n",
      "epoch 553\n",
      "epoch 554\n",
      "epoch 555\n",
      "epoch 556\n",
      "epoch 557\n",
      "epoch 558\n",
      "epoch 559\n",
      "epoch 560\n",
      "epoch 561\n",
      "epoch 562\n",
      "epoch 563\n",
      "epoch 564\n",
      "epoch 565\n",
      "epoch 566\n",
      "epoch 567\n",
      "epoch 568\n",
      "epoch 569\n",
      "epoch 570\n",
      "epoch 571\n",
      "epoch 572\n",
      "epoch 573\n",
      "epoch 574\n",
      "epoch 575\n",
      "epoch 576\n",
      "epoch 577\n",
      "epoch 578\n",
      "epoch 579\n",
      "epoch 580\n",
      "epoch 581\n",
      "epoch 582\n",
      "epoch 583\n",
      "epoch 584\n",
      "epoch 585\n",
      "epoch 586\n",
      "epoch 587\n",
      "epoch 588\n",
      "epoch 589\n",
      "epoch 590\n",
      "epoch 591\n",
      "epoch 592\n",
      "epoch 593\n",
      "epoch 594\n",
      "epoch 595\n",
      "epoch 596\n",
      "epoch 597\n",
      "epoch 598\n",
      "epoch 599\n",
      "training time 335.85153675079346\n"
     ]
    }
   ],
   "source": [
    "# Our Model\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torch.nn import Parameter\n",
    "from scipy.special import comb\n",
    "from scipy import interpolate\n",
    "import scipy\n",
    "\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "class PolyNet(nn.Module):\n",
    "    def __init__(self,L=7,H=200):\n",
    "        super(PolyNet,self).__init__()\n",
    "        self.L = L\n",
    "        self.H = H\n",
    "        self.weight = Parameter(torch.Tensor(self.H))\n",
    "        self.weight.requires_grad = True\n",
    "        self.fc1 = nn.Linear(self.L*10, self.L*10//2)\n",
    "        self.fc2 = nn.Linear(self.L*10//2, self.L)\n",
    "        self.fc3 = nn.Linear(self.L, self.L)\n",
    "        self.fc1.requires_grad = True\n",
    "        self.fc2.requires_grad = True\n",
    "        self.fc3.requires_grad = True\n",
    "        self.i = torch.Tensor( [ i+1 for i in range(self.H) ] ).float().to(device)\n",
    "    def set_paras(self):\n",
    "        self.weight.data = torch.Tensor([1.0 for i in range(self.H)])\n",
    "    def forward(self,fi,M_tensor,N,n):\n",
    "        bt_temp= F.relu(self.fc1( fi) )\n",
    "        bt=self.fc2(bt_temp)\n",
    "        f0_loss = torch.einsum('ij,ijk->ik', bt,M_tensor)\n",
    "        f0_loss=((1.0-f0_loss)*self.weight/self.weight.sum()).sum(1)\n",
    "        f0 =    torch.relu(  (bt*fi[:,0:self.L]).sum(1) )\n",
    "        return f0,f0_loss\n",
    "pnet = PolyNet().to(device)\n",
    "pnet.set_paras()\n",
    "for param in pnet.parameters():\n",
    "    param.data = param.data.to(device)\n",
    "\n",
    "optimizer = torch.optim.Adam(pnet.parameters(),lr=0.001)\n",
    "\n",
    "# data prepare\n",
    "\n",
    "fea = torch.tensor([[0.0 for j in range(70)]  for i in range(f_s.shape[0])]).to(device)\n",
    "entropy_plus=[]\n",
    "for i in range(f_s.shape[0]):\n",
    "    f_s_n=sum(f_s[i][:,0]*f_s[i][:,1])\n",
    "    entropy_p=0.0\n",
    "    for item in f_s[i]:\n",
    "        if item[0]>70:\n",
    "            entropy_p+= item[1]*item[0]/f_s_n*np.log(f_s_n/item[0])\n",
    "        else:\n",
    "            fea[i][item[0]-1]=item[1]\n",
    "    entropy_plus.append(entropy_p)\n",
    "entropy_plus = torch.Tensor(entropy_plus).to(device)\n",
    "\n",
    "d = torch.tensor([ sum(f[:,1]) for f in f_s]).to(device)\n",
    "n = torch.tensor([ sum(f[:,0]*f[:,1]) for f in f_s] ).to(device)\n",
    "N = torch.tensor(F_info[:,0]).to(device)\n",
    "\n",
    "#Entropy = torch.tensor( F_info[:,2]).to(device)\n",
    "import time\n",
    "begin1=time.time()\n",
    "M_tensor=[]\n",
    "size = F_info.shape[0]\n",
    "L=7\n",
    "H=200\n",
    "data_n=np.array([sum(te[:,1]) for te in f_s])\n",
    "data_N=F_info[:,0]\n",
    "for s in range(size):\n",
    "    M=[]\n",
    "    for t in range(1,L+1):\n",
    "        M_temp=[]\n",
    "        for j in range(1,H+1):\n",
    "            M_temp.append( comb(data_n[s],t)*(j/(data_N[s]-j))**t)\n",
    "        M.append(M_temp)\n",
    "    M=np.array(M)\n",
    "    M[M<0.001]=0\n",
    "    M_tensor.append(M)    \n",
    "M_tensor = torch.tensor(M_tensor).float().to(device)\n",
    "begin=time.time()\n",
    "print(\"prepare time\", begin-begin1)\n",
    "def lossfun(y,y_pred,d,n,N,l0_loss):\n",
    "    l_log_unreduced = torch.square( torch.log( (y_pred+d)/(y+d)))\n",
    "    l_bound = torch.sqrt(( (1 / (n.true_divide(N) + 1e-10) - 1) / 2 * 0.2554128 + d) / (d + 1e-6))\n",
    "    l_log_bound = torch.square(torch.log(l_bound))\n",
    "    l = torch.mean(torch.abs(l_log_unreduced - l_log_bound) + l_log_bound + 5*torch.log(l0_loss.abs()+1))  \n",
    "    return torch.log(l)\n",
    "\n",
    "batch_size=200\n",
    "train_loss=[]\n",
    "for i in range(600):\n",
    "    print(\"epoch\",i)\n",
    "    for epoch in range(200):\n",
    "        predict,l0_loss=pnet(\n",
    "            fea[epoch*batch_size:((epoch+1)*batch_size)] , \n",
    "            M_tensor[epoch*batch_size:((epoch+1)*batch_size)],\n",
    "            N[epoch*batch_size:((epoch+1)*batch_size)],\n",
    "            n[epoch*batch_size:((epoch+1)*batch_size)])\n",
    "        #print(\"finish predict\")\n",
    "        y=torch.tensor(F_info[:,1]).to(device)-d\n",
    "        loss=lossfun(y[epoch*batch_size:((epoch+1)*batch_size)],predict,\n",
    "                     d[epoch*batch_size:((epoch+1)*batch_size)],\n",
    "                     #torch.tensor(F_info[epoch*batch_size:((epoch+1)*batch_size),1]),\n",
    "                     n[epoch*batch_size:((epoch+1)*batch_size)],\n",
    "                     torch.tensor(F_info[epoch*batch_size:((epoch+1)*batch_size),0]).to(device),\n",
    "                     l0_loss)\n",
    "        #print(sum(l0_loss))\n",
    "        loss.requires_grad_(True)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "    train_loss.append(loss)\n",
    "end=time.time()\n",
    "print(\"training time\",end-begin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEGCAYAAACNaZVuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5b0/8M8zM9kXsrKEJWECyCoSAoIWqiW2VmvdUNp7u7jUoPa23i4/tdrbvnqrVbDX7rXGttpaFwT13lZcalBE0QpJEAFZEwhrSMi+Z5bn98ecMzkzmT2ZTM6cz/v14sXMOWcmz5DhM898n+c8R0gpQURExmKKdQOIiGj0MfyJiAyI4U9EZEAMfyIiA2L4ExEZEMOfiMiALLFuQKjy8vJkUVFRrJtBRKQb1dXV56SU+b726Sb8i4qKUFVVFetmEBHphhCi3t8+ln2IiAyI4U9EZEAMfyIiA2L4ExEZEMOfiMiAGP5ERAYU1fAXQpQpf9b52b9O+bs8mu0gIiJPUQt/IUQZgBuklJUASoQQJT4OKxdC1AKoi1Y7nvmwHkX3bobd4YzWjyAi0p2oneSlhH6lctcqpazxcdhtUspN0WoDAKx77QAAoMfmQKaZVS4iImAUav5CiLsBrPWz26qUhe6O1s/ndcqIiIaKevhLKdcDWCuEyPK1T/mGkKuUiaLXDlZ9iIjcolnz19b56wCUe+0vF0KsVu42A7D6eI5yIUSVEKKqqakpsoYoXX+7k+lPRKSKZs+/DECOcjsLyqCu5htAFQbHBIqV+x6klBVSylIpZWl+vs+F6ULmcLIARESkimb4V8BV0y8HAM3A7hblfg2AG5Xef62fAeER45AMfyIiVTRn+7TB9QHgvX2x5vaQ/dFidzD8iYhUhpn76JQS9c3deOi1/ZD8FkBEBhf34a/GvN0psfbpajz+Th1qm7pi2iYioliL+/BXOZ0SdmXQN9yx344+G278wwc40dIThZYREY0+w4S/3SkhInzs63sasONYC3615fCItomIKFYME/4Op4RQ0j/ckr/J5Hqgk9NFiShOxH34q4O7DqeEUPr+zjDTX10SiNNFiShexH34qxxysOcfLpPyQJ4oRkTxIu7DX41rbXCHXfZRwp8dfyKKF3Ef/iq7Q0KoIR7mWp9mE3v+RBRfDBP+2jp/pD1/1vyJKF4YJvy1Uz3DH/DlbB8iii+GCX+nZqqnPcwQV2f7hPuhQUQ0VsV9+Ev3ev6D4R9uD164yz4j2TIiotiJ+/BXaef5hztw6y4XsexDRHHCWOEvBm+HQ/32wNk+RBQvjBP+cnDAN9xZO2qtn7N9iChexH34q3P6HU4n1K5/uAO+6uG8DgARxYu4D3+Vwxl57d6pWR+IiCgexH34D9brncOo+atln5FsGRFR7MR9+Ku0J3mFG/7q4ZztQ0TxwjDh7zrJK7JlGtQPC5Z9iCheGCb8tYO84ff8pcffRER6F/fhr13SOdKyj5r5DH8iihdxH/6q4Zzkxdk+RBRvDBP+nX12d3iH24N3D/j6edjPXt2P2f/12nCaR0Q0qiyxbsBo+e3bR9y3wz/JK/CHRsW2usgbRkQUA/Hf8/eR1+FO2ZQs+xBRnIn78L9+8ZQh29QQl1KisbMv6HNwnj8RxZu4D/+HrlswZJta9nly+zEsfXALjjR2BXwOLuxGRPEm7sPfFzXM3ztyDgBQ39wd5HjX3w4n8OqeM3j8ndqoto+IKNoMM+Cr5XC6/h68pm/g47XlnjufqQEArP10cRRaRkQ0OgzZ83c4XemvzvsPtlTz4Cwfln2IKD4YNPxdf6tr/QSL9MH1/KPXJiKi0WSI8LeYhMf9Zz6sR3V9i7vsE6znr+4Plv282AsR6YUhwj/R4vkyGzv7cf1jH7jLPsGoZZ/g5aGImkdENOoMEf4pCWaf24XS9w/WYXeXfYL8HJ4ERkR6YYjw/+utS31uV3v+/jK7tqkLRfduxs6jLQBC+ZBg+BORPkQ1/IUQZcqfdX72r1b23x3NdswrGIclRdlDtpvcA76+Q7v6WCsAYMuBRgDBw53ZT0R6EbXwF0KUAbhBSlkJoEQIUeK1vwQAlP1t3vtHmtnko8AfpOefn5Hkcb93wBHwZ/j6cJh1/2v49ZbDIbWRiGi0RC38pZSVUsq1yl2rlLLG65A1ANqU23UAyqLVFgCwmIa+1GCzfRLMno/xXg30V5WHcdVv3nPf97X8w4DDiUffPBRma4mIoivqZ/gqJZ21PnZlAWjR3M+NZjtMPnr+Ish0H7tyMpg/v6j0DHUZ+HAiojEj6gO+Usr1ANYKIbLCfawQolwIUSWEqGpqahpWO7zn+gPanr/vx0R6xS9/OvpseH3vmbCek4goGqJZ89fW+esAlHsd0gYgR7mdBaDZ+zmklBVSylIpZWl+fv6w2uOr5q9u8jfgG+lFX9z3vR7/vRd24/a/1aCuqQtSSrR2D4T1/EREIyWaPf8yeIZ7HQBovgFsAGBVblsBVEaxLb57/krZx+GnXBNuz9+75u/9YXCipQcA0Gtz4Kn3j2HRT9/EsXOBVxQlIoqGaIZ/BQCrEKIcAKSUm5TtW5T7NYB7VlCbjwHhEeWr569ucfip7Yfb8/eu+vhb/19AYMt+1/TR48oHAhHRaIragK+Usg2uDwDv7Ys1t4fsj5ZAPX9/Ie/vQ8GfoWWf4Mf6nIJKRBRlhjjDFwDMvqZ6Krnrr7xjd/jv+fuaHur9NP4GgM919eP92maPNhARjSbDhL+vnr9KG/J7T7Wj6N7NqG3qCljz97XLe4DXX9nn5qd2um+bmP5EFAOGCX+z2f9sH20P/e+7TwMAKj85G7Dm7+uDwbun72/ev/axDH8iigXjhL+PkFVX9dSGvLsUJGXAeftf//OOIdu8Pw9CueA7S/5EFAvGCX8lZW8snYLvlM0CMDi/X9sTVz8knE4ZsOb/Qd2Q0xKGfFh4fzvwdUZxsLOMiYiiwTDhr9b8TUIgweLZ49eGvPoh4ZThz/P3HgQOdp+IKFYME/5qzd9kEkhUFmzrt7uK8topnWpP3Cll2PP8vU8WC6Xsww8EIooFw4T/YM9/8LKO/TZXWtv9lH2GO88/lG8OvPoXEcWCYcJfnedvEsK9VHOfzbU+v0fNX/kXcUTQ8x8y28fr4b7q+8x+IooFw4S/Wl4Zl5KA5ATXy+4ZsAPwDP/Bsk8kNX/P+6E8npd+JKJYMEz4N3X2AwAmZCYjJcG1qkV3v6vn71H2MWlm+4Rd8w+8sJsvgY7Ze6od1fWtYbWBiCgUUb+Yy1hxpr0PADBpXDIsSm2nq39oz9/sXunT30LP/g1Z2yeE8A/07eALylXCjj18ZZgtISIKzDDh39bjWjt/4rhkd4+/Wyn72J0SA3YnXtt7xn2Sl1MCzgDz/H0ZcpJXCOPF/j4gXvn4dFg/m4goHIYJ/0duWIhn/lWP2RMz8cnpDgBAW48NgGuq52Nba/GLykNYXJgNILJavPe0zZDKPn4+IP7j2V1h/3wiolAZJvxnTcjAT66eDwBISTR77LM7JZq6XGWhBqU85Cr7DK/mH9JUTw74ElEMGCb8tbzD3+GU7umfvcr0T2eQtX18CXVJZy2e5EVEsWCY2T5aKQme4d9nc7jP+u0dGAx/u0NiYmZyyM87tOwT/DGhjAsQEY00hj+A1m7bkJ6/wyld3wgsoS+8NnTAl2UfIhqbDFn2UU/yUp3r7ncv+aCyOVzhb/FxBTB/vIPc+5uAr48Rln2IKBYM2fP3XmahuWtgyLr6L+86hb/vPg2zSeC1u1bgu5fNCvq8wdb28RXzXNuHiGLBkD1/b+29NnQrtX5vZiEwZ1ImWroHgj6Pdy8+lJIOs5+IYsGQPX9fHtta63O7WXMdgGC85+wPWdjN52OY/kQ0+gwb/q9861MhHWdRrgNgDuF6ixEt6TzMJSCIiCJh2PDPTU8M6Tg19M0h/EsNCf9hLuymsnE+KBGNMMOGf5JlcLrnptuX+z1Onf8fStnn9b0NHmMDoczkCaXsw/AnopFm2PDXTu3MS0/ye1xGsmtMPJSyz/9+dBrfe+Ej9/1QMvuTM51Bjwl0IXkiokgYN/w1dZy0JP+TntR9ofT8AdfMIVWgks6TNy8BADy34zg+PtmG3Sfa/B5rC/NykkREwRg2/BPMg2GeHiD81TN/LebQwr/meBsONLhWDdWWdN4+2IjjLT3u++rqoQDw6y1HcPXvtuONfQ0+n5M9fyIaaYYNf+2JXt5n/GqpHxIZyQkhP/flv3wXgOeA781P7nRfPAYYvGgMAHT1u74tbNl/1ufzseZPRCPNsOGv5evC6iq13JOdGnr4A67B3kBjudoykvrN460DjfjmMzXo6LN5HGtjz5+IRhjP8A3Cogz0ei8GF0xXvz3gTB7tkkF9NlfP/lzXADbvOYMLpmZ5HGtnzZ+IRhh7/kGo1/sN9O3Al6bO/oAnZ2l7/v12z6UlvM8PsNnZ8yeikWXo8F8xMw93rZoJAPjDV0pw1cKCIceEMsXTl6bO/oCzfcwe4e/Zs/d+nDrb52BDJ96vPRdRe4iItAxd9nn61gvdty+fPwk9Aw78Y7fnhdMzAswECqSpK3D4mzQfKt2agWAAWP/6QY/76myfz/1yGwDg2MNXRtQmIiKVoXv+3lK9Lu+YaDbhtpXWiJ7LVfYJ7djOPnvA/XbO9iGiEcbw15jgdcnGR9csRLJmoPet730a935+dkjPdehsZ8jXAA4W/gMMfyIaYYYu+3ibVzDOffvFOy5CyTTPWTfW/HQsmNwX9HkKc1Ox42gLNladDOnnqpeO9Md7TICIaLgi7vkLITJDOKZc+bPOz/516nGRtmMkJVpMWDEzD19fXojFhdk+Z/j4m8FTrikPLS3KQW1TN+wjtBRzX5APByKicIUc/kKInUKITCHEdCFEC4CNQojHAhxfBqBSSlkBwKrc91YuhKgFUBd2y6Pk6VsvxE+unu93f06aaynowtxUj+33XD5YDpqUleKxL8IJQ249fq4yRkQUqXB6/llSyg4A5QAeklJ+DoCvQFdZNfvrlPvebpNSFkspK8NoR0zNnzwOz5cvww8+P8e9LTs1wWNK6ESvsQPvi8P7kpbo/yQyhj8RjbRwwv+oEOI6AGsBbFS2tfs7WEpZofT6AaAEQJWPw6xCiDIhxN1htCPmlllzkaSsB7RyVj52/eizHvvzMzyXiB6XEnxpiPRk/8MvLPsQ0UgLJ/zXAlgK4AYp5TEhxPkAHg/2ICFECYAaKWWN9z4p5Xql15/rqyykjBdUCSGqmpqawmhq9KnLPmiXcLhu0WQAQJbXOkBP3bw06PMFWla6ZyDwbCAionCFE/4vAPgZgDohRDOAR+Dq0QdTJqW8x3ujEuyrlbvN8FEWUr49lEopS/Pz88NoavSpZ+hq19155IaFOPDTy4f09GdNyMB1JZMDPl+gk8l+93btkBPBiIiGI9Ka/8Mh1PwhhCiXUq5Xbpcpf6vzJ6sAqLX+YvguC41ZJnfPf3Cb2SSQnGBGpmb55w3ly2A2CSRoVnL70RfmDnm+lAA1fwD4/dYjw2wxEdGgqNX8lbBfJ4SoFUK0anZtAQClDHSj0vuv9VUWGsvUAV5fF2nX9vxTE109+gSL6/ivLS/ELZ+aPuQxiZbA4c8LuhDRSArnJK+1yh+15j8dAWr+Si0/28f2xZrbFd779cKalwYA+PLSaUP2aS8Oo/boLUrPP9PPRWG0l5X0uT+EGUNERKEKOVGklEcBPARguhDiG65N8omotWyMy01PwrGHr8TqxVOG7NOeHKauF6Qu2+xv5o/2vIGFXuv5A4D2C4YMcdkIIiJ/wjnJaxVcJZssADMAvCmEuDZaDYsXavirc/X9hf+U7METw2ZPyBiyv613wH17hE4cJiIDC6eW8LAy8+bnUsp7ASwBcF+U2qV7k8a5TvRSyz7d/a7wz0wZrLR9eN8qvHjHRZgzKdPjWgLqOQQLp2bhVmV8oK1n8NKOgS4SQ0QUinBq/h6LFEgp24QQQ2r65PLSnRdhx9EWJCkDub0211RN7Xz+CZnJmJCZjNfuWuHx2CSlvn/lgokoX1mMTdUnPcI/1NVCiYj8CSf8K4QQb2BwkHctQjjJy6gmjUvB1RcMzu1Xe/7e1wzwRf3AGFBW80xJMKO5e7Dsw54/EQ1XOAO+FQB+ANdZvp8FsF5K+Ui0GhZv1LN01amfgVw623VC2zJrrvIYM1q6+937fU0vJSIKR1jr+Stz8XU1H3+sWHf9+Xj0zUMozk8PcMwCpCRasLgwB3U/u8J9IllKohnHW3rcxznZ8yeiYQoY/kKIIwD8JY2Aa7rnzBFvVRxaNC3b45rBvqxZMnjOgPYav6mJZo9rA7DsQ0TDFaznvzjIfhoFKV6lIpZ9iGi4Aoa/lNLv8g00eqbleF4cxsmrOhLRMHHNAB24dpHniqDs+RPRcDH8dWB6nucgMQd8iWi4GP46kO11cRgO+BLRcDH8dUC7UBzAsg8RDR/DX4dY9iGi4WL46xB7/kQ0XAx/naj+YRl+fJXr8o8HGzqx7dDYuqA9EelLWMs7UOzkpidhSrbrgi93Pf8RAODYw1fGsklEpGPs+etIkCs9EhGFjHGiIyavWT9ERJFi+OuI2cTwJ6KRwfDXETN7/kQ0Qhj+OmJiz5+IRgjDX0csDH8iGiEMfx1hz5+IRgrDX0dY8yeikcLw1xHv2T6SyzwQUYQY/jriPc+f67sRUaQY/jri3fPnuv5EFCmGv44M2D0v3svwJ6JIMfx1ZFJWssd9Lu1MRJFi+OtIXnoSHll9vvs+e/5EFCmGv86kJQ2uws3wJ6JIMfx1JktzMXeGPxFFiuGvM8utubDmpwEAnKz5E1GEGP46I4TA2pVWAICdPX8iihDDX4fUk72cDH8iilBUr+ErhChXbhZLKe/xsX81gDYAJVLK9dFsSzyxmF3hz5o/EUUqaj1/IUQZgEopZQUAq3Jfu78EAKSUlQDa1PsUnNrzZ9mHiCIVzbKPFYAa+HXKfa01cPX61f1loJCoyzxwwJeIIhW1so/S41eVANjgdUgWgBbN/dxotSXeqBd1YdmHiCIV9QFfpZxTI6WsieCx5UKIKiFEVVNTUxRap09q2YfhT0SRGo3ZPmW+BnvhKvnkKLezADR7HyClrJBSlkopS/Pz86PZRl0xs+dPRMMU1fAXQpSrs3jUAV8hRJayewMGxwGsACqj2ZZ44g5/1vyJKELRnu2zTghRK4Ro1ezaAgBqGUg5ri2SspBRsedPRMMVzQHfSgDZPrYv1tyu8N5PwZlZ8yeiYeIZvjrknurJ8CeiCDH8dUgNf57kRUSRYvjrkIkDvkQ0TAx/HTJzYTciGiaGvw6x7ENEw8Xw1yEO+BLRcDH8dYgneRHRcDH8dYgneRHRcDH8dYgneRHRcDH8dYg9fyIaLoa/DqmXcbQ5GP5EFBmGvw4lW8wAgH67I8YtISK9YvjrUHKCK/z7bM4Yt4SI9Irhr0NJFtevrc/Gnj8RRYbhr0Mmk0CixYQ+ln2IKEIMf51KtpjQz7IPEUWI4a9TyQlm9A6w509EkWH461RygpllHyKKGMNfp5ITTBzwJaKIMfx1KiXBzKmeRBQxhr9OJSWY2fMnoogx/HXKVfNnz5+IIsPw1ynXVE/Pnr/d4UTRvZvx27cOx6hVRKQXDH+dSvZR9jnXNQAA+PWWI7FoEtGYIKXkirchYPjrVFqSBa09No9LOZ7t6HPdEDFqFNEYULGtDsX3vYr2XlusmzKmMfx1apk1B+29NvzmrSOo2FYLKSUalPBn9pORPb/zBACgqbM/xi0Z2yyxbgBFZtWcCZiQmYRfVB4CALT22FAwLhkAIJj+RBQEe/46lZ5kwfrVC933H9tai//6v30AALtD4ui57lg1jYh0gOGvYytn5uGXay7AjvtW4ZLz8gEAk7NSkJmSgBv+8MHgGACRgahffKXkoG8gDH8dE0LgmkWTMT4zGX/6+hL89Op52Hj7cjxfvgzd/XZ8+7ldsDt4LgAZEy9zGhjDP06YTQJfXV6EgqwUzJqQgQeumY8Pj7bgodcOxLppRDFhd7LjEwjDP05dv3gKvr68EH967yj+ua8h1s0hGnU2fusNiOEfx+6/ci7mTMrEfS/v5bQ3MhyWfQJj+MexRIsJv1izEB29Nvx+K8/6JWNhzz8whn+cmz0xEyWFWXhy+zFU17fEujlEo8bOnn9AUQ9/IURJgH3rlL/Lo90OIyvMSQMAXP/YB1wGmgxjgD3/gKIa/kKIMgAbAxxSLoSoBVAXzXYY3Z2XFsOa5/oAqG/uiXFriEYHe/6BRTX8pZSVCBzst0kpi5XjKEoKc9PwPze6zgY+0cLwJ2NgzT+wWNf8rUKIMiHE3TFuR9ybmpMKADjRyvAnY2D4BxbT8JdSrld6/blKiYiiJDctEamJZvzkH59gw87jONfVj90n2mLdLKKo4VTPwGK2qqcyyNsipdwEoBmANVZtMQIhBAqyUnCksQv3vLgHwB4AwMpZ+XjwmvnIz0hCcoI5to0kGkFj7QzfPpsD/7vrFNYsmQoxBpbeHfWevxAiS7lZBUCt9Rcr972PLRdCVAkhqpqamkariXHL7OMNt+1QE1asfxulD1Ri94k22BxO99flM+29UV8cy+5wcgEuioqBMXaN619UHsK9L+3BPz85G+umAIhyz18IsRpAqRBitdLDB4AtABZLKWuUcG8BUCulrPF+vJSyAkAFAJSWljIhhumOS4rxnxs+AuDq8T/xtcWobezGH9+rw6t7zuDq320HABTmpuIrFxbiwVf3IzPZginZqTCZgM/OnYjJWSlYODULM8anB/15fTaH+9vE3lPt2HuqHTMnZGBeQSZ6Bxxo7u7H9zd+DJMA7rl8Nho6+mASAlctLPD5fKfaenHnMzW4btFkfHnpNCRahvZdpJQevar2XhtSE81IMMd6eItGmz2CSzlKKfHGvrNYNWf8iL9nGjtcZ9l39tlH9HkjFdXwVwJ/k9e2xZrbFdH8+eTpmkWTcc2iyTjQ0IHctCQkWcyYW5CJR2+8AF9ZVogHN+9Hc1c/jrf04MFX9wMA5k8eh+4BB8629+HRNw+5n2tCZhJy05LwjRXTkZpoxqm2PuSkJWDB5Cw8++Fx/Hn7Ufdx5SuL8bNX97uvq5poNg2Zg72m4l/u2wsmj8OHR5sxd9I43LVhF8609eHy+ROx+0Qb6s51Y/eJNry86xT+s2wmLjlvPADgYEMnHtj8CeqaunHzxUWYNC4Fv996BPtOd+C6ksl49MYLYHM43f+hHU4Jsyn4V+9+uwNPbT+Gqy+YjInKxXLGuiONnXip5hS+/9nzYArhNcYrWwQ9/22Hz+H2v1Xj26tm4ruXzYpCq8YOXsnLgGZPzByyrWRaNl684yIAQH1zN9a/fhC3fGo6Fhdmu4+pOtaC/Q2d+Oe+BrR0D+B4Sw+++8LugD/rbEc/fvrKJ0hNNOO2FVaYhMBLu06ipWsAa5ZMhd0pccWCSajYVofK/a6vw5f8fOuQ53l51ykAQHZqAi4qzsPWg4246cmd+OqyQrT0DODDumb3Bewf2Lzf47Ev1ZzCtJxU/LLyMC6ekYui3DS8WHMSJdOycehsF6Zkp2DV7PFo7OzH5j1nMCU7BSYh0NDe57405kOvHcCF03Mwr2AchADmFWTiiwsLcKSpC0W5afjkTAfe2t+I71w2y+eHypn2XuSkJSLJ4jmuMmB3onL/WVw+b+KIBfUdf6vB4cYurF48Bdb84N/Q4o3a37dF0PM/p6yBVd8cvYshjZVZSAx/GqIwNw2/+/ehJ2aXFuWgtCgHX11WCADoGbDjlY/PIC3Rgt0n25CTlog/vluHb6yw4stLp6Hf5kB6sgUnW3uRnZqI/IwkAMC3PjMDDik9vlYvnZ4DAPjBS3vw3I7jAIDl1lw8fP0CTM1OxXtHzuFrf96B+6+ci9WLp6B3wIHvb9yNp/9V736Omy8uQvlKK8xCYO/pdkgJvHWgERt2nsAvKw8DALYfacb2I80AgPdrmzFzfDo6em34H+VbTUayBa09AxiXkgCL2RXGQgB56Un48GgLPjw6uESGrw++1/c1YNWc8ZiWk4p5BePw/I7jONLYhar6VkzMTMbTty7F1JxUvFB1Amc7+vC7t2sBAF9fXoiVs/Lx8q5T+OalM2B3SKQkmpGXnoj0JAvsTonGjn6kJplhMQn0252YkOn6JtLdb8eT24+iqr4VV51fgMONXe623HnJDPcxx1t6MC0nFWlJo/ff/oltdZg5Id39DW00qOEaKGQH7E4IgSGlHUcEHxjh6uwbGxeWF3oZbCstLZVVVUPGhGmMcTrlsHqwH59swxd/ux1/uWUpPj0r32NfU2c/8tIT3TV9p1NiY/UJPPvhcZQW5eDbn5mJcakJQ57zoxNt2HqwEX969yiWF+di9eIpmDMpE398tw7fuWwWslIT8cnpDvTZHSiZlj3k8U6nRHuvDZf8fCvae22YNSEdy625+MsHgx88F8/IRV56Et7Y14A+29DQueXi6fj77tMAJCZnpWD3yfaw/l1MAtDmktkksLpkCkwmgdrGLuw45n/dpqsvKMA/dp+GUwJJFhPSkyy445JiLC/OxbyCcTh0thN3PlODLy+dhgun52D+5HHoHXCgZ8CO7n4HctMT8e3ndqHX5kB7rw03lk7FVQsL8FLNSSycmoUBuxPzCjJRXd+KC6Zm4Z1DTfjuC7vxfPkyfEkp5+3+8WdhNgmkJ1nQ3NWP//voNL60dCpSEz0/iA6f7cSM8enDmg1z4c8qcbajH7etmI77r5zr85jSByoxPiMJr961wmP7Y1trse71A7j6ggL86kuLIm6DL996bhf+sfs0/uPSGfj+584b0ef2RwhRLaUs9bmP4U9jTb/dMaQ8MhZ09dvR3W/HhMxkSClxsrUX/7vrFJITzLhtpWumspQSDR19+MzP30Gvso7Si3dchMWF2dj88Rl881nXvIZrLihA94ADD1wzH5X7z+L+l/e6f861iyZjmTVHmZILlEzLQs1x1zkZEzOTYRLAIuX5VMutuXjy5iU4fLYLbx9sxL7T7Xhj3+CskuQEE+69fDaq6lux7VATOpRBR+8PFQBIMJq9tBcAAAmnSURBVAtIOThgajYJOJwS4zOS0BhkafCJmcnuUpkvy6w5aOzsR11TN86bkAFrfhp2HmvBtYsmY+WsfHz1Tztw7+dnY8XMPLT32LBkeg4SzCbsPdWOf+5rwCWzx2PPyXZMy03F9sPnMOBw4uIZebj0vPF4seYkxmck4e5NH6O5ewBzJ2XiqVuWYMDuxLuHz+Gi4lwkWcywOZxYsf5tAMCxh6/0aN8Dr3yCP77nGq/62vJC/PfV8wO+Xl/ae2x490gTvnC+58SFm57cga0HmyJ+3kgw/IlG2YmWHuw73Y73a5vxoy/MhUUpLzy34zim5aTi4hl5HsefbO1BzfE2mATcoVFd34q2ngGsmjMBlzzyNo419+DoQ1fA7pTuQMxNT8Sp1l7MLcgc0ouurm/FkcZO3PPiHiwpysbG211jOgN2J3afbMOT24/i1T2uC/18fv5EHGjoxNFzQ2vdSRYTlk7PwV9vWYrGzn4cPtuF946cw5xJGTjV1ov1rx8E4CqNnesa+uHwuXkTUHWsFRJAS7drXOYry6bhg9pm1DYFr60/8bVSfH/jbrT3+i+XpCaa0TMQ2aKF910xG7d+yoq3DzTix3/fh1Ntve59B356OR7bWou8jCQst+Ziz6k2LJichTf2NSAvPRGTxqWgMDcVj22tRWFuGm7/tBVf/O127DnVjnXXL8CaJdPQ1W9Hc1c/yv9ajYNnOwG4rrW99tNWCCGQkWTBxTPy0Gdz4L6X9yA/PQk3XzwdMyekD/vcG4Y/kc519tnQM+Bw1/nD8a+6ZkzNScXkrBSP7X02B97Y1wCLyYTL5k5AosWErn47PjrehgWTx+GZHfX4t6XTkJZkgcUk/JZiXt1zBk+9fwxP3rQEWw824Ux7L3Ydb8PmPa5vJq/dtQKzJ2bgzmdq8NreBvzwyjn4xgornE6J2qYu1Df3oL3Xhu9t9D95wGwS2Hj7cuw73YGmjj78+q0juG7RZKxbfT5+/s+DePwd1xJiN11UhCnZKbiuZApuenIHPj7Zjotn5MLhlJiel4765m68X9vsft6UBDP67A7MGp/hDmatGePTcUQZQ/FH+23nh1fO8Zhw8MA18/GHd2pxsrXX38MD+vcLp+HBaxdE9FiA4U9EMfD8juN460AjfvNvi5BkMeMfu0/jW8/twpvfWYmZEzKGHF/b1AWHU2JGfrp7QkDRvZsBwN2LVlXXt+C8iZlIT7Jg1/FWXPv79wF4lnHOdvThxZqTKF9hdX/zUkkp0djZD4tJ4JvP1uB0Wx+OK4se7v/vy/HKx6fx/zZ9HPJrDfYhsXrxFNx0URHmTsrEr7YcRmZKAmaMT8e3nq1xl+BUb35nJf747lFsqDoBAKj6YRny0pNCbosWw5+IxoSufjvSw5ht9PKuk9h/phP3XTHH7zG9Aw7M+dHrmFeQic3fXuH3uGCe2FaH9l6bezD2mQ/rcf/Le/G5eRNwyXnjse1QE25bacV1v38f4zOSsOP+MjicEn95/xgunT0eh852Yu3T1R7P+aUlU/GF8wvwqZl5vn4kznX1o/QBz0WN1Q+wD2qb8eUn/oX7rpiN8pXFEb0mhj8RxbWdx1pQnJ+OnLTEqP6cfrsDN/zhA3z3slk+p6/WN3fjxeqT+PaqmTAJEdLMt94BB378973YfaIdF8/Iw4+ucs1QklLi5qd2Ys/Jdmy/9zMR1f8Z/kREOnS6rRd2h8S03NSIHh8o/HmSFxHRGFXgNUg/krjaFRGRATH8iYgMiOFPRGRADH8iIgNi+BMRGRDDn4jIgBj+REQGpJuTvIQQTQDqgx7oWx6AcyPYnFiKl9cSL68D4GsZq/hagEIpZb6vHboJ/+EQQlT5O8tNb+LltcTL6wD4WsYqvpbAWPYhIjIghj8RkQEZJfwrYt2AERQvryVeXgfA1zJW8bUEYIiaPxEReTJKz193hBAlXvdXCyHKhBB3B9pGFI5g7ye+x+JXXIe/Xt+4QogyABs190sAQEpZCaBNCFHia1tMGhuEEKJc+bNOs02XIaO0ryweXgvgfp9dptzW5XtM/V0IIco123T5O1H+zVcLIVZrtkXttcRt+OvhjeuP0uY6zaY1ANqU23UAyvxsG1OUcKmUUlYAsCpvWL2GTBmAG5Q2lvhrtx5eix+6fI8BKBdC1EL5/6Lz38kPpJSb4Pq/EvX3V9yGP/Txxg1VFoAWzf1cP9vGGisG/93rlPu6DBkpZaWUcq1y1yqlrIFOXwvgCkklQFR6fY/dJqUs1rwWXf5OlN7+TgCQUq4fjfdXPIe/Ht64cU1KWaH0+gGgBEAV9BsyANw1cvVDQM+vJSfWDRghVq8SiF5/J0sA5Cq9+1F5LfEc/vGkDYP/WbMANPvZNiYpX01rlN6Mrkkp1wNYK4TIinVbIuWj1w/o9D2m9JIr4QrOMdejD1Oz+n9EW/ePlni+hu+Yf+OGYQMA9dRuKwD1P66vbWNRmZTyHuW2v9/LmP5daWqtNXB93S6HTl8LXL1lK1ztzFFem+7eY8ogb4tSJ2+Gq416/Z00Y3Ccrw2ubwJRfS3x3PPfANebARiDb9xAlE/9UvXTX9MbKAPQJqWs8bUtZg0OQAhRrvSW1bb6+r3o4XdVBs//dHXQ6WuRUm5SAhNwvRa9vseqMPjvW6zc1+XvBMAmDLYxC676f1RfS1yf5KX0DOrgGqCLp7P9dEEzZbUFruC8QUpZ6ev3MtZ/V0qZ50bl7mJ18FePryWeqL1/uP6t12u26e53onktS9RvytF8LXEd/kRE5Fs8l32IiMgPhj8RkQEx/ImIDIjhTxRlQogsZQkCojGD4U9EZEAMfyIiA2L4EymEEHcLIWqFEG8qpZos5f7j6t+aY8uFENXKnzKv5ylXjq/1Wv75ce9tRLHCef5EcC/f8AMp5Q3K0gf3KH+OApgupWwTQlQr21oArJNSqmvhVwNYpRxTAuAJKeViZZ9VOb4VQLGUsk4I0SqlzB71F0mkEc9r+xCFYw1ca9686bW9RUqpLqH7OJSLnyi3VRvgOvu3QnmeDeoOJeyz4FoeQV27pUUIkaV5XqJRx/AnGvSQZs0bdUkHb81wrSOjlYvBNdb9aQmyn2hUseZP5LIBg+v0a6+hbNV8CNwD10JaG9VjlX2rMbjA1ga4ev/ez0M0pjD8ieBe1XKjOoiLwaWM2wA8oczT36SsdlkJ4E1l2xYA96glHeV51IHdWozBq0YRARzwJfJL6dVXSym9yzxEuseePxGRATH8iYgMiGUfIiIDYs+fiMiAGP5ERAbE8CciMiCGPxGRATH8iYgMiOFPRGRA/x++W6lxD9a+7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(train_loss)\n",
    "plt.xlabel(r'epoch',fontdict={'size':12},y=0.1)\n",
    "plt.ylabel(r'loss',fontdict={'size':12})\n",
    "plt.savefig(\"loss.pdf\",bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3385619430048774\n",
      "1.8486078153821617\n",
      "1.3312275608205664\n",
      "4.67299425586396\n",
      "3.3245961721403825\n",
      "1.6384946555278073\n",
      "1.3836944325163298\n"
     ]
    }
   ],
   "source": [
    "def Fi2fi(Fi,q):\n",
    "    fs = []\n",
    "    for item in Fi:\n",
    "        f_i=np.random.binomial(item[0],q,size=item[1])\n",
    "        fs.append(f_i)\n",
    "    fs = np.concatenate(fs)\n",
    "    fs = fs[fs>0]\n",
    "    unique, counts = np.unique(fs,return_counts=True)\n",
    "    f_sparse = np.array([unique,counts]).T\n",
    "    return f_sparse\n",
    "renEsti=ndvEstimator(\"learned_ndv_estimator/model_training/model_paras.npy\")\n",
    "import os\n",
    "Dir=\"kasandr\"\n",
    "q=0.01\n",
    "kasResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(5):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        r4.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r5.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    kasResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "kasResult=np.array(kasResult)\n",
    "for i in range(kasResult.shape[1]):\n",
    "    print(np.mean(kasResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.4986996050782586\n",
      "2.215426108400742\n",
      "7.746348827979548\n",
      "4.364261793293312\n",
      "1.67873266605254\n",
      "1.6944794950701172\n",
      "1.2933506012672242\n"
     ]
    }
   ],
   "source": [
    "def Fi2fi(Fi,q):\n",
    "    fs = []\n",
    "    for item in Fi:\n",
    "        f_i=np.random.binomial(item[0],q,size=item[1])\n",
    "        fs.append(f_i)\n",
    "    fs = np.concatenate(fs)\n",
    "    fs = fs[fs>0]\n",
    "    unique, counts = np.unique(fs,return_counts=True)\n",
    "    f_sparse = np.array([unique,counts]).T\n",
    "    return f_sparse\n",
    "renEsti=ndvEstimator(\"learned_ndv_estimator/model_training/model_paras.npy\")\n",
    "import os\n",
    "Dir=\"kasandr\"\n",
    "q=0.005\n",
    "kasResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(5):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    kasResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "        \n",
    "\n",
    "kasResult=np.array(kasResult)\n",
    "for i in range(kasResult.shape[1]):\n",
    "    print(np.mean(kasResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.4464644170645795\n",
      "3.7264714943998016\n",
      "30.492334252481466\n",
      "7.56339028107088\n",
      "3.864506283795677\n",
      "2.2479982056866126\n",
      "2.15549028228203\n"
     ]
    }
   ],
   "source": [
    "def Fi2fi(Fi,q):\n",
    "    fs = []\n",
    "    for item in Fi:\n",
    "        f_i=np.random.binomial(item[0],q,size=item[1])\n",
    "        fs.append(f_i)\n",
    "    fs = np.concatenate(fs)\n",
    "    fs = fs[fs>0]\n",
    "    unique, counts = np.unique(fs,return_counts=True)\n",
    "    f_sparse = np.array([unique,counts]).T\n",
    "    return f_sparse\n",
    "renEsti=ndvEstimator(\"learned_ndv_estimator/model_training/model_paras.npy\")\n",
    "import os\n",
    "Dir=\"kasandr\"\n",
    "q=0.001\n",
    "kasResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(5):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    kasResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "        \n",
    "\n",
    "kasResult=np.array(kasResult)\n",
    "for i in range(kasResult.shape[1]):\n",
    "    print(np.mean(kasResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.2358821838854206\n",
      "1.1672712291220557\n",
      "1.2565958701992543\n",
      "1.0772652796877626\n",
      "1.2565672860583708\n",
      "1.2678685416723021\n",
      "1.203259003162384\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"Airlines\"\n",
    "q=0.01\n",
    "AirResult=[]\n",
    "order_D = []\n",
    "\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(5):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        \n",
    "        fi=np.array(fi)\n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    AirResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "    \n",
    "AirResult=np.array(AirResult)\n",
    "for i in range(AirResult.shape[1]):\n",
    "    print(np.mean(AirResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3721122722298031\n",
      "1.2533745410144042\n",
      "1.3343687892461717\n",
      "1.1403129986690168\n",
      "1.1831470820596173\n",
      "1.2553913576140947\n",
      "1.271225130558014\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"Airlines\"\n",
    "q=0.005\n",
    "AirResult=[]\n",
    "order_D = []\n",
    "\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(5):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        \n",
    "        fi=np.array(fi)\n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    AirResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5),\n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "    \n",
    "AirResult=np.array(AirResult)\n",
    "for i in range(AirResult.shape[1]):\n",
    "    print(np.mean(AirResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7832722908252263\n",
      "1.4473927456413942\n",
      "1.6150182167561862\n",
      "5.635739767169734\n",
      "1.3955923639503738\n",
      "1.5085251543774922\n",
      "1.4197829604148864\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"Airlines\"\n",
    "q=0.001\n",
    "AirResult=[]\n",
    "order_D = []\n",
    "\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(5):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        \n",
    "        fi=np.array(fi)\n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    AirResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "    \n",
    "AirResult=np.array(AirResult)\n",
    "for i in range(AirResult.shape[1]):\n",
    "    print(np.mean(AirResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5789030603429637\n",
      "1.0468932433629379\n",
      "4.361799022401808\n",
      "5.461873079809291\n",
      "1.2592537658240268\n",
      "1.304217646743238\n",
      "1.7761901259896486\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"SSB\"\n",
    "q=0.01\n",
    "SSBResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(10):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        #GEE\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    SSBResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "SSBResult=np.array(SSBResult)\n",
    "for i in range(SSBResult.shape[1]):\n",
    "    print(np.mean(SSBResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8250919168908823\n",
      "1.054498351384307\n",
      "7.863474963202286\n",
      "8.337910799643115\n",
      "1.5498386432008269\n",
      "1.467215689591704\n",
      "1.637346257218074\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"SSB\"\n",
    "q=0.005\n",
    "SSBResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(10):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        #GEE\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    SSBResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "SSBResult=np.array(SSBResult)\n",
    "for i in range(SSBResult.shape[1]):\n",
    "    print(np.mean(SSBResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.773466806492016\n",
      "1.0709174997442172\n",
      "35.944947526873676\n",
      "25.740329495700912\n",
      "3.984417438738141\n",
      "1.5444374362901205\n",
      "2.5350985329014706\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"SSB\"\n",
    "q=0.001\n",
    "SSBResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(10):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        #GEE\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    SSBResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "SSBResult=np.array(SSBResult)\n",
    "for i in range(SSBResult.shape[1]):\n",
    "    print(np.mean(SSBResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8785342448652207\n",
      "7.586876009287286\n",
      "8.969416825686945\n",
      "1.2760317176292197\n",
      "1.917025881131156\n",
      "1.7172723161236918\n",
      "1.359943182635345\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"NCVoter\"\n",
    "q=0.01\n",
    "NCVResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    if filename[0]=='.':\n",
    "        continue\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(10):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    NCVResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "    \n",
    "NCVResult=np.array(NCVResult)\n",
    "for i in range(NCVResult.shape[1]):\n",
    "    print(np.mean(NCVResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.4634740854741475\n",
      "4.348021320097957\n",
      "16.681485230631655\n",
      "2.0315868294775754\n",
      "3.037868106156685\n",
      "2.0048561608910944\n",
      "1.7746150741171953\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"NCVoter\"\n",
    "q=0.005\n",
    "NCVResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    if filename[0]=='.':\n",
    "        continue\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(10):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    NCVResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "    \n",
    "NCVResult=np.array(NCVResult)\n",
    "for i in range(NCVResult.shape[1]):\n",
    "    print(np.mean(NCVResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.36798012269765\n",
      "13.064870484809221\n",
      "68.57140469803\n",
      "13.356480174343645\n",
      "8.594464256319311\n",
      "3.6171609926823725\n",
      "2.8458315647839707\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "Dir=\"NCVoter\"\n",
    "renEsti=ndvEstimator(\"learned_ndv_estimator/src/estndv/model_paras.npy\")\n",
    "q=0.001\n",
    "NCVResult=[]\n",
    "order_D = []\n",
    "for filename in os.listdir(Dir):\n",
    "    if filename[0]=='.':\n",
    "        continue\n",
    "    pathname=os.path.join(Dir,filename)\n",
    "    Fi=np.loadtxt(pathname,delimiter=\",\",dtype=int)\n",
    "    if len(Fi.shape)==1:\n",
    "        Fi=np.array([Fi])\n",
    "    r1=[]\n",
    "    r2=[]\n",
    "    r3=[]\n",
    "    r4=[]\n",
    "    r5=[]\n",
    "    r6=[]\n",
    "    r7=[]\n",
    "    r8=[]\n",
    "    r9=[]\n",
    "    order_D.append(sum(Fi[:,1]))\n",
    "    for r in range(10):\n",
    "        fii=Fi2fi(Fi,q)\n",
    "        fi_temp=fii\n",
    "        max_i=max(max(fi_temp[:,0]),L)\n",
    "        fi_dict={ item[0]:item[1] for item in fi_temp }\n",
    "        fi=[]\n",
    "        n=0\n",
    "        for i in range(1,max_i+1):\n",
    "            if i in fi_dict.keys():\n",
    "                fi.append([i,fi_dict[i]])\n",
    "                n+=i*fi_dict[i]\n",
    "            else:\n",
    "                fi.append([i,0])\n",
    "        fi=np.array(fi)\n",
    "        \n",
    "        N=sum(Fi[:,0]*Fi[:,1])\n",
    "        D=sum(Fi[:,1])\n",
    "        r1.append( max(GEE(fi,q)/D,D/GEE(fi,q)))\n",
    "        r2.append( max(Chao(fi)/D,D/ Chao(fi) ) )\n",
    "        r3.append( max(GT(fi)/D,D/GT(fi)))\n",
    "        r4.append( max( Shlosser(fi,q)/D,D/Shlosser(fi,q)) )\n",
    "        r5.append( max(Cheby4f0(fi,1/N)/D,D/Cheby4f0(fi,1/N) ))\n",
    "        L=7\n",
    "        H=200\n",
    "        data_n=n\n",
    "        data_N=N\n",
    "        M=[]\n",
    "        for t in range(1,L+1):\n",
    "            M_temp=[]\n",
    "            for j in range(1,H+1):\n",
    "                M_temp.append( comb(n,t)*(j/(N-j))**t)\n",
    "            M.append(M_temp)\n",
    "        M=np.array(M)\n",
    "        M[M<0.001]=0\n",
    "        feat=[0.0 for i in range(70)]\n",
    "        for item in fii:\n",
    "            if item[0]>70:\n",
    "                continue\n",
    "            else:\n",
    "                feat[item[0]-1]=item[1]\n",
    "        temp=pnet(torch.tensor([feat]).float().to(device),\n",
    "                  torch.tensor([M]).float().to(device),\n",
    "                  torch.tensor([N]).float().to(device),\n",
    "                  torch.tensor([n]).float().to(device),\n",
    "                 )[0].cpu().detach().numpy()+sum(fi[:,1])\n",
    "        r8.append( max(temp/D,D/temp))\n",
    "        temp=renEsti.profile_predict( f=fi[:,1],N=N )\n",
    "        r7.append( max(temp/D,D/temp) )\n",
    "    NCVResult.append([np.median(r1), \n",
    "                      np.median(r2),\n",
    "                      np.median(r3),\n",
    "                      np.median(r4), \n",
    "                      np.median(r5), \n",
    "                      np.median(r7),\n",
    "                      np.median(r8)])\n",
    "    \n",
    "NCVResult=np.array(NCVResult)\n",
    "for i in range(NCVResult.shape[1]):\n",
    "    print(np.mean(NCVResult[:,i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py37",
   "language": "python",
   "name": "py37"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}