Merge pull request #117 from pcubillos/hotfixes

Hotfixes

mc3/VERSION.py

@@ -4,6 +4,6 @@
 # MC3 Version:
 MC3_VER =  3  # Major version
 MC3_MIN =  0  # Minor version
-MC3_REV =  '0rc2'  # Revision
+MC3_REV =  0  # Revision
 
 __version__ = '{}.{}.{}'.format(MC3_VER, MC3_MIN, MC3_REV)

mc3/plots/plots.py

@@ -34,7 +34,7 @@
 
 
 def trace(posterior, zchain=None, pnames=None, thinning=1,
-          burnin=0, fignum=100, savefile=None, fmt=".", ms=2.5, fs=11):
+          burnin=0, fignum=1000, savefile=None, fmt=".", ms=2.5, fs=11):
   """
   Plot parameter trace MCMC sampling.
 
@@ -98,7 +98,7 @@ def trace(posterior, zchain=None, pnames=None, thinning=1,
   axes = []
   ipar = 0
   for page in range(npages):
-      fig = plt.figure(page, figsize=(8.5,11.0))
+      fig = plt.figure(fignum+page, figsize=(8.5,11.0))
       plt.clf()
       plt.subplots_adjust(left=0.15, right=0.95, bottom=0.05, top=0.97,
                           hspace=0.15)
@@ -143,7 +143,7 @@ def trace(posterior, zchain=None, pnames=None, thinning=1,
   return axes
 
 
-def histogram(posterior, pnames=None, thinning=1, fignum=300,
+def histogram(posterior, pnames=None, thinning=1, fignum=1100,
               savefile=None, bestp=None, quantile=None, pdf=None,
               xpdf=None, ranges=None, axes=None, lw=2.0, fs=11,
               # Deprecated: Remove by 2020-07-01
@@ -285,18 +285,18 @@ def histogram(posterior, pnames=None, thinning=1, fignum=300,
       ax.set_ylim(0, maxylim)
 
   if savefile is not None:
-      for j, fig in enumerate(figs):
+      for page, fig in enumerate(figs):
           if npages > 1:
               sf = os.path.splitext(savefile)
-              fig.savefig("{:s}_page{:02d}{:s}".format(sf[0], j+1, sf[1]),
+              fig.savefig("{:s}_page{:02d}{:s}".format(sf[0], page, sf[1]),
                           bbox_inches='tight')
           else:
               fig.savefig(savefile, bbox_inches='tight')
 
   return axes
 
 
-def pairwise(posterior, pnames=None, thinning=1, fignum=200,
+def pairwise(posterior, pnames=None, thinning=1, fignum=1200,
              savefile=None, bestp=None, nbins=35, nlevels=20,
              absolute_dens=False, ranges=None, fs=11, rect=None, margin=0.01):
   """
@@ -453,7 +453,7 @@ def pairwise(posterior, pnames=None, thinning=1, fignum=200,
 
 
 def rms(binsz, rms, stderr, rmslo, rmshi, cadence=None, binstep=1,
-        timepoints=[], ratio=False, fignum=410,
+        timepoints=[], ratio=False, fignum=1300,
         yran=None, xran=None, savefile=None):
   """
   Plot the RMS vs binsize curve.
@@ -543,7 +543,7 @@ def rms(binsz, rms, stderr, rmslo, rmshi, cadence=None, binstep=1,
 
 
 def modelfit(data, uncert, indparams, model, nbins=75,
-             fignum=411, savefile=None, fmt="."):
+             fignum=1400, savefile=None, fmt="."):
   """
   Plot the binned dataset with given uncertainties and model curves
   as a function of indparams.

mc3/sampler_driver.py

@@ -494,7 +494,8 @@ def sample(data=None, uncert=None, func=None, params=None, indparams=[],
           thinning, resume, log, **kwargs)
 
   if leastsq is not None:
-      if output['best_log_post']-fit_output['best_log_post'] > 3.0e-8:
+      if (output['best_log_post'] - fit_output['best_log_post'] > 5.0e-8
+          and np.any((output['bestp'] - fit_output['bestp'])!=0.0)):
           log.warning("MCMC found a better fit than the minimizer:\n"
               "MCMC best-fitting parameters:        (chisq={:.8g})\n{}\n"
               "Minimizer best-fitting parameters:   (chisq={:.8g})\n{}".

mc3/stats/stats.py

@@ -179,6 +179,7 @@ def chisq(model, data, uncert,
     Examples
     --------
     >>> import mc3.stats as ms
+    >>> import numpy as np
     >>> # Compute chi-squared for a given model fitting a data set:
     >>> data   = np.array([1.1, 1.2, 0.9, 1.0])
     >>> model  = np.array([1.0, 1.0, 1.0, 1.0])
@@ -243,20 +244,33 @@ def dwt_chisq(model, data, params, priors=None, priorlow=None, priorup=None):
     --------
     >>> import mc3.stats as ms
     >>> import numpy as np
-
+    >>> # Compute chi-squared for a given model fitting a data set:
     >>> data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0])
     >>> model = np.ones(8)
     >>> params = np.array([1.0, 0.1, 0.1])
-    >>> chisq = ms.chisq(model, data, params)
+    >>> chisq = ms.dwt_chisq(model, data, params)
     >>> print(chisq)
     1693.22308882
+    >>> # Now, say this is a three-parameter model, with a Gaussian prior
+    >>> # on the last parameter:
+    >>> priors = np.array([1.0, 0.2, 0.3])
+    >>> plow   = np.array([0.0, 0.0, 0.1])
+    >>> pup    = np.array([0.0, 0.0, 0.1])
+    >>> chisq = ms.dwt_chisq(model, data, params, priors, plow, pup)
+    >>> print(chisq)
+    1697.2230888243134
     """
+    if len(params) < 3:
+        with mu.Log() as log:
+            log.error('Wavelet chisq should have at least three parameters.')
+
     if priors is None or priorlow is None or priorup is None:
         return dwt.chisq(params, model, data)
 
     iprior = (priorlow > 0) & (priorup > 0)
     dprior = (params - priors)[iprior]
-    return dwt.chisq(params, model, data, dprior, priorlow, priorup)
+    return dwt.chisq(params, model, data, dprior,
+                     priorlow[iprior], priorup[iprior])
 
 
 def log_prior(posterior, prior, priorlow, priorup, pstep):

tests/test_stats.py

@@ -107,15 +107,15 @@ def test_residuals():
     np.testing.assert_allclose(res, np.array([-1.0, -2.0, 1.0, 0.0, 0.5]))
 
 
-def test_chisq_no_priors():
+def test_chisq():
     data   = np.array([1.1, 1.2, 0.9, 1.0])
     model  = np.array([1.0, 1.0, 1.0, 1.0])
     uncert = np.array([0.1, 0.1, 0.1, 0.1])
     chisq  = ms.chisq(model, data, uncert)
     assert chisq == 6.0
 
 
-def test_chisq():
+def test_chisq_priors():
     data   = np.array([1.1, 1.2, 0.9, 1.0])
     model  = np.array([1.0, 1.0, 1.0, 1.0])
     uncert = np.array([0.1, 0.1, 0.1, 0.1])
@@ -127,6 +127,33 @@ def test_chisq():
     assert chisq == 6.25
 
 
+def test_dwt_chisq():
+    data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0])
+    model = np.ones(8)
+    params = np.array([1.0, 0.1, 0.1])
+    chisq = ms.dwt_chisq(model, data, params)
+    np.testing.assert_allclose(chisq, 1693.22308882)
+
+
+def test_dwt_chisq_priors():
+    data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0])
+    model = np.ones(8)
+    params = np.array([1.0, 0.1, 0.1])
+    priors = np.array([1.0, 0.2, 0.3])
+    plow   = np.array([0.0, 0.0, 0.1])
+    pup    = np.array([0.0, 0.0, 0.1])
+    chisq = ms.dwt_chisq(model, data, params, priors, plow, pup)
+    np.testing.assert_allclose(chisq, 1697.2230888243134)
+
+
+def test_dwt_chisq_params_error():
+    data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0])
+    model = np.ones(8)
+    params = np.array([1.0, 0.1])
+    with pytest.raises(SystemExit):
+        chisq = ms.dwt_chisq(model, data, params)
+
+
 def test_log_prior_uniform():
     post = np.array([[3.0, 2.0], [3.1, 1.0], [3.6, 1.5]])
     prior    = np.array([3.5, 0.0])