rooted/notebooks/hcv.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dataManipulation import *\n",
    "from utils import tree_summary, summary, summary_raw, get_support_info\n",
    "from vbpi import VBPI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the sequence data and estimate the subsplit support\n",
    "data, taxa = loadData('../data/HCV/HCV.nexus', 'nexus')\n",
    "mcmc_support_trees_dict, mcmc_support_trees_ids = summary('../data/HCV/HCV_skyride_support_short_run', 'nexus', burnin=251)\n",
    "rootsplit_supp_dict, subsplit_supp_dict = get_support_info(taxa, mcmc_support_trees_dict)\n",
    "del mcmc_support_trees_dict, mcmc_support_trees_ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up the model\n",
    "model = VBPI(taxa, rootsplit_supp_dict, subsplit_supp_dict, data, pden=np.ones(4)/4., subModel=('JC', 1.0),\n",
    "             root_height_offset=5.0, clock_rate=7.9e-4, psp=True, coalescent_type='skyride')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.float32 torch.Size([3673])\n",
      "torch.float32 torch.Size([9690, 2])\n",
      "torch.float32 torch.Size([62, 2])\n"
     ]
    }
   ],
   "source": [
    "# variational parameters\n",
    "for param in model.parameters():\n",
    "    print(param.dtype, param.size())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1000:(178.0s) Lower Bound: -8303.1516 | Logprior: -709.4425 | Logll: -7886.3276 | Root Age: 276.3767\n",
      "Iter 2000:(170.7s) Lower Bound: -8229.5442 | Logprior: -683.7360 | Logll: -8194.3271 | Root Age: 323.8568\n",
      "Iter 3000:(169.0s) Lower Bound: -8154.3307 | Logprior: -643.6879 | Logll: -8077.0947 | Root Age: 331.0407\n",
      "Iter 4000:(168.6s) Lower Bound: -8082.2354 | Logprior: -629.9238 | Logll: -7835.6763 | Root Age: 333.3037\n",
      "Iter 5000:(168.4s) Lower Bound: -8005.8367 | Logprior: -654.0965 | Logll: -8902.6484 | Root Age: 321.8905\n",
      ">>> Iter 5000:(6.5s) Test Lower Bound: -8581.5615\n",
      "Iter 6000:(167.9s) Lower Bound: -7947.3722 | Logprior: -602.7484 | Logll: -7665.4434 | Root Age: 323.5475\n",
      "Iter 7000:(167.9s) Lower Bound: -7899.1694 | Logprior: -570.8593 | Logll: -7841.6689 | Root Age: 320.6855\n",
      "Iter 8000:(168.9s) Lower Bound: -7861.2569 | Logprior: -576.4910 | Logll: -7574.8428 | Root Age: 313.6620\n",
      "Iter 9000:(171.5s) Lower Bound: -7829.8090 | Logprior: -589.8684 | Logll: -7835.7935 | Root Age: 308.0958\n",
      "Iter 10000:(168.2s) Lower Bound: -7800.1156 | Logprior: -555.5772 | Logll: -7535.9434 | Root Age: 310.6952\n",
      ">>> Iter 10000:(6.5s) Test Lower Bound: -8129.6250\n",
      "Iter 11000:(167.0s) Lower Bound: -7777.7134 | Logprior: -680.0777 | Logll: -9497.7529 | Root Age: 307.5589\n",
      "Iter 12000:(168.1s) Lower Bound: -7755.8785 | Logprior: -582.4579 | Logll: -7848.2358 | Root Age: 294.3237\n",
      "Iter 13000:(167.0s) Lower Bound: -7740.0610 | Logprior: -557.8148 | Logll: -7543.7549 | Root Age: 286.7356\n",
      "Iter 14000:(166.6s) Lower Bound: -7725.6438 | Logprior: -614.6996 | Logll: -8570.1211 | Root Age: 282.0554\n",
      "Iter 15000:(166.7s) Lower Bound: -7713.0668 | Logprior: -491.9824 | Logll: -7610.9521 | Root Age: 278.9356\n",
      ">>> Iter 15000:(6.8s) Test Lower Bound: -7904.6167\n",
      "Iter 16000:(165.7s) Lower Bound: -7701.2382 | Logprior: -514.2464 | Logll: -7457.7227 | Root Age: 268.9487\n",
      "Iter 17000:(165.0s) Lower Bound: -7693.0884 | Logprior: -504.9359 | Logll: -7479.0679 | Root Age: 265.4723\n",
      "Iter 18000:(164.2s) Lower Bound: -7686.4725 | Logprior: -497.5490 | Logll: -7420.6851 | Root Age: 259.4631\n",
      "Iter 19000:(164.4s) Lower Bound: -7679.3082 | Logprior: -518.8069 | Logll: -7441.8794 | Root Age: 255.8956\n",
      "Iter 20000:(163.8s) Lower Bound: -7674.8442 | Logprior: -516.1630 | Logll: -7415.2124 | Root Age: 251.2778\n",
      ">>> Iter 20000:(6.6s) Test Lower Bound: -7754.7524\n",
      "Iter 21000:(164.0s) Lower Bound: -7670.0725 | Logprior: -508.4063 | Logll: -7430.9658 | Root Age: 246.7732\n",
      "Iter 22000:(164.0s) Lower Bound: -7666.6450 | Logprior: -499.8655 | Logll: -7432.2285 | Root Age: 244.8114\n",
      "Iter 23000:(163.4s) Lower Bound: -7663.8109 | Logprior: -503.1346 | Logll: -7393.5332 | Root Age: 240.2775\n",
      "Iter 24000:(163.3s) Lower Bound: -7661.0921 | Logprior: -502.2778 | Logll: -7391.4268 | Root Age: 239.8238\n",
      "Iter 25000:(163.5s) Lower Bound: -7658.5769 | Logprior: -496.1500 | Logll: -7399.5488 | Root Age: 238.1252\n",
      ">>> Iter 25000:(6.8s) Test Lower Bound: -7705.6230\n",
      "Iter 26000:(162.9s) Lower Bound: -7656.6082 | Logprior: -519.0089 | Logll: -7420.7949 | Root Age: 237.2367\n",
      "Iter 27000:(163.0s) Lower Bound: -7655.1070 | Logprior: -501.6746 | Logll: -7425.8335 | Root Age: 235.3040\n",
      "Iter 28000:(163.4s) Lower Bound: -7652.7400 | Logprior: -524.6695 | Logll: -7407.0469 | Root Age: 234.5384\n",
      "Iter 29000:(163.1s) Lower Bound: -7651.3074 | Logprior: -497.9752 | Logll: -7404.0723 | Root Age: 233.7842\n",
      "Iter 30000:(163.0s) Lower Bound: -7649.9600 | Logprior: -499.9448 | Logll: -7384.2510 | Root Age: 233.6265\n",
      ">>> Iter 30000:(6.5s) Test Lower Bound: -7687.8857\n",
      "Iter 31000:(162.9s) Lower Bound: -7648.2748 | Logprior: -491.1624 | Logll: -7384.3462 | Root Age: 232.5979\n",
      "Iter 32000:(162.6s) Lower Bound: -7647.0258 | Logprior: -502.4512 | Logll: -7387.7466 | Root Age: 233.2384\n",
      "Iter 33000:(162.7s) Lower Bound: -7645.6354 | Logprior: -494.4209 | Logll: -7384.2490 | Root Age: 231.9039\n",
      "Iter 34000:(162.6s) Lower Bound: -7644.5336 | Logprior: -504.7370 | Logll: -7380.4355 | Root Age: 232.1502\n",
      "Iter 35000:(163.1s) Lower Bound: -7643.4187 | Logprior: -508.8334 | Logll: -7379.5986 | Root Age: 231.7607\n",
      ">>> Iter 35000:(6.6s) Test Lower Bound: -7670.8701\n",
      "Iter 36000:(163.0s) Lower Bound: -7642.7292 | Logprior: -501.4060 | Logll: -7374.4204 | Root Age: 231.5032\n",
      "Iter 37000:(172.2s) Lower Bound: -7641.4392 | Logprior: -503.6097 | Logll: -7383.9194 | Root Age: 232.1591\n",
      "Iter 38000:(163.1s) Lower Bound: -7640.8885 | Logprior: -496.6266 | Logll: -7364.3560 | Root Age: 231.1580\n",
      "Iter 39000:(162.3s) Lower Bound: -7640.2248 | Logprior: -485.0327 | Logll: -7376.9487 | Root Age: 232.5778\n",
      "Iter 40000:(163.1s) Lower Bound: -7639.5070 | Logprior: -492.3342 | Logll: -7378.2803 | Root Age: 231.5388\n",
      ">>> Iter 40000:(6.5s) Test Lower Bound: -7664.3076\n",
      "Iter 41000:(162.9s) Lower Bound: -7638.2405 | Logprior: -503.5095 | Logll: -7398.3208 | Root Age: 232.1204\n",
      "Iter 42000:(162.1s) Lower Bound: -7638.0529 | Logprior: -502.8067 | Logll: -7414.6475 | Root Age: 232.2223\n",
      "Iter 43000:(162.8s) Lower Bound: -7637.5075 | Logprior: -492.2388 | Logll: -7354.2100 | Root Age: 231.3662\n",
      "Iter 44000:(162.7s) Lower Bound: -7637.0355 | Logprior: -489.1970 | Logll: -7370.4136 | Root Age: 230.7169\n",
      "Iter 45000:(162.8s) Lower Bound: -7636.5505 | Logprior: -517.5406 | Logll: -7377.5566 | Root Age: 231.5634\n",
      ">>> Iter 45000:(6.7s) Test Lower Bound: -7654.6123\n",
      "Iter 46000:(162.4s) Lower Bound: -7635.9848 | Logprior: -510.0287 | Logll: -7387.6772 | Root Age: 232.1068\n",
      "Iter 47000:(162.7s) Lower Bound: -7635.4222 | Logprior: -488.2578 | Logll: -7374.9531 | Root Age: 231.5099\n",
      "Iter 48000:(162.9s) Lower Bound: -7634.9827 | Logprior: -489.0120 | Logll: -7369.0674 | Root Age: 232.5177\n",
      "Iter 49000:(162.6s) Lower Bound: -7634.4686 | Logprior: -490.1121 | Logll: -7357.1348 | Root Age: 232.2520\n",
      "Iter 50000:(163.4s) Lower Bound: -7634.1411 | Logprior: -497.5703 | Logll: -7357.8394 | Root Age: 232.0987\n",
      ">>> Iter 50000:(6.6s) Test Lower Bound: -7651.4434\n",
      "Iter 51000:(161.9s) Lower Bound: -7634.0593 | Logprior: -481.1544 | Logll: -7368.9097 | Root Age: 231.2213\n",
      "Iter 52000:(162.4s) Lower Bound: -7633.4962 | Logprior: -506.8669 | Logll: -7371.8550 | Root Age: 231.0899\n",
      "Iter 53000:(162.0s) Lower Bound: -7633.4510 | Logprior: -499.9985 | Logll: -7383.7012 | Root Age: 231.3570\n",
      "Iter 54000:(162.6s) Lower Bound: -7632.9465 | Logprior: -500.9464 | Logll: -7374.6128 | Root Age: 231.4786\n",
      "Iter 55000:(163.3s) Lower Bound: -7632.8179 | Logprior: -493.9654 | Logll: -7359.8315 | Root Age: 231.2230\n",
      ">>> Iter 55000:(6.5s) Test Lower Bound: -7649.8726\n",
      "Iter 56000:(162.0s) Lower Bound: -7632.4588 | Logprior: -492.2606 | Logll: -7409.5825 | Root Age: 231.8178\n",
      "Iter 57000:(161.7s) Lower Bound: -7632.3899 | Logprior: -500.6789 | Logll: -7349.6763 | Root Age: 230.7872\n",
      "Iter 58000:(162.5s) Lower Bound: -7631.9281 | Logprior: -488.3924 | Logll: -7358.5991 | Root Age: 231.5512\n",
      "Iter 59000:(161.8s) Lower Bound: -7631.9098 | Logprior: -489.8987 | Logll: -7366.6401 | Root Age: 232.3044\n",
      "Iter 60000:(162.0s) Lower Bound: -7631.8772 | Logprior: -500.0750 | Logll: -7379.0332 | Root Age: 231.1616\n",
      ">>> Iter 60000:(6.5s) Test Lower Bound: -7645.9771\n",
      "Iter 61000:(162.0s) Lower Bound: -7631.6205 | Logprior: -492.5102 | Logll: -7352.5713 | Root Age: 232.3910\n",
      "Iter 62000:(163.2s) Lower Bound: -7631.4616 | Logprior: -518.9067 | Logll: -7378.4146 | Root Age: 230.6365\n",
      "Iter 63000:(162.0s) Lower Bound: -7631.0983 | Logprior: -509.8586 | Logll: -7354.7856 | Root Age: 230.9677\n",
      "Iter 64000:(162.2s) Lower Bound: -7631.0469 | Logprior: -500.5643 | Logll: -7376.7603 | Root Age: 231.9910\n",
      "Iter 65000:(162.0s) Lower Bound: -7630.9080 | Logprior: -475.1968 | Logll: -7353.5771 | Root Age: 231.1657\n",
      ">>> Iter 65000:(6.6s) Test Lower Bound: -7644.4893\n",
      "Iter 66000:(161.7s) Lower Bound: -7630.6351 | Logprior: -483.4478 | Logll: -7353.2974 | Root Age: 230.5550\n",
      "Iter 67000:(161.7s) Lower Bound: -7630.5116 | Logprior: -496.7332 | Logll: -7362.9385 | Root Age: 231.2090\n",
      "Iter 68000:(161.9s) Lower Bound: -7630.3581 | Logprior: -503.4282 | Logll: -7352.5977 | Root Age: 231.3127\n",
      "Iter 69000:(162.4s) Lower Bound: -7630.3258 | Logprior: -491.6397 | Logll: -7355.5376 | Root Age: 231.2108\n",
      "Iter 70000:(162.3s) Lower Bound: -7630.2578 | Logprior: -491.0009 | Logll: -7375.2095 | Root Age: 230.8981\n",
      ">>> Iter 70000:(6.5s) Test Lower Bound: -7643.8135\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 71000:(161.9s) Lower Bound: -7629.9038 | Logprior: -493.2608 | Logll: -7359.3008 | Root Age: 231.2656\n",
      "Iter 72000:(162.2s) Lower Bound: -7630.2013 | Logprior: -505.1701 | Logll: -7344.5562 | Root Age: 231.1418\n",
      "Iter 73000:(162.3s) Lower Bound: -7630.0860 | Logprior: -487.2344 | Logll: -7342.3594 | Root Age: 230.8121\n",
      "Iter 74000:(162.7s) Lower Bound: -7629.8805 | Logprior: -510.3177 | Logll: -7361.2998 | Root Age: 231.6403\n",
      "Iter 75000:(162.3s) Lower Bound: -7629.4953 | Logprior: -496.5406 | Logll: -7357.2671 | Root Age: 231.7545\n",
      ">>> Iter 75000:(6.5s) Test Lower Bound: -7644.0215\n",
      "Iter 76000:(162.6s) Lower Bound: -7629.6981 | Logprior: -498.8295 | Logll: -7349.1025 | Root Age: 230.5377\n",
      "Iter 77000:(162.6s) Lower Bound: -7629.4963 | Logprior: -498.1212 | Logll: -7392.4292 | Root Age: 231.4242\n",
      "Iter 78000:(161.4s) Lower Bound: -7629.5855 | Logprior: -505.1074 | Logll: -7372.7119 | Root Age: 231.7245\n",
      "Iter 79000:(161.7s) Lower Bound: -7629.5403 | Logprior: -494.3390 | Logll: -7367.4458 | Root Age: 230.7300\n",
      "Iter 80000:(162.1s) Lower Bound: -7629.3567 | Logprior: -508.5142 | Logll: -7366.0298 | Root Age: 231.2303\n",
      ">>> Iter 80000:(6.6s) Test Lower Bound: -7641.2729\n",
      "Iter 81000:(162.6s) Lower Bound: -7629.0357 | Logprior: -488.4519 | Logll: -7357.3447 | Root Age: 231.4041\n",
      "Iter 82000:(161.7s) Lower Bound: -7628.9538 | Logprior: -485.0546 | Logll: -7380.2080 | Root Age: 230.7546\n",
      "Iter 83000:(161.5s) Lower Bound: -7628.9445 | Logprior: -490.8369 | Logll: -7345.5449 | Root Age: 230.9799\n",
      "Iter 84000:(161.4s) Lower Bound: -7628.9300 | Logprior: -495.5671 | Logll: -7373.8984 | Root Age: 231.5220\n",
      "Iter 85000:(161.8s) Lower Bound: -7628.9543 | Logprior: -500.4485 | Logll: -7373.9819 | Root Age: 231.0343\n",
      ">>> Iter 85000:(6.5s) Test Lower Bound: -7640.3867\n",
      "Iter 86000:(161.4s) Lower Bound: -7628.7643 | Logprior: -502.8125 | Logll: -7353.9014 | Root Age: 231.1172\n",
      "Iter 87000:(161.6s) Lower Bound: -7628.9279 | Logprior: -529.0870 | Logll: -7352.4019 | Root Age: 231.1706\n",
      "Iter 88000:(161.7s) Lower Bound: -7628.8170 | Logprior: -488.8431 | Logll: -7365.2314 | Root Age: 231.1032\n",
      "Iter 89000:(162.5s) Lower Bound: -7628.8057 | Logprior: -496.3645 | Logll: -7359.3604 | Root Age: 230.9261\n",
      "Iter 90000:(161.8s) Lower Bound: -7628.6498 | Logprior: -497.7448 | Logll: -7366.6230 | Root Age: 231.2025\n",
      ">>> Iter 90000:(6.5s) Test Lower Bound: -7640.7183\n",
      "Iter 91000:(162.6s) Lower Bound: -7628.5673 | Logprior: -483.4451 | Logll: -7371.4956 | Root Age: 229.7249\n",
      "Iter 92000:(162.6s) Lower Bound: -7628.5344 | Logprior: -489.5488 | Logll: -7346.9688 | Root Age: 231.4923\n",
      "Iter 93000:(162.0s) Lower Bound: -7628.6590 | Logprior: -505.7679 | Logll: -7352.1147 | Root Age: 231.4105\n",
      "Iter 94000:(161.4s) Lower Bound: -7628.5748 | Logprior: -496.6418 | Logll: -7349.6104 | Root Age: 232.1599\n",
      "Iter 95000:(161.9s) Lower Bound: -7628.8532 | Logprior: -509.7205 | Logll: -7340.2285 | Root Age: 231.2498\n",
      ">>> Iter 95000:(6.6s) Test Lower Bound: -7639.8770\n",
      "Iter 96000:(161.6s) Lower Bound: -7628.5048 | Logprior: -495.5580 | Logll: -7363.8706 | Root Age: 230.7173\n",
      "Iter 97000:(162.6s) Lower Bound: -7628.4235 | Logprior: -504.0695 | Logll: -7351.2637 | Root Age: 231.0357\n",
      "Iter 98000:(162.6s) Lower Bound: -7628.4892 | Logprior: -505.1912 | Logll: -7359.5791 | Root Age: 231.6480\n",
      "Iter 99000:(162.0s) Lower Bound: -7628.4889 | Logprior: -488.6393 | Logll: -7391.6943 | Root Age: 231.6745\n",
      "Iter 100000:(161.5s) Lower Bound: -7628.5817 | Logprior: -494.1156 | Logll: -7356.8369 | Root Age: 231.4134\n",
      ">>> Iter 100000:(6.5s) Test Lower Bound: -7640.3301\n",
      "Iter 101000:(161.9s) Lower Bound: -7628.3369 | Logprior: -502.0077 | Logll: -7350.0630 | Root Age: 231.9255\n",
      "Iter 102000:(162.2s) Lower Bound: -7628.2603 | Logprior: -503.5062 | Logll: -7355.7627 | Root Age: 231.2727\n",
      "Iter 103000:(162.4s) Lower Bound: -7628.3066 | Logprior: -494.1783 | Logll: -7363.4736 | Root Age: 230.5928\n",
      "Iter 104000:(162.2s) Lower Bound: -7628.2598 | Logprior: -499.8183 | Logll: -7383.0918 | Root Age: 231.0602\n",
      "Iter 105000:(162.0s) Lower Bound: -7628.0845 | Logprior: -494.2985 | Logll: -7359.5513 | Root Age: 230.6636\n",
      ">>> Iter 105000:(6.5s) Test Lower Bound: -7639.6655\n",
      "Iter 106000:(162.1s) Lower Bound: -7628.1364 | Logprior: -510.9932 | Logll: -7354.1597 | Root Age: 230.4902\n",
      "Iter 107000:(161.9s) Lower Bound: -7627.8790 | Logprior: -492.3604 | Logll: -7357.1699 | Root Age: 231.2729\n",
      "Iter 108000:(162.1s) Lower Bound: -7627.9676 | Logprior: -522.7899 | Logll: -7337.2173 | Root Age: 231.5221\n",
      "Iter 109000:(161.8s) Lower Bound: -7628.1719 | Logprior: -485.0158 | Logll: -7344.2988 | Root Age: 230.6819\n",
      "Iter 110000:(162.6s) Lower Bound: -7627.9440 | Logprior: -489.5177 | Logll: -7348.7280 | Root Age: 230.7314\n",
      ">>> Iter 110000:(6.6s) Test Lower Bound: -7638.7622\n",
      "Iter 111000:(161.7s) Lower Bound: -7628.0525 | Logprior: -499.1592 | Logll: -7348.4033 | Root Age: 231.2722\n",
      "Iter 112000:(161.9s) Lower Bound: -7627.8582 | Logprior: -486.5005 | Logll: -7362.3350 | Root Age: 230.8777\n",
      "Iter 113000:(162.1s) Lower Bound: -7627.8873 | Logprior: -491.9319 | Logll: -7363.3955 | Root Age: 230.7625\n",
      "Iter 114000:(162.2s) Lower Bound: -7627.8883 | Logprior: -495.8082 | Logll: -7351.8813 | Root Age: 230.7015\n",
      "Iter 115000:(162.1s) Lower Bound: -7627.9494 | Logprior: -506.6639 | Logll: -7357.0532 | Root Age: 230.3958\n",
      ">>> Iter 115000:(6.5s) Test Lower Bound: -7638.2271\n",
      "Iter 116000:(161.3s) Lower Bound: -7627.7922 | Logprior: -502.0563 | Logll: -7370.9609 | Root Age: 230.9988\n",
      "Iter 117000:(161.6s) Lower Bound: -7628.0760 | Logprior: -509.2075 | Logll: -7360.6118 | Root Age: 230.2220\n",
      "Iter 118000:(161.8s) Lower Bound: -7627.8182 | Logprior: -491.8719 | Logll: -7346.8389 | Root Age: 231.7009\n",
      "Iter 119000:(162.1s) Lower Bound: -7627.9284 | Logprior: -506.5088 | Logll: -7354.3062 | Root Age: 230.8821\n",
      "Iter 120000:(163.0s) Lower Bound: -7627.7596 | Logprior: -487.4134 | Logll: -7343.1602 | Root Age: 229.8440\n",
      ">>> Iter 120000:(6.5s) Test Lower Bound: -7638.0415\n",
      "Iter 121000:(161.3s) Lower Bound: -7627.6560 | Logprior: -503.9926 | Logll: -7355.0977 | Root Age: 231.0961\n",
      "Iter 122000:(160.5s) Lower Bound: -7627.7255 | Logprior: -513.1376 | Logll: -7367.2969 | Root Age: 231.3725\n",
      "Iter 123000:(160.5s) Lower Bound: -7627.7869 | Logprior: -496.1857 | Logll: -7359.7720 | Root Age: 230.8231\n",
      "Iter 124000:(161.5s) Lower Bound: -7627.6605 | Logprior: -486.9229 | Logll: -7351.0986 | Root Age: 231.5294\n",
      "Iter 125000:(160.5s) Lower Bound: -7627.8471 | Logprior: -495.1885 | Logll: -7351.5112 | Root Age: 231.3646\n",
      ">>> Iter 125000:(6.6s) Test Lower Bound: -7638.7656\n",
      "Iter 126000:(161.1s) Lower Bound: -7627.6825 | Logprior: -498.7136 | Logll: -7348.7544 | Root Age: 231.4678\n",
      "Iter 127000:(167.2s) Lower Bound: -7627.8578 | Logprior: -493.1040 | Logll: -7356.5864 | Root Age: 231.0021\n",
      "Iter 128000:(165.1s) Lower Bound: -7627.6056 | Logprior: -511.3414 | Logll: -7377.4189 | Root Age: 231.5651\n",
      "Iter 129000:(169.0s) Lower Bound: -7627.5317 | Logprior: -506.6848 | Logll: -7361.8828 | Root Age: 231.5844\n",
      "Iter 130000:(166.2s) Lower Bound: -7627.7357 | Logprior: -493.5832 | Logll: -7338.9565 | Root Age: 231.6446\n",
      ">>> Iter 130000:(6.5s) Test Lower Bound: -7637.5078\n",
      "Iter 131000:(160.9s) Lower Bound: -7627.7238 | Logprior: -524.0216 | Logll: -7366.5615 | Root Age: 231.4284\n",
      "Iter 132000:(160.4s) Lower Bound: -7627.7696 | Logprior: -507.7968 | Logll: -7352.8896 | Root Age: 231.4687\n",
      "Iter 133000:(160.4s) Lower Bound: -7627.5792 | Logprior: -501.7781 | Logll: -7346.0933 | Root Age: 231.4707\n",
      "Iter 134000:(160.4s) Lower Bound: -7627.5235 | Logprior: -486.7023 | Logll: -7351.4170 | Root Age: 230.5879\n",
      "Iter 135000:(160.5s) Lower Bound: -7627.7114 | Logprior: -501.3450 | Logll: -7374.2202 | Root Age: 231.6397\n",
      ">>> Iter 135000:(6.4s) Test Lower Bound: -7637.6177\n",
      "Iter 136000:(161.2s) Lower Bound: -7627.6737 | Logprior: -481.8359 | Logll: -7371.6538 | Root Age: 231.3264\n",
      "Iter 137000:(160.8s) Lower Bound: -7627.7374 | Logprior: -498.6835 | Logll: -7368.5713 | Root Age: 231.1383\n",
      "Iter 138000:(160.6s) Lower Bound: -7627.6741 | Logprior: -508.1761 | Logll: -7363.9668 | Root Age: 230.8047\n",
      "Iter 139000:(160.4s) Lower Bound: -7627.6203 | Logprior: -512.3735 | Logll: -7369.1558 | Root Age: 231.4910\n",
      "Iter 140000:(160.7s) Lower Bound: -7627.6145 | Logprior: -503.8876 | Logll: -7347.0034 | Root Age: 230.6426\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">>> Iter 140000:(6.5s) Test Lower Bound: -7637.9097\n",
      "Iter 141000:(160.4s) Lower Bound: -7627.5459 | Logprior: -500.2725 | Logll: -7366.3091 | Root Age: 231.1342\n",
      "Iter 142000:(160.8s) Lower Bound: -7627.4465 | Logprior: -495.6278 | Logll: -7344.1426 | Root Age: 231.4919\n",
      "Iter 143000:(160.9s) Lower Bound: -7627.5472 | Logprior: -497.7082 | Logll: -7347.4722 | Root Age: 230.9996\n",
      "Iter 144000:(160.5s) Lower Bound: -7627.3580 | Logprior: -524.3526 | Logll: -7353.8760 | Root Age: 229.9201\n",
      "Iter 145000:(161.4s) Lower Bound: -7627.4496 | Logprior: -506.0348 | Logll: -7355.9546 | Root Age: 230.5397\n",
      ">>> Iter 145000:(6.4s) Test Lower Bound: -7637.2168\n",
      "Iter 146000:(160.5s) Lower Bound: -7627.4511 | Logprior: -503.9289 | Logll: -7346.7021 | Root Age: 230.6755\n",
      "Iter 147000:(160.5s) Lower Bound: -7627.3757 | Logprior: -499.1236 | Logll: -7364.0601 | Root Age: 230.2330\n",
      "Iter 148000:(160.5s) Lower Bound: -7627.3714 | Logprior: -483.7567 | Logll: -7367.0562 | Root Age: 230.9007\n",
      "Iter 149000:(160.5s) Lower Bound: -7627.3567 | Logprior: -490.3896 | Logll: -7354.0288 | Root Age: 231.8833\n",
      "Iter 150000:(160.8s) Lower Bound: -7627.2801 | Logprior: -514.0757 | Logll: -7356.1416 | Root Age: 231.6502\n",
      ">>> Iter 150000:(6.4s) Test Lower Bound: -7637.3037\n",
      "Iter 151000:(160.4s) Lower Bound: -7627.2889 | Logprior: -498.1958 | Logll: -7348.6470 | Root Age: 230.8508\n",
      "Iter 152000:(160.1s) Lower Bound: -7627.5559 | Logprior: -493.3813 | Logll: -7353.6108 | Root Age: 230.8690\n",
      "Iter 153000:(160.6s) Lower Bound: -7627.6168 | Logprior: -499.0013 | Logll: -7363.1138 | Root Age: 230.8227\n",
      "Iter 154000:(160.6s) Lower Bound: -7627.4377 | Logprior: -490.1387 | Logll: -7358.4609 | Root Age: 230.9146\n",
      "Iter 155000:(160.4s) Lower Bound: -7627.3589 | Logprior: -485.2886 | Logll: -7364.8433 | Root Age: 231.4036\n",
      ">>> Iter 155000:(6.5s) Test Lower Bound: -7637.6006\n",
      "Iter 156000:(160.4s) Lower Bound: -7627.2828 | Logprior: -498.3828 | Logll: -7356.8086 | Root Age: 231.6706\n",
      "Iter 157000:(161.7s) Lower Bound: -7627.4534 | Logprior: -495.5741 | Logll: -7359.2114 | Root Age: 231.5013\n",
      "Iter 158000:(160.2s) Lower Bound: -7627.4787 | Logprior: -502.7847 | Logll: -7340.2334 | Root Age: 231.1860\n",
      "Iter 159000:(160.7s) Lower Bound: -7627.3662 | Logprior: -509.8567 | Logll: -7338.9126 | Root Age: 231.2672\n",
      "Iter 160000:(161.0s) Lower Bound: -7627.4665 | Logprior: -498.0924 | Logll: -7346.5137 | Root Age: 231.3480\n",
      ">>> Iter 160000:(6.4s) Test Lower Bound: -7636.8027\n",
      "Iter 161000:(160.3s) Lower Bound: -7627.5277 | Logprior: -502.1515 | Logll: -7343.6250 | Root Age: 230.8031\n",
      "Iter 162000:(160.3s) Lower Bound: -7627.2648 | Logprior: -498.5085 | Logll: -7356.7778 | Root Age: 230.9580\n",
      "Iter 163000:(160.7s) Lower Bound: -7627.3254 | Logprior: -497.6039 | Logll: -7347.6104 | Root Age: 231.2908\n",
      "Iter 164000:(160.6s) Lower Bound: -7627.1393 | Logprior: -489.0295 | Logll: -7353.4956 | Root Age: 231.5462\n",
      "Iter 165000:(160.9s) Lower Bound: -7627.5328 | Logprior: -495.3406 | Logll: -7357.7500 | Root Age: 231.1848\n",
      ">>> Iter 165000:(6.4s) Test Lower Bound: -7637.1899\n",
      "Iter 166000:(160.5s) Lower Bound: -7627.2522 | Logprior: -506.4167 | Logll: -7350.6689 | Root Age: 231.1133\n",
      "Iter 167000:(160.4s) Lower Bound: -7627.1511 | Logprior: -498.8877 | Logll: -7366.4883 | Root Age: 231.4417\n",
      "Iter 168000:(160.6s) Lower Bound: -7627.5784 | Logprior: -505.3448 | Logll: -7352.1309 | Root Age: 231.4577\n",
      "Iter 169000:(160.9s) Lower Bound: -7627.4038 | Logprior: -497.4060 | Logll: -7340.7983 | Root Age: 231.5778\n",
      "Iter 170000:(160.8s) Lower Bound: -7627.2128 | Logprior: -506.6912 | Logll: -7350.3071 | Root Age: 231.7252\n",
      ">>> Iter 170000:(6.5s) Test Lower Bound: -7637.2476\n",
      "Iter 171000:(160.3s) Lower Bound: -7627.2601 | Logprior: -499.1617 | Logll: -7356.5757 | Root Age: 231.5216\n",
      "Iter 172000:(161.1s) Lower Bound: -7627.4979 | Logprior: -494.9306 | Logll: -7346.1123 | Root Age: 231.0991\n",
      "Iter 173000:(160.6s) Lower Bound: -7627.3394 | Logprior: -497.0512 | Logll: -7356.4771 | Root Age: 230.6694\n",
      "Iter 174000:(160.6s) Lower Bound: -7627.2045 | Logprior: -497.9583 | Logll: -7356.1431 | Root Age: 230.7339\n",
      "Iter 175000:(160.7s) Lower Bound: -7627.3739 | Logprior: -490.7696 | Logll: -7349.2578 | Root Age: 230.6340\n",
      ">>> Iter 175000:(6.4s) Test Lower Bound: -7637.7993\n",
      "Iter 176000:(160.6s) Lower Bound: -7627.0822 | Logprior: -509.5099 | Logll: -7367.8032 | Root Age: 230.9112\n",
      "Iter 177000:(160.5s) Lower Bound: -7627.3227 | Logprior: -502.7585 | Logll: -7349.9209 | Root Age: 231.0007\n",
      "Iter 178000:(160.8s) Lower Bound: -7627.3420 | Logprior: -511.6732 | Logll: -7347.0879 | Root Age: 230.8444\n",
      "Iter 179000:(160.5s) Lower Bound: -7627.2116 | Logprior: -499.3580 | Logll: -7339.3579 | Root Age: 231.1621\n",
      "Iter 180000:(160.5s) Lower Bound: -7627.2526 | Logprior: -506.2283 | Logll: -7360.5249 | Root Age: 231.1475\n",
      ">>> Iter 180000:(6.4s) Test Lower Bound: -7637.0322\n",
      "Iter 181000:(160.5s) Lower Bound: -7627.0898 | Logprior: -492.4859 | Logll: -7362.1572 | Root Age: 230.9484\n",
      "Iter 182000:(160.5s) Lower Bound: -7627.3225 | Logprior: -510.9291 | Logll: -7357.9507 | Root Age: 231.1711\n",
      "Iter 183000:(161.1s) Lower Bound: -7627.3996 | Logprior: -501.4386 | Logll: -7364.9053 | Root Age: 231.2108\n",
      "Iter 184000:(161.2s) Lower Bound: -7627.4110 | Logprior: -521.5304 | Logll: -7343.3862 | Root Age: 231.1247\n",
      "Iter 185000:(160.5s) Lower Bound: -7627.3474 | Logprior: -505.3718 | Logll: -7332.8228 | Root Age: 231.4922\n",
      ">>> Iter 185000:(6.5s) Test Lower Bound: -7636.5654\n",
      "Iter 186000:(161.4s) Lower Bound: -7627.0863 | Logprior: -493.3018 | Logll: -7341.0005 | Root Age: 231.2377\n",
      "Iter 187000:(162.1s) Lower Bound: -7627.3750 | Logprior: -487.4860 | Logll: -7349.0903 | Root Age: 231.2500\n",
      "Iter 188000:(160.3s) Lower Bound: -7627.1621 | Logprior: -500.3050 | Logll: -7348.3682 | Root Age: 231.0225\n",
      "Iter 189000:(160.5s) Lower Bound: -7627.2630 | Logprior: -506.8087 | Logll: -7338.8271 | Root Age: 231.1641\n",
      "Iter 190000:(160.6s) Lower Bound: -7627.0025 | Logprior: -499.6472 | Logll: -7344.1367 | Root Age: 231.5906\n",
      ">>> Iter 190000:(6.4s) Test Lower Bound: -7636.9873\n",
      "Iter 191000:(160.5s) Lower Bound: -7627.1112 | Logprior: -506.8625 | Logll: -7361.4971 | Root Age: 231.1296\n",
      "Iter 192000:(160.3s) Lower Bound: -7627.3077 | Logprior: -494.5694 | Logll: -7347.4326 | Root Age: 230.7306\n",
      "Iter 193000:(160.7s) Lower Bound: -7627.0872 | Logprior: -509.5860 | Logll: -7379.5942 | Root Age: 230.8663\n",
      "Iter 194000:(161.2s) Lower Bound: -7626.8573 | Logprior: -504.9041 | Logll: -7350.0527 | Root Age: 231.2547\n",
      "Iter 195000:(160.9s) Lower Bound: -7627.2693 | Logprior: -500.4894 | Logll: -7333.4541 | Root Age: 230.6217\n",
      ">>> Iter 195000:(6.4s) Test Lower Bound: -7636.5454\n",
      "Iter 196000:(160.5s) Lower Bound: -7626.9893 | Logprior: -499.8200 | Logll: -7343.4980 | Root Age: 230.8232\n",
      "Iter 197000:(160.5s) Lower Bound: -7627.1152 | Logprior: -497.8738 | Logll: -7352.3159 | Root Age: 231.1077\n",
      "Iter 198000:(160.8s) Lower Bound: -7627.1511 | Logprior: -490.3508 | Logll: -7354.1245 | Root Age: 230.9384\n",
      "Iter 199000:(162.2s) Lower Bound: -7627.1789 | Logprior: -501.2636 | Logll: -7337.5752 | Root Age: 230.9669\n",
      "Iter 200000:(160.3s) Lower Bound: -7627.1004 | Logprior: -511.0046 | Logll: -7345.6250 | Root Age: 230.9669\n",
      ">>> Iter 200000:(6.5s) Test Lower Bound: -7636.8232\n"
     ]
    }
   ],
   "source": [
    "# run vbpi.\n",
    "# the trained model will be saved to 'results/HCV/mcmc_vimco_10_psp_skyride_fixed_rate.pt'.\n",
    "test_lb, _ = model.learn(0.001, maxiter=200000, n_particles=10, warm_start_interval=50000, method='vimco',\n",
    "                         save_to_path='../results/HCV/mcmc_vimco_10_psp_skyride_fixed_rate.pt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample (log) effective population trajectory \n",
    "vbpi_skyride_pop_traj = model.sample_pop_traj(cut_off=250, n_traj=75000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[ 138.8319,  138.8319,  138.8319,  ..., 3861.5222, 3861.5222,\n",
       "         4277.7002],\n",
       "        [ 113.7947,  113.7947,  113.7947,  ..., 4056.6475, 4056.6475,\n",
       "         9845.8291],\n",
       "        [  87.8090,   87.8090,   87.8090,  ..., 2786.1008, 7535.8096,\n",
       "         4744.2803],\n",
       "        ...,\n",
       "        [ 145.4557,  145.4557,  145.4557,  ..., 4610.2266, 4610.2266,\n",
       "         5466.1299],\n",
       "        [  39.2343,   39.2343,   39.2343,  ..., 5445.3076, 4209.6978,\n",
       "         4070.5056],\n",
       "        [ 134.4723,  134.4723,  134.4723,  ..., 6724.9614, 2584.2993,\n",
       "         2584.2993]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vbpi_skyride_pop_traj.exp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample tree height\n",
    "vbpi_tree_height = model.sample_tree_height(75000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([232.4026, 231.2287, 231.0143,  ..., 228.3037, 231.1365, 241.8365])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vbpi_tree_height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample tree loglikelihood\n",
    "vbpi_tree_logll = model.sample_tree_loglikelihood(75000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([-7353.8950, -7355.3555, -7345.7427,  ..., -7341.6108, -7349.6147,\n",
       "        -7344.8442])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vbpi_tree_logll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}