Merge pull request #127 from joezuntz/fix-list-thin-burn

Fix thinning and burning in list sampler

cosmosis/__init__.py

@@ -9,3 +9,4 @@
 from .gaussian_likelihood import GaussianLikelihood, \
                                  SingleValueGaussianLikelihood, \
                                  WindowedGaussianLikelihood
+from .output import TextColumnOutput, InMemoryOutput, NullOutput, AstropyOutput, CosmoMCOutput, FitsOutput

cosmosis/main.py

@@ -430,7 +430,7 @@ def run_cosmosis(ini, pool=None, pipeline=None, values=None, priors=None, overri
             logs.overview(f"Successful posterior evaluations = {run_count_ok_total} across all processes")
             if output:
                 output.final("evaluations", run_count_total)
-                output.final("successes", run_count_total)
+                output.final("successes", run_count_ok_total)
                 output.final("complete", "1")
 
 

cosmosis/samplers/list/list_sampler.py

@@ -26,7 +26,13 @@ def config(self):
         self.converged = False
         self.filename = self.read_ini("filename", str)
         self.save_name = self.read_ini("save", str, "")
-        self.burn = self.read_ini("burn", int, 0)
+        # The burn parameter can be an int or a float
+        # so we read it as a string and then convert it
+        burn = self.read_ini("burn", str, "0")
+        try:
+            self.burn = int(burn)
+        except ValueError:
+            self.burn = float(burn)
         self.thin = self.read_ini("thin", int, 1)
         limits = self.read_ini("limits", bool, False)
         self.chunk_size = self.read_ini("chunk_size", int, 100)
@@ -57,6 +63,12 @@ def execute(self):
         file_options = {"filename":self.filename}
         column_names, samples, _, _, _ = TextColumnOutput.load_from_options(file_options)
         samples = samples[0]
+        if (self.burn != 0) or (self.thin != 1):
+            if isinstance(self.burn, float):
+                burn = int(self.burn*len(samples))
+            else:
+                burn = self.burn
+            samples = samples[burn::self.thin]
         # find where in the parameter vector of the pipeline
         # each of the table parameters can be found
         replaced_params = []

cosmosis/samplers/minuit/minuit_sampler.py

@@ -28,15 +28,16 @@ class MinuitOptionsStruct(ct.Structure):
 class MinuitSampler(ParallelSampler):
     needs_output = True
     needs_parallel_output = False
+    libminuit_name = libname
     libminuit = None
     sampler_outputs = [("prior", float), ("post", float)]
 
     def config(self):
         self.converged = False
         if MinuitSampler.libminuit is None:
-            if not os.path.exists(libname):
+            if not os.path.exists(self.libminuit_name):
                 raise ValueError("The CosmoSIS minuit2 wrapper was not compiled.  If you installed CosmoSIS manually you need the library minuit2 to use this sampler. See the wiki for more details.")
-            MinuitSampler.libminuit = ct.cdll.LoadLibrary(libname)
+            MinuitSampler.libminuit = ct.cdll.LoadLibrary(self.libminuit_name)
         self.maxiter = self.read_ini("maxiter", int, 1000)
         self.save_dir = self.read_ini("save_dir", str, "")
         self.save_cov = self.read_ini("save_cov", str, "")

cosmosis/samplers/sampler.py

@@ -154,7 +154,7 @@ def read_resume_info(self):
     @classmethod
     def get_sampler(cls, name):
         try:
-            cls.registry[name.lower()]
+            return cls.registry[name.lower()]
         except KeyError:
             raise KeyError(f"Unknown sampler {name}")
 

cosmosis/test/test_samplers.py

@@ -1,8 +1,4 @@
-from cosmosis.runtime.config import Inifile
-from cosmosis.samplers.sampler import Sampler
-import cosmosis.samplers.minuit.minuit_sampler
-from cosmosis.runtime.pipeline import LikelihoodPipeline
-from cosmosis.output.in_memory_output import InMemoryOutput
+from cosmosis import Inifile, Sampler, LikelihoodPipeline, InMemoryOutput, TextColumnOutput
 from cosmosis.postprocessing import postprocessor_for_sampler
 import tempfile
 import os
@@ -11,7 +7,8 @@
 import numpy as np
 from astropy.table import Table
 
-minuit_compiled = os.path.exists(cosmosis.samplers.minuit.minuit_sampler.libname)
+
+minuit_compiled = os.path.exists(Sampler.get_sampler("minuit").libminuit_name)
 
 def run(name, check_prior, check_extra=True, can_postprocess=True, do_truth=False, no_extra=False, pp_extra=True, pp_2d=True, **options):
 
@@ -175,6 +172,37 @@ def test_star():
 def test_test():
     run('test', False, can_postprocess=False)
 
+def test_list_sampler():
+    # test that the burn and thin parameters work
+    # make a mock output file with some samples
+    with tempfile.TemporaryDirectory() as dirname:
+
+        input_chain_filename = os.path.join(dirname, "input_chain.txt")
+        input_chain = TextColumnOutput(input_chain_filename)
+        input_chain.add_column('parameters--p1', float)
+        input_chain.add_column('parameters--p2', float)
+        nrow = 100
+        data = np.random.normal(size=(nrow, 2)).clip(-2.99, 2.99)
+        for i in range(nrow):
+            input_chain.parameters(data[i, 0], data[i, 1])
+        input_chain.close()
+
+        burn = 0.2
+        thin = 2
+
+        data = data[20::thin]
+
+        output = run("list", True, can_postprocess=False, filename=input_chain_filename, burn=burn, thin=thin)
+        p1 = output['parameters--p1']
+        p2 = output['parameters--p2']
+        p3 = p1 + p2
+        expected_like = -(p1**2 + p2**2)/2.
+        assert p1.size == 40
+        assert np.allclose(p1, data[:, 0])
+        assert np.allclose(p2, data[:, 1])
+        assert np.allclose(output['PARAMETERS--P3'], p3)
+        assert np.allclose(output["post"] - output["prior"], expected_like)
+
 @pytest.mark.skipif(sys.version_info < (3,7), reason="pocomc requires python3.6+")
 def test_poco():
     try: