Merge branch 'v3' into rebound-bugfix

docs/requirements.txt

@@ -1,3 +1,4 @@
 nbsphinx
 pandoc
-ipykernel
+ipykernel
+sphinx_rtd_theme

tests/test_OFTI.py

@@ -126,7 +126,7 @@ def test_run_sampler():
 
     # test to make sure outputs are reasonable
     start = time.time()
-    orbits = s.run_sampler(1000, num_cores=4)
+    orbits = s.run_sampler(100, num_cores=4)
     end = time.time()
 
     print()
@@ -155,15 +155,15 @@ def test_run_sampler():
     ecc_exp = 0.19
     inc_exp = np.radians(140)
 
-    # test to make sure OFTI values are within 20% of expectations
-    assert sma == pytest.approx(sma_exp, abs=0.2 * sma_exp)
-    assert ecc == pytest.approx(ecc_exp, abs=0.2 * ecc_exp)
-    assert inc == pytest.approx(inc_exp, abs=0.2 * inc_exp)
+    # test to make sure OFTI values are consistent with expectations
+    assert sma == pytest.approx(sma_exp, abs=5)
+    assert ecc == pytest.approx(ecc_exp, abs=0.1)
+    assert inc == pytest.approx(inc_exp, abs=np.radians(10))
 
     sma_seppa = sma  # use for covarinaces test
 
     # test with only one core
-    orbits = s.run_sampler(100, num_cores=1)
+    orbits = s.run_sampler(10, num_cores=1)
 
     # test with only one epoch
     myDriver = orbitize.driver.Driver(
@@ -410,12 +410,13 @@ def test_OFTI_pan_priors():
 
 
 if __name__ == "__main__":
-    test_scale_and_rotate()
-    test_run_sampler()
-    test_OFTI_covariances()
-    test_OFTI_multiplanet()
-    test_not_implemented()
-    test_fixed_sys_params_sampling()
+    # test_scale_and_rotate()
+    # test_run_sampler()
+
+    # test_OFTI_covariances()
+    # test_OFTI_multiplanet()
+    # test_not_implemented()
+    # test_fixed_sys_params_sampling()
     test_OFTI_pan_priors()
-    # profile_system()
-    print("Done!")
+    # # profile_system()
+    # print("Done!")

tests/test_chi2_lnlike.py

@@ -1,8 +1,9 @@
-import orbitize.driver 
+import orbitize.driver
 import numpy as np
 import orbitize.lnlike as lnlike
 import pytest
 
+
 def test_chi2lnlike():
     """
     Test the ability of ``orbitize.lnlike.chi2_lnlike()``
@@ -19,7 +20,7 @@ def test_chi2lnlike():
     chi2 = lnlike.chi2_lnlike(data, errors, None, model, jitter, seppa_indices)
     assert chi2.shape == (3, 2)
     assert chi2 == pytest.approx(
-        -0.5 * np.ones((3, 2)) - np.log(np.sqrt(2*np.pi*np.ones((3, 2))))
+        -0.5 * np.ones((3, 2)) - np.log(np.sqrt(2 * np.pi * np.ones((3, 2))))
     )
 
     # test with multiple models
@@ -33,7 +34,7 @@ def test_chi2lnlike():
     chi2 = lnlike.chi2_lnlike(data, errors, None, model, jitter, seppa_indices)
     assert chi2.shape == (3, 2, 5)
     assert chi2 == pytest.approx(
-        -0.5 * np.ones((3, 2, 5)) - np.log(np.sqrt(2*np.pi*np.ones((3, 2, 5))))
+        -0.5 * np.ones((3, 2, 5)) - np.log(np.sqrt(2 * np.pi * np.ones((3, 2, 5))))
     )
 
 
@@ -42,29 +43,29 @@ def test_chi2lnlike_withcov():
     Test the ability of ``orbitize.lnlike.chi2_lnlike()`` to work with some or all data having covariances
     """
     ### all all covariances
-    data = np.array([[5,-4], [3,-2], [1,0] ])
+    data = np.array([[5, -4], [3, -2], [1, 0]])
     model = np.zeros(data.shape)
     jitter = np.zeros(data.shape)
-    errs = np.array([[2,2], [2,2], [2,2]])
+    errs = np.array([[2, 2], [2, 2], [2, 2]])
     covs = np.array([1, 0.25, 0.25])
-    corrs = covs/errs[:,0]/errs[:,1]
+    corrs = covs / errs[:, 0] / errs[:, 1]
 
     chi2s = lnlike.chi2_lnlike(data, errs, corrs, model, jitter, [])
 
     residuals = data - model
     for res, err, cov, chi2 in zip(residuals, errs, covs, chi2s):
-        cov_matrix = np.array([[err[0]**2, cov], [cov, err[1]**2]])
+        cov_matrix = np.array([[err[0] ** 2, cov], [cov, err[1] ** 2]])
         cov_inv = np.linalg.inv(cov_matrix)
         cov_inv_dot_diff = np.dot(cov_inv, res)
         logdet = np.linalg.slogdet(cov_matrix)[1]
         res_cov_res = res.dot(cov_inv_dot_diff)
-        numpy_chi2 = -0.5 * (res_cov_res + logdet + 2 * np.log(2 * np.pi)) 
+        numpy_chi2 = -0.5 * (res_cov_res + logdet + 2 * np.log(2 * np.pi))
 
         assert np.sum(chi2) == pytest.approx(numpy_chi2)
 
     ### only one covariance term
     covs = np.array([1, np.nan, np.nan])
-    corrs = covs/errs[:,0]/errs[:,1]
+    corrs = covs / errs[:, 0] / errs[:, 1]
     new_chi2s = lnlike.chi2_lnlike(data, errs, corrs, model, jitter, [])
 
     assert chi2s[0] == pytest.approx(new_chi2s[0])
@@ -74,91 +75,114 @@ def test_2x2_analytical_solution():
     """
     Tests that our analytical solution to the 2x2 covariance matrix is correct
     """
-    residuals = np.array([[5,-4], [3,-2], [1,0] ])
+    residuals = np.array([[5, -4], [3, -2], [1, 0]])
 
-    errs = np.array([[2,2], [2,2], [2,2]])
+    errs = np.array([[2, 2], [2, 2], [2, 2]])
     covs = np.array([1, 0.25, 0.25])
-    corrs = covs/errs[:,0]/errs[:,1]
+    corrs = covs / errs[:, 0] / errs[:, 1]
 
     chi2s = lnlike._chi2_2x2cov(np.array([residuals]), np.array([errs**2]), corrs)
 
     # compare to numpy solution
     for res, err, cov, chi2 in zip(residuals, errs, covs, chi2s[0]):
-        cov_matrix = np.array([[err[0]**2, cov], [cov, err[1]**2]])
+        cov_matrix = np.array([[err[0] ** 2, cov], [cov, err[1] ** 2]])
         cov_inv = np.linalg.inv(cov_matrix)
         cov_inv_dot_diff = np.dot(cov_inv, res)
         logdet = np.linalg.slogdet(cov_matrix)[1]
         res_cov_res = res.dot(cov_inv_dot_diff)
-        numpy_chi2 = -0.5 * (res_cov_res + logdet + 2 * np.log(2 * np.pi)) 
+        numpy_chi2 = -0.5 * (res_cov_res + logdet + 2 * np.log(2 * np.pi))
 
         assert np.sum(chi2) == pytest.approx(numpy_chi2)
 
 
 def test_chi2_log():
-    #initiate OFTI driver with chi2 log
+    # initiate OFTI driver with chi2 log
     myDriver = orbitize.driver.Driver(
-    '{}/GJ504.csv'.format(orbitize.DATADIR), 'OFTI', 1, 1.22, 56.95, mass_err=0.08, plx_err=0.26, chi2_type='log')
+        "{}/GJ504.csv".format(orbitize.DATADIR),
+        "OFTI",
+        1,
+        1.22,
+        56.95,
+        mass_err=0.08,
+        plx_err=0.26,
+        chi2_type="log",
+    )
     s = myDriver.sampler
-    params = [44, 0, 45*np.pi/180, 0, 325*np.pi/180, 0, 56.95, 1.22]
+    params = [44, 0, 45 * np.pi / 180, 0, 325 * np.pi / 180, 0, 56.95, 1.22]
     log_chi2 = s._logl(params)
 
     sys = myDriver.system
-    data = np.array([sys.data_table['quant1'], sys.data_table['quant2']]).T
-    errors = np.array([sys.data_table['quant1_err'], sys.data_table['quant2_err']]).T
+    data = np.array([sys.data_table["quant1"], sys.data_table["quant2"]]).T
+    errors = np.array([sys.data_table["quant1_err"], sys.data_table["quant2_err"]]).T
     model, jitter = sys.compute_model(params)
 
-    sep_data = data[:,0]
+    sep_data = data[:, 0]
     sep_model = model[:, 0]
-    sep_error = errors[:,0]
-    pa_data = data[:,1]
+    sep_error = errors[:, 0]
+    pa_data = data[:, 1]
     pa_model = model[:, 1]
-    pa_error = errors[:,1]*np.pi/180
+    pa_error = errors[:, 1] * np.pi / 180
 
+    # calculating sep chi squared
+    sep_chi2_log = (np.log(sep_data) - np.log(sep_model)) ** 2 / (
+        sep_error / sep_data
+    ) ** 2
 
-    #calculating sep chi squared
-    sep_chi2_log = (np.log(sep_data)-np.log(sep_model))**2/(sep_error/sep_data)**2
+    # calculting pa chi squared Log
+    pa_resid = (pa_model - pa_data + 180.0) % 360.0 - 180.0
+    pa_chi2_log = 2 * (1 - np.cos(pa_resid * np.pi / 180)) / pa_error**2
 
-    #calculting pa chi squared Log
-    pa_resid = (pa_model-pa_data +180.) % 360. - 180.
-    pa_chi2_log = 2*(1-np.cos(pa_resid*np.pi/180))/pa_error**2
+    chi2 = np.zeros((len(sep_data), 2))
 
-    chi2 = np.zeros((len(sep_data),2))
+    sigma2 = errors**2 + jitter**2
 
-    sigma2 = errors**2 + jitter**2 
+    chi2[:, 0] = sep_chi2_log
+    chi2[:, 1] = pa_chi2_log
 
-    chi2[:,0] = sep_chi2_log
-    chi2[:,1] = pa_chi2_log
-
-    chi2 = -0.5 * chi2 - np.log(np.sqrt(2*np.pi*sigma2))
+    chi2 = -0.5 * chi2 - np.log(np.sqrt(2 * np.pi * sigma2))
 
     lnlike = np.sum(chi2)
 
     assert lnlike == pytest.approx(log_chi2)
 
 
 def test_log_vs_standard():
-
-    #initiate driver with standard chi2
-    myDriver_standard = orbitize.driver.Driver('{}/GJ504.csv'.format(orbitize.DATADIR),'OFTI', 1, 1.22, 56.95, mass_err=0.08, plx_err=0.26)
+    # initiate driver with standard chi2
+    myDriver_standard = orbitize.driver.Driver(
+        "{}/GJ504.csv".format(orbitize.DATADIR),
+        "OFTI",
+        1,
+        1.22,
+        56.95,
+    )
     s_standard = myDriver_standard.sampler
-    orbits = s_standard.run_sampler(3000)
-
-    #initiate driver with log chi2
-    myDriver_log = orbitize.driver.Driver('{}/GJ504.csv'.format(orbitize.DATADIR), 'OFTI', 1, 1.22, 56.95, mass_err=0.08, plx_err=0.26, chi2_type = 'log')
+    _ = s_standard.run_sampler(500)
+
+    # initiate driver with log chi2
+    myDriver_log = orbitize.driver.Driver(
+        "{}/GJ504.csv".format(orbitize.DATADIR),
+        "OFTI",
+        1,
+        1.22,
+        56.95,
+        chi2_type="log",
+    )
     s_log = myDriver_log.sampler
-    orbits = s_log.run_sampler(10000)   
+    _ = s_log.run_sampler(500)
 
-    #take mean of result objects
-    myResults_standard = np.mean(s_standard.results.post,axis=0)
-    myResults_log = np.mean(s_log.results.post,axis=0)
+    # take mean of result objects
+    myResults_standard = np.mean(s_standard.results.post, axis=0)
+    myResults_log = np.mean(s_log.results.post, axis=0)
 
-    assert myResults_log == pytest.approx(myResults_standard, rel=0.05) 
+    print(myResults_log[1], myResults_standard[1])
 
+    # check that the eccentricity means are about the same
+    assert myResults_log[1] == pytest.approx(myResults_standard[1], rel=0.1)
 
-if __name__ == "__main__":
-    test_chi2lnlike()
-    test_chi2_log()
-    test_chi2lnlike_withcov()
-    test_chi2lnlike_withcov()
-    test_2x2_analytical_solution()
 
+if __name__ == "__main__":
+    # test_chi2lnlike()
+    # test_chi2_log()
+    # test_chi2lnlike_withcov()
+    # test_2x2_analytical_solution()
+    test_log_vs_standard()