docs/examples/Tutorials/Prep_Rest_Post_Ana/posteriorsampler.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PosteriorSampler\n",
    "==============\n",
    "\n",
    "<hr style=\"height:2.5px\">\n",
    "\n",
    "In this tutorial, we will perform sampling when given an ensemble. Previously, in the [Ensemble](https://biceps.readthedocs.io/en/latest/examples/Tutorials/Prep_Rest_Post_Ana/ensemble.html) tutorial we constructed ensembles for two lambda values and saved them to pickle files.  Please read the section [lambda values](https://biceps.readthedocs.io/en/latest/examples/Tutorials/Prep_Rest_Post_Ana/restraint.ipynb#Quick-note-on-lambda-values:) in the previous tutorial for more information.\n",
    "\n",
    "<hr style=\"height:2.5px\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle, os\n",
    "import numpy as np\n",
    "import biceps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = 'results'\n",
    "n_lambdas = 2\n",
    "lambda_values = np.linspace(0.0, 1.0, n_lambdas)\n",
    "ensembles = []\n",
    "for lam in lambda_values:\n",
    "    with open(outdir+\"/ensemble_%s.pkl\"%lam,'rb') as file:\n",
    "        ensemble = pickle.load(file)\n",
    "    ensembles.append(ensemble)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to specify number of steps for BICePs sampling. We recommend to run at least 1M steps for converged Monte Carlo samplings.  Checking the convergence of the MCMC trajectory can be done simply using the submodule, `biceps.Convergence`.\n",
    "\n",
    "\n",
    "Using the first ensemble of the list of ensembles...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsteps = 100000 # number of steps of MCMC simulation\n",
    "sampler = biceps.PosteriorSampler(ensembles[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we sample for 100 k steps..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step\t\tState\tPara Indices\t\tAvg Energy\tAcceptance (%)\n",
      "0\t\t[87]\t[196, 127, 195]\t\t7.150\t\t56.94\tFalse\n",
      "10000\t\t[67]\t[235, 150, 184]\t\t8.713\t\t69.42\tFalse\n",
      "20000\t\t[67]\t[267, 141, 202]\t\t8.858\t\t73.08\tTrue\n",
      "30000\t\t[68]\t[243, 137, 196]\t\t9.303\t\t72.81\tTrue\n",
      "40000\t\t[27]\t[229, 150, 199]\t\t7.451\t\t72.53\tTrue\n",
      "50000\t\t[25]\t[230, 134, 197]\t\t9.089\t\t72.14\tTrue\n",
      "60000\t\t[65]\t[276, 129, 182]\t\t9.437\t\t72.21\tFalse\n",
      "70000\t\t[18]\t[224, 147, 198]\t\t11.604\t\t72.58\tTrue\n",
      "80000\t\t[21]\t[302, 117, 182]\t\t11.992\t\t73.06\tTrue\n",
      "90000\t\t[79]\t[241, 135, 199]\t\t7.123\t\t73.16\tTrue\n",
      "\n",
      "Accepted 72.88811881188118 % \n",
      "\n",
      "\n",
      "Accepted [ ...Nuisance paramters..., state] %\n",
      "Accepted [32.68811881 31.30891089 31.30891089  8.89108911] % \n",
      "\n"
     ]
    }
   ],
   "source": [
    "sampler.sample(nsteps=nsteps, burn=1000, print_freq=10000, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, process and save the results for the next step, [Analysis](https://biceps.readthedocs.io/en/latest/examples/Tutorials/Prep_Rest_Post_Ana/analysis.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "lam = lambda_values[0]\n",
    "filename = os.path.join(outdir,'traj_lambda%2.2f.npz'%(lam))\n",
    "sampler.traj.process_results(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we will do the same thing for the second ensemble in the list of ensembles.  **Note that we separated them for simplicity of the tutorial.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 99%|██████████████████████████████████████████████████████████▍| 100000/101000 [00:01<00:00, 61468.68it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Accepted 65.55940594059406 % \n",
      "\n",
      "\n",
      "Accepted [ ...Nuisance paramters..., state] %\n",
      "Accepted [32.73762376 31.24950495 31.24950495  1.57227723] % \n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "lam = lambda_values[1]\n",
    "sampler = biceps.PosteriorSampler(ensembles[1])\n",
    "sampler.sample(nsteps=nsteps, burn=1000, verbose=False)\n",
    "filename = os.path.join(outdir,'traj_lambda%2.2f.npz'%(lam))\n",
    "sampler.traj.process_results(filename)\n",
    "biceps.toolbox.save_object(sampler, filename.replace(\".npz\", \".pkl\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 style=\"text-align: center;font-size: 18pt;\">Conclusion</h1>\n",
    "\n",
    "In this tutorial, we used the `biceps.PosteriorSampler` class to perform MCMC sampling given the ensemble Python object. In the [next tutorial](https://biceps.readthedocs.io/en/latest/examples/Tutorials/Prep_Rest_Post_Ana/analysis.html), we will analyze our trajectory data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h6 style=\"align: justify;font-size: 12pt\"># <span style=\"color:red;\">NOTE</span>: The following cell is for pretty notebook rendering</h6>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    #notebook { /* centers the page content */\n",
       "    padding-left:8%;\n",
       "    padding-right:9%;\n",
       "    }\n",
       "    #ipython-main-app {\n",
       "        position: relative;\n",
       "    }\n",
       "    #jupyter-main-app {\n",
       "        position: relative;\n",
       "    }\n",
       "    div.cell{\n",
       "        width:98%;\n",
       "        margin-left:4% !important;*/\n",
       "        margin-right:4%;\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: \"Charis SIL\", Palatino, serif;\n",
       "        text-align: center;\n",
       "    }\n",
       "    h4{\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "    }\n",
       "    p {\n",
       "        display: block;\n",
       "        text-align: justify;\n",
       "        word-wrap: break-word;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 145%;\n",
       "        font-size: 12pt;\n",
       "        width:99%;\n",
       "        margin-left:auto;\n",
       "        margin-right:auto;\n",
       "\ttext-align: justify;\n",
       "        word-wrap: break-word;\n",
       "    }\n",
       "    .output_png {\n",
       "        display: block;/*table-cell;*/\n",
       "        text-align: center;\n",
       "        align: center;\n",
       "        vertical-align: middle;\n",
       "        margin-left: auto;\n",
       "        margin-right: auto;\n",
       "        width: auto; /* width: 100%;*/\n",
       "    }\n",
       "    div.output_area pre {font-weight: normal; text-align: justify; }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\", source-code-pro,Consolas, monospace;\n",
       "    }\n",
       "    .prompt{\n",
       "        display: None;\n",
       "    }\n",
       "    .text_cell_render p {\n",
       "        display: block;\n",
       "        text-align: justify;\n",
       "        word-wrap: break-word;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #4057A1;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    /* override table width restrictions as found on https://github.com/getpelican/pelican/issues/1311 */\n",
       "    .wy-table-responsive table td, .wy-table-responsive table th {\n",
       "        /* !important prevents the common CSS stylesheets from\n",
       "           overriding this as on RTD they are loaded after this stylesheet */\n",
       "        white-space: normal !important;\n",
       "    }\n",
       "\n",
       "    .wy-table-responsive {\n",
       "        overflow: visible !important;\n",
       "    }\n",
       "</style>\n",
       "\n",
       "\n",
       "<!-- MathJax configuration -->\n",
       "<script type=\"text/x-mathjax-config\">\n",
       "    src=\"http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML\">\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}},\n",
       "                    linebreaks: { automatic: true }}\n",
       "        });\n",
       "</script>\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"../../../theme.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}