Merge pull request #133 from pcubillos/fork

Fork

.travis.yml

@@ -1,15 +1,16 @@
 language: python
+os: linux
+dist: focal
 
-matrix:
+jobs:
   include:
   - python: 3.6
   - python: 3.7
-    sudo: required
-    dist: xenial
+  - python: 3.8
+  - python: 3.9
 
 install:
   - pip install -r requirements.txt
-  - pip install dynesty>=0.9.5
   - pip install -e .
 
 script: pytest tests -v

mc3/VERSION.py

@@ -4,6 +4,6 @@
 # mc3 Version:
 MC3_VER = 3  # Major version
 MC3_MIN = 0  # Minor version
-MC3_REV = 8  # Revision
+MC3_REV = 9  # Revision
 
 __version__ = f'{MC3_VER}.{MC3_MIN}.{MC3_REV}'

mc3/chain.py

@@ -15,7 +15,7 @@
 warnings.simplefilter("ignore", RuntimeWarning)
 
 
-class Chain(mp.Process):
+class Chain(mp.get_context('fork').Process):
   """
   Background process.  This guy evaluates the model and calculates chisq.
   """

mc3/mcmc_driver.py

@@ -160,7 +160,7 @@ def mcmc(data, uncert, func, params, indparams, pmin, pmax, pstep,
     sm_log_post = mpr.Array(ctypes.c_double, zlen)
     log_post = np.ctypeslib.as_array(sm_log_post.get_obj())
     # Chain index for given state in the Z array:
-    sm_zchain = mpr.Array(ctypes.c_int, -np.ones(zlen, np.int))
+    sm_zchain = mpr.Array(ctypes.c_int, -np.ones(zlen, int))
     zchain = np.ctypeslib.as_array(sm_zchain.get_obj())
     # Current number of samples in the Z array:
     zsize = mpr.Value(ctypes.c_int, M0)
@@ -197,10 +197,11 @@ def mcmc(data, uncert, func, params, indparams, pmin, pmax, pstep,
     ncpp[0:nchains % ncpu] += 1
 
     # Launch Chains:
+    mp_context = mpr.get_context('fork')
     pipes  = []
     chains = []
     for i in range(ncpu):
-        p = mpr.Pipe()
+        p = mp_context.Pipe()
         pipes.append(p[0])
         chains.append(
             ch.Chain(func, indparams, p[1], data, uncert,

mc3/ns_driver.py

@@ -25,7 +25,7 @@ def resample_equal(weights, rstate=None):
     nsamples = len(weights)
     positions = (rstate.random() + np.arange(nsamples)) / nsamples
 
-    idx = np.zeros(nsamples, dtype=np.int)
+    idx = np.zeros(nsamples, dtype=int)
     cumulative_sum = np.cumsum(weights)
     i, j = 0, 0
     while i < nsamples:

mc3/stats/stats.py

@@ -24,8 +24,8 @@
 from .. import utils as mu
 sys.path.append(mu.ROOT + 'mc3/lib/')
 import _binarray as ba
-import _chisq    as cs
-import _dwt      as dwt
+import _chisq as cs
+import _dwt as dwt
 
 
 def bin_array(data, binsize, uncert=None):
@@ -561,7 +561,7 @@ def dwt_daub4(array, inverse=False):
     Examples
     --------
     >>> import numpy as np
-    >>> improt matplotlib.pyplot as plt
+    >>> import matplotlib.pyplot as plt
     >>> import mc3.stats as ms
 
     >>> # Calculate the inverse DWT for a unit vector:

requirements.txt

@@ -1,4 +1,4 @@
 numpy>=1.15.0
 scipy>=0.17.1
 matplotlib>=2.0
-pytest>=3.9
+pytest>=6.0

setup.py

@@ -51,7 +51,7 @@
         'matplotlib>=2.0',
     ],
     tests_require = [
-        'pytest>=3.9',
+        'pytest>=6.0',
         'dynesty>=0.9.5',
     ],
     include_package_data=True,

src_c/_dwt.c

@@ -138,7 +138,7 @@ output will zero-padded to the next size of the form 2**M.             \n\
 Examples                                                               \n\
 --------                                                               \n\
 >>> import numpy as np                                                 \n\
->>> improt matplotlib.pyplot as plt                                    \n\
+>>> import matplotlib.pyplot as plt                                    \n\
 >>> import mc3.stats as ms                                             \n\
                                                                        \n\
 >>> # Calculate the inverse DWT for a unit vector:                     \n\
@@ -157,14 +157,15 @@ static PyObject *daub4(PyObject *self, PyObject *args){
       int isign,  /* Sign indicator to perform DWT or inverse DWT */
           asize, dwt_size,
           M, j;
+    npy_intp size[1];
 
       if (!PyArg_ParseTuple(args, "Oi", &array, &isign))
           return NULL;
 
       /* Expand to a size proportional to 2^M (zero padding) */
       asize = (int)PyArray_DIM(array, 0);
       M = ceil(1.0*log2(asize));
-      dwt_size = (int)pow(2, M);
+      size[0] = dwt_size = (int)pow(2, M);
 
       /* Allocate memory for pointer with the data: */
       ptr_array = malloc(dwt_size*sizeof(double));
@@ -176,7 +177,7 @@ static PyObject *daub4(PyObject *self, PyObject *args){
       /* Calculate the discrete wavelet transform: */
       dwt(ptr_array, dwt_size, isign);
 
-      dwt_array = (PyArrayObject *) PyArray_FromDims(1, &dwt_size, NPY_DOUBLE);
+      dwt_array = (PyArrayObject *) PyArray_SimpleNew(1, size, NPY_DOUBLE);
       for(j=0; j<dwt_size; j++)
           INDd(dwt_array,j) = ptr_array[j];
 

tests/test_dynesty.py

@@ -48,6 +48,7 @@ def quad(p, x):
 sampler = 'dynesty'
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_minimal():
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
         indparams=[x], pmin=pmin, pmax=pmax,
@@ -56,13 +57,15 @@ def test_dynesty_minimal():
     assert output is not None
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_ncpu():
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
         indparams=[x], pmin=pmin, pmax=pmax, ncpu=8,
         pstep=pstep, sampler=sampler, maxiter=5000)
     assert output is not None
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_shared():
     output = mc3.sample(data1, uncert1, func=quad, params=np.copy(params),
         sampler=sampler, indparams=[x], pmin=pmin, pmax=pmax,
@@ -71,6 +74,7 @@ def test_dynesty_shared():
     assert output['bestp'][1] == output['bestp'][0]
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_fixed():
     pars = np.copy(params)
     pars[0] = p0[0]
@@ -85,6 +89,7 @@ def test_dynesty_fixed():
     assert output['stdp'][0] == 0
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 @pytest.mark.parametrize('leastsq', ['lm', 'trf'])
 def test_dynesty_optimize(capsys, leastsq):
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
@@ -98,6 +103,7 @@ def test_dynesty_optimize(capsys, leastsq):
         np.array([4.28263253, -2.40781859, 0.49534411]), rtol=1e-7)
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_priors_gauss():
     prior    = np.array([ 4.5,  0.0,   0.0])
     priorlow = np.array([ 0.1,  0.0,   0.0])
@@ -111,6 +117,7 @@ def test_dynesty_priors_gauss():
     assert np.all(-2*output['log_post'] > output['chisq'])
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_log(capsys, tmp_path):
     os.chdir(str(tmp_path))
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
@@ -123,6 +130,7 @@ def test_dynesty_log(capsys, tmp_path):
     assert "NS.log" in os.listdir(".")
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_savefile(capsys, tmp_path):
     os.chdir(str(tmp_path))
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
@@ -137,6 +145,7 @@ def test_dynesty_savefile(capsys, tmp_path):
 
 
 # Trigger errors:
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_pmin_error(capsys):
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
         sampler=sampler, indparams=[x], pstep=pstep, pmax=pmax)
@@ -146,6 +155,7 @@ def test_dynesty_pmin_error(capsys):
         in captured.out
 
 
+@pytest.mark.skip(reason="Stop dynesty support")
 def test_dynesty_pmax_error(capsys):
     output = mc3.sample(data, uncert, func=quad, params=np.copy(params),
         sampler=sampler, indparams=[x], pstep=pstep, pmin=pmin)

tests/test_logging.py

@@ -74,5 +74,5 @@ def test_tracklev(tmp_path):
 
     with open(log_file, 'r') as f:
         content = f.read()
-    assert "Error in module: '__init__.py', function: 'main', line:" in content
+    assert "Error in module: '__init__.py', function: 'console_main'" in content
 

tests/test_stats.py

@@ -67,6 +67,24 @@
      0.73813656, 0.70672705, 0.67847008, 0.65371869, 0.63189428,
      0.6125697 , 0.59543317, 0.58026767, 0.56694288, 0.55198132])
 
+expected_daub4_inverse = np.array([
+    -0.0301851821, -0.0522822690, -0.0662912607, -0.0824674511, -0.0905555462,
+    -0.1008108399, -0.1132333322, -0.1250751254,  0.1325825215,  0.3180280110,
+     0.4312613433,  0.5638438647,  0.1412513157, -0.1325825215, -0.2576576469,
+    -0.4225925490, -0.1671021007, -0.0242642855,  0.0059208966,  0.0662912607,
+     0.0140089918, -0.0080880952,  0.0000000000,  0.0000000000,  0.0000000000,
+     0.0000000000,  0.0000000000,  0.0000000000,  0.0000000000,  0.0000000000,
+     0.0000000000,  0.0000000000,])
+
+expected_daub4_forward = np.array([
+     0.1625300592, 0.0874699408, -0.0463140877,  0.2795672632, -0.0905555462,
+     0.0000000000, 0.0140089918,  0.1412513157,  0.3537658774, -0.0625000000,
+     0.0000000000, 0.0000000000,  0.0000000000,  0.0000000000,  0.0000000000,
+    -0.1082531755, 0.0000000000,  0.8365163037, -0.1294095226,  0.0000000000,
+     0.0000000000, 0.0000000000,  0.0000000000,  0.0000000000,  0.0000000000,
+     0.0000000000, 0.0000000000,  0.0000000000,  0.0000000000,  0.0000000000,
+     0.0000000000, 0.0000000000,])
+
 
 def test_bin_array_unweighted():
     data = np.array([0,1,2, 3,3,3, 3,3,4])
@@ -193,7 +211,7 @@ def test_log_prior_single_sample():
     log_prior = ms.log_prior(params, prior, priorlow, priorup, pstep)
     np.testing.assert_allclose(log_prior, -12.5)
 
-    
+
 
 def test_cred_region():
     np.random.seed(2)
@@ -220,8 +238,8 @@ def test_ppf_uniform_array():
         np.array([-10.0, 0.0, 10.]))
 
 
-@pytest.mark.parametrize('u, result', 
-    [(1e-10, -6.361340902404056), 
+@pytest.mark.parametrize('u, result',
+    [(1e-10, -6.361340902404056),
      (0.5,   0.0),
      (1.0-1e-10, 6.361340889697422)])
 def test_ppf_gaussian_scalar(u, result):
@@ -244,6 +262,32 @@ def test_ppf_gaussian_two_sided():
     np.testing.assert_allclose(np.array(ppf_g(u)), result)
 
 
+def test_dwt_daub4_inverse():
+    nx = 32
+    e4 = np.zeros(nx)
+    e4[4] = 1.0
+    ie4 = ms.dwt_daub4(e4, True)
+    np.testing.assert_allclose(ie4, expected_daub4_inverse)
+
+
+def test_dwt_daub4_forward():
+    nx = 32
+    ie4 = np.zeros(nx)
+    ie4[4] = 1.0
+    fe4 = ms.dwt_daub4(ie4)
+    np.testing.assert_allclose(fe4, expected_daub4_forward)
+
+
+def test_dwt_daub4_inverse_forward():
+    nx = 32
+    e4 = np.zeros(nx)
+    e4[4] = 1.0
+    ie4 = ms.dwt_daub4(e4, True)
+    fe4 = ms.dwt_daub4(ie4)
+    np.testing.assert_allclose(fe4, e4, atol=1e-8)
+
+
+
 def test_timeavg_values_red():
     rms, rmslo, rmshi, stderr, binsz = ms.time_avg(data, len(data)/10, 5)
     np.testing.assert_almost_equal(rms,    expected_red_rms)