Minor improvements.

- Adding new sampler to the notebook with sampler comparions.
- Added support for weight samples.
- Fixed bug in ncm_stats_vec.c (first element with zero weight resulted in a nan).

.pylintrc

@@ -61,3 +61,4 @@ ignore-none=false
 
 # Expected format of line ending, e.g. empty (any line ending), LF or CRLF.
 expected-line-ending-format=LF
+max-line-length=88

notebooks/apes_tests/rosenbrock_simple.ipynb

@@ -70,7 +70,8 @@
     "import emcee\n",
     "import zeus\n",
     "from pyhmc import hmc\n",
-    "import pocomc"
+    "import pocomc\n",
+    "import nautilus"
    ]
   },
   {
@@ -95,9 +96,9 @@
     "Ncm.cfg_init()\n",
     "Ncm.cfg_set_log_handler(lambda msg: sys.stdout.write(msg) and sys.stdout.flush())\n",
     "\n",
-    "ssize = 600000\n",
+    "ssize = 3000000\n",
     "nwalkers = 300\n",
-    "burin_steps = 600\n",
+    "burin_steps = 1200\n",
     "verbose = False"
    ]
   },
@@ -171,7 +172,7 @@
    "id": "cb1b0cfa",
    "metadata": {},
    "source": [
-    "### Running Emcee\n",
+    "## Running Emcee\n",
     "\n",
     "In the cell below, we will run the Emcee algorithm with the same initial point `p0` generated previously. We will generate a chain of samples using Emcee and store the resulting chain in a `Catalog` object. This will allow us to apply the same tests to all algorithms and compare their performance on an even footing."
    ]
@@ -197,7 +198,7 @@
    "id": "f85f7bfc",
    "metadata": {},
    "source": [
-    "### Running Zeus\n",
+    "## Running Zeus\n",
     "\n",
     "In the cell below, we will run the Zeus algorithm with the same initial point `p0` generated previously. We will generate a chain of samples using Zeus and store the resulting chain in a `Catalog` object. One difference from Emcee is that Zeus output chains are not interweaved, so we need to inform the `Catalog` object of this fact to ensure that autocorrelation estimates are correct."
    ]
@@ -223,7 +224,7 @@
    "id": "5b046365-3f85-41ed-a0a6-7f6afc59300d",
    "metadata": {},
    "source": [
-    "### Running pyhmc: Hamiltonian Monte Carlo\n",
+    "## Running pyhmc: Hamiltonian Monte Carlo\n",
     "\n",
     "In the cell below, we will run the pyhmc algorithm with the same initial point `p0` generated previously, however since pyhmc is not an ensemble sampler we use only the first point. We will generate a chain of samples using pyhmc and store the resulting chain in a `Catalog` object. Moreover, here we need to remove a longer burn-in period of 1000 samples, since pyhmc is not an ensemble sampler."
    ]
@@ -243,10 +244,21 @@
     "    return_logp=True,\n",
     ")\n",
     "mcat_pyhmc = Catalog(ndim=ndim, nwalkers=1, run_type=\"PyHMC\")\n",
-    "mcat_pyhmc.add_points_m2lnp(chain_pyhmc, -2.0 * log_prob_zeus)\n",
+    "mcat_pyhmc.add_points_m2lnp(chain_pyhmc, -2.0 * log_prob_pyhmc)\n",
     "mcat_pyhmc.trim(ssize // 2)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "3c8f6efa-8c27-4e01-be7f-d14dee45bfce",
+   "metadata": {},
+   "source": [
+    "## Running PocoMC\n",
+    "\n",
+    "In the cell below, we will run the pocomc algorithm with the same initial point `p0` generated previously. \n",
+    "We will generate a chain of samples using PocoMC and store the resulting chain in a `Catalog` object. "
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -260,19 +272,15 @@
     "    else:\n",
     "        return 0.0\n",
     "\n",
-    "n_dim = 2\n",
-    "n_particles = nwalkers\n",
-    "prior_samples = np.random.uniform(size=(n_particles, n_dim), low=-100.0, high=100.0)\n",
-    "\n",
     "sampler = pocomc.Sampler(\n",
-    "    n_particles,\n",
-    "    n_dim,\n",
+    "    nwalkers,\n",
+    "    ndim,\n",
     "    lambda x: log_prob(x,()),\n",
     "    log_prior,\n",
     "    vectorize_likelihood=True,\n",
     "    bounds=(-100.0, 100.0)\n",
     ")\n",
-    "sampler.run(prior_samples)\n",
+    "sampler.run(p0)\n",
     "\n",
     "results = sampler.results\n"
    ]
@@ -345,7 +353,59 @@
    "outputs": [],
    "source": [
     "pocomc.plotting.trace(results, dims = [0, 1])\n",
-    "plt.show()\n"
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2cce7d77-c95b-4df4-b4f4-50cf1af83a08",
+   "metadata": {},
+   "source": [
+    "## Running Nautilus\n",
+    "\n",
+    "In the cell below, we will run the pocomc algorithm with the same initial point `p0` generated previously. \n",
+    "We will generate a chain of samples using Nautilus and store the resulting chain in a `Catalog` object. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9927f655-8e71-4e01-8c6e-10b582874b44",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prior = nautilus.Prior()\n",
+    "prior.add_parameter('x0', dist=(-100, +100))\n",
+    "prior.add_parameter('x1', dist=(-100, +100))\n",
+    "\n",
+    "neval=0\n",
+    "def nautilus_likelihood(param_dict):\n",
+    "    x = np.array([param_dict['x0'], param_dict['x1']])\n",
+    "    global neval\n",
+    "    neval = neval + 1\n",
+    "    return log_prob(x,())\n",
+    "\n",
+    "nautilus_sampler = nautilus.Sampler(prior, nautilus_likelihood, n_live=3000)\n",
+    "nautilus_sampler.run(verbose=False, n_eff=ssize/10.0)\n",
+    "print(f\"Total number of likelihood evaluations {neval}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7c27c292-290f-449d-989a-1207a83e2c7e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chain_nautilus, log_w_nautilus, log_prob_naitulus = nautilus_sampler.posterior()\n",
+    "\n",
+    "mcat_nautilus = Catalog(ndim=ndim, nwalkers=nwalkers, run_type=\"Nautilus\", weighted=True)\n",
+    "\n",
+    "weights=np.exp(log_w_nautilus - np.max(log_w_nautilus))\n",
+    "mcat_nautilus.add_points_m2lnp(chain_nautilus, -2.0 * log_prob_naitulus, weights=weights)\n",
+    "#mcat_nautilus.trim(ssize // 2)\n",
+    "\n",
+    "nautilus_sampler.print_status()"
    ]
   },
   {
@@ -379,7 +439,9 @@
     "Ncm.cfg_msg_sepa()\n",
     "mcat_pyhmc.print_status()\n",
     "Ncm.cfg_msg_sepa()\n",
-    "mcat_poco.print_status()"
+    "mcat_poco.print_status()\n",
+    "Ncm.cfg_msg_sepa()\n",
+    "mcat_nautilus.print_status()\n"
    ]
   },
   {
@@ -403,7 +465,8 @@
     "print(np.abs(mcat_emcee.get_mean() / mean - 1.0))\n",
     "print(np.abs(mcat_zeus.get_mean() / mean - 1.0))\n",
     "print(np.abs(mcat_pyhmc.get_mean() / mean - 1.0))\n",
-    "print(np.abs(mcat_poco.get_mean() / mean - 1.0))"
+    "print(np.abs(mcat_poco.get_mean() / mean - 1.0))\n",
+    "print(np.abs(mcat_nautilus.get_mean() / mean - 1.0))"
    ]
   },
   {
@@ -427,7 +490,8 @@
     "print(np.abs(np.sqrt(np.diagonal(mcat_emcee.get_covar())) / sigma - 1.0))\n",
     "print(np.abs(np.sqrt(np.diagonal(mcat_zeus.get_covar())) / sigma - 1.0))\n",
     "print(np.abs(np.sqrt(np.diagonal(mcat_pyhmc.get_covar())) / sigma - 1.0))\n",
-    "print(np.abs(np.sqrt(np.diagonal(mcat_poco.get_covar())) / sigma - 1.0))"
+    "print(np.abs(np.sqrt(np.diagonal(mcat_poco.get_covar())) / sigma - 1.0))\n",
+    "print(np.abs(np.sqrt(np.diagonal(mcat_nautilus.get_covar())) / sigma - 1.0))"
    ]
   },
   {
@@ -451,7 +515,8 @@
     "mcs_emcee = mcat_emcee.get_mcsamples(\"EMCEE\")\n",
     "mcs_zeus = mcat_zeus.get_mcsamples(\"ZEUS\")\n",
     "mcs_pyhmc = mcat_pyhmc.get_mcsamples(\"PyHMC\")\n",
-    "mcs_poco = mcat_poco.get_mcsamples(\"POCO\")"
+    "mcs_poco = mcat_poco.get_mcsamples(\"POCO\")\n",
+    "mcs_nautilus = mcat_nautilus.get_mcsamples(\"NAUTILUS\")"
    ]
   },
   {
@@ -475,8 +540,8 @@
     "\n",
     "g = getdist.plots.get_subplot_plotter(width_inch=plt.rcParams[\"figure.figsize\"][0])\n",
     "g.settings.linewidth = 0.01\n",
-    "g.triangle_plot([mcs_apes, mcs_emcee, mcs_zeus, mcs_pyhmc, mcs_poco], shaded=True)\n",
-    "#g.triangle_plot([mcs_apes, mcs_poco], shaded=True)"
+    "g.triangle_plot([mcs_apes, mcs_emcee, mcs_zeus, mcs_pyhmc, mcs_poco, mcs_nautilus], shaded=True)\n",
+    "g.triangle_plot([mcs_apes, mcs_nautilus], shaded=True)"
    ]
   }
  ],