add tests of pykima.keplerian and of default prior example

pykima/display.py

@@ -46,10 +46,6 @@ def percentile68_ranges_latex(a):
     lp, median, up = np.percentile(a, [16, 50, 84])
     return r'$%.2f ^{+%.2f} _{-%.2f}$' % (median, up-median, median-lp)
 
-# def get_aliases(Preal):
-#     fs = np.array([0.0027381631, 1.0, 1.0027]) #, 0.018472000025212765])
-#     return np.array([abs(1 / (1./Preal + i*fs)) for i in range(-2, 2)]).T
-
 
 def get_planet_mass(P, K, e, star_mass=1.0, full_output=False, verbose=False):
     if verbose: print('Using star mass = %s solar mass' % star_mass)
@@ -97,7 +93,7 @@ def get_skip(line):
             load_args = re.findall(r'\((.*?)\)', line, re.DOTALL)[1]
             load_args = load_args.split(',')
             if len(load_args) == 3:
-                # use gave 'skip' option
+                # user gave 'skip' option
                 return int(load_args[2])
             else:
                 # default is skip=2
@@ -179,7 +175,6 @@ def get_skip(line):
             n_trend = 0
 
 
-
         # find fiber offset in the compiled model
         if fiber_offset is None:
             try:

pykima/keplerian.py

@@ -1,3 +1,5 @@
+__all__ = ['keplerian', 'true_anomaly', 'ecc_anomaly']
+
 import numpy as np 
 pi = np.pi
 

tests/conftest.py

@@ -0,0 +1,13 @@
+import pytest
+def pytest_addoption(parser):
+    parser.addoption("--slow", action="store_true",
+                     default=False, help="run slow tests")
+
+def pytest_collection_modifyitems(config, items):
+    if config.getoption("--slow"):
+        # --runslow given in cli: do not skip slow tests
+        return
+    skip_slow = pytest.mark.skip(reason="need --slow option to run")
+    for item in items:
+        if "slow" in item.keywords:
+            item.add_marker(skip_slow)

tests/test_kima.py

@@ -0,0 +1,68 @@
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+def test_import():
+    import pykima
+
+def test_true_ecc_anomaly():
+    from pykima.keplerian import true_anomaly, ecc_anomaly
+
+    # true anomaly f
+    npt.assert_allclose(true_anomaly(0., 0.), 0.0)
+    npt.assert_allclose(true_anomaly(np.pi, 0.), np.pi)
+
+    # eccentric anomaly E
+    npt.assert_allclose(ecc_anomaly(0., 0.), 0.)
+    npt.assert_allclose(ecc_anomaly(np.pi, 0.), np.pi)
+    npt.assert_allclose(ecc_anomaly(2*np.pi, 0.), 2*np.pi)
+    npt.assert_allclose(ecc_anomaly(0., 0.2), 0.)
+    # checked with http://orbitsimulator.com/sheela/kepler.htm
+    E = np.rad2deg(ecc_anomaly(np.deg2rad(80), 0.2))
+    npt.assert_allclose(E, 91.45545886486815, rtol=1e-8)
+
+def test_keplerian():
+    from pykima.keplerian import keplerian
+    # keplerian(time, p, k, ecc, omega, t0, vsys)
+
+    times1 = [0.,]
+    times2 = [0., 1., 2.]
+
+    with pytest.warns(RuntimeWarning):
+        result = keplerian(times1, 0., 0., 0., 0., 0., 0.)
+        assert np.isnan(result)
+
+        result = keplerian(times2, 0., 0., 0., 0., 0., 0.)
+        assert np.all(np.isnan(result))
+
+    npt.assert_allclose(keplerian(times1, 100., 1., 0., 0., 0., 0.), 1.)
+    npt.assert_allclose(keplerian(times2, 1., 1., 0., 0., 0., 0.), 1.)
+
+
+
+def test_percentiles():
+    from pykima.display import percentile68_ranges, percentile68_ranges_latex
+
+    a = np.random.uniform(size=int(1e6))
+    b = np.random.randn(int(1e6))
+
+    ra = percentile68_ranges(a)
+    rb = percentile68_ranges(b)
+
+    npt.assert_allclose(ra[0], np.median(a))
+    npt.assert_allclose(ra[0]+ra[1], np.percentile(a, 84))
+    npt.assert_allclose(ra[0]-ra[2], np.percentile(a, 16))
+    npt.assert_allclose(ra[0], 0.5, rtol=1e-2)
+
+    npt.assert_allclose(rb[0], np.median(b))
+    npt.assert_allclose(rb[1], 0.9944578832097531677397, rtol=1e-2)
+    npt.assert_allclose(rb[2], 0.9944578832097531677397, rtol=1e-2)
+    npt.assert_allclose(rb[0], 0.0, atol=1e-2)
+
+    a = np.linspace(0, 10, 100)
+    npt.assert_allclose(percentile68_ranges(a), [5.0, 3.4, 3.4])
+    assert percentile68_ranges_latex(a) == '$5.00 ^{+3.40} _{-3.40}$'
+
+
+def test_KimaResults():
+    pass

tests/test_priors.py

@@ -0,0 +1,54 @@
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+import os
+
+def example_priors_run():
+	""" 
+	Run kima with maxlevels=1 to sample from the default priors
+	"""
+	os.system('make -S' + ' >/dev/null 2>&1') # compile
+	os.system('./run' + ' >/dev/null 2>&1') # run
+
+	assert os.path.exists('sample.txt')
+	assert os.path.exists('sample_info.txt')
+
+
+def cleanup():
+	os.system('make -S cleanout' + ' >/dev/null 2>&1')
+	assert not os.path.exists('sample.txt')
+
+
+@pytest.mark.slow
+def test_priors():
+	""" 
+	After running kima with maxlevels=1 to sample from the default priors,
+	test some basic statistics of the samples
+	"""
+	topdir = os.getcwd()
+	os.chdir('examples/default_priors')
+
+	example_priors_run()
+
+	extra_sigma, vsys, P, ecc = np.loadtxt('sample.txt', unpack=True,
+		                                   usecols=(0, -1, 4, 7))
+
+	npt.assert_allclose(extra_sigma.min(), 0., rtol=0, atol=1e-1)
+	npt.assert_allclose(extra_sigma.max(), 99., rtol=1e-1, atol=0)
+
+	npt.assert_allclose(vsys.min(), -1000., rtol=0, atol=1)
+	npt.assert_allclose(vsys.max(), 1000., rtol=0, atol=1)	
+
+	npt.assert_allclose(P.min(), 1., rtol=1e-2, atol=0)
+	npt.assert_allclose(P.max(), 1e5, rtol=1e-2, atol=0)
+	npt.assert_allclose(P.mean(), (1e5 - 1.)/(np.log(1e5) - np.log(1.)),
+	                    rtol=1, atol=0)
+
+	npt.assert_allclose(ecc.min(), 0., rtol=0, atol=1e-4)
+	npt.assert_allclose(ecc.max(), 1., rtol=0, atol=1e-4)
+	npt.assert_allclose(ecc.mean(), 0.5, rtol=0, atol=1e-1)
+
+
+	cleanup()
+	os.chdir(topdir)