Merge pull request #84 from sblunt/mcmc-bugfixes

Mcmc bugfixes
sblunt · Dec 4, 2018 · e291cd8 · e291cd8
2 parents 68567f3 + 38552fa
commit e291cd8
Show file tree

Hide file tree

Showing 10 changed files with 175 additions and 93 deletions.
diff --git a/docs/tutorials/MCMC_tutorial.ipynb b/docs/tutorials/MCMC_tutorial.ipynb
diff --git a/docs/tutorials/MCMC_vs_OFTI.ipynb b/docs/tutorials/MCMC_vs_OFTI.ipynb
diff --git a/docs/tutorials/Modifying_Priors.ipynb b/docs/tutorials/Modifying_Priors.ipynb
diff --git a/docs/tutorials/OFTI_tutorial.ipynb b/docs/tutorials/OFTI_tutorial.ipynb
diff --git a/docs/tutorials/Plotting_tutorial.ipynb b/docs/tutorials/Plotting_tutorial.ipynb
diff --git a/orbitize/lnlike.py b/orbitize/lnlike.py
@@ -19,7 +19,7 @@ def chi2_lnlike(data, errors, model, seppa_indices):
             given in sep/PA. This list is located in System.seppa.
 
     Returns:
-    	np.array: Nobsx2xM array of chi-squared values. 
+    	np.array: Nobsx2xM array of chi-squared values.  
     
     .. note::
 
@@ -40,9 +40,7 @@ def chi2_lnlike(data, errors, model, seppa_indices):
 
     # if there are PA values, we should take the difference modulo angle wrapping
     if np.size(seppa_indices) > 0:
-        residual[:, seppa_indices, 1] = np.arctan2(
-            np.sin(data[seppa_indices, 1] - model[:, seppa_indices, 1]), 
-            np.cos(data[seppa_indices, 1] - model[:, seppa_indices, 1]))
+        residual[:, seppa_indices, 1] = (residual[:, seppa_indices, 1] + 180.) % 360. - 180.
 
     chi2 = -0.5 * residual**2 / errors**2
 

diff --git a/orbitize/results.py b/orbitize/results.py
@@ -9,6 +9,7 @@
 import orbitize.kepler as kepler
 import h5py
 from astropy.io import fits
+from astropy.time import Time
 
 # define modified color map for default use in orbit plots
 cmap = mpl.cm.Purples_r
@@ -228,25 +229,25 @@ def plot_corner(self, param_list=[], **corner_kwargs):
         default_labels = [
             'a [au]',
             'ecc',
-            'inc [deg]',
-            '$\omega$ [deg]',
-            '$\Omega$ [deg]',
+            'inc [rad]',
+            '$\omega$ [rad]',
+            '$\Omega$ [rad]',
             '$\\tau$',
-            '$M_T$ [Msol]',
-            '$\pi$ [mas]'
+            '$\pi$ [mas]',
+            '$M_T$ [Msol]'
         ]
         if len(param_list)>0: # user chose to plot specific parameters only
             num_orb_param = self.post.shape[1] # number of orbital parameters (+ mass, parallax)
             num_objects,remainder = np.divmod(num_orb_param,6)
             have_mtot_and_plx = remainder == 2
             param_indices = []
             for param in param_list:
-                if param=='mtot':
+                if param=='plx':
                     if have_mtot_and_plx:
                         param_indices.append(num_orb_param-2) # the 2nd last index
-                elif param=='plx':
+                elif param=='mtot':
                     if have_mtot_and_plx:
-                        param_indices.append(num_orb_param-2) # the last index
+                        param_indices.append(num_orb_param-1) # the last index
                 elif len(param)==4: # to prevent invalid, short param names breaking
                     if param[0:3] in dict_of_indices:
                         object_id = np.int(param[3])
@@ -359,7 +360,10 @@ def plot_orbits(self, parallax=None, total_mass=None, object_mass=0,
 
         # Create a linearly increasing colormap for our range of epochs
         norm = mpl.colors.Normalize(vmin=np.min(epochs), vmax=np.max(epochs[-1,:]))
-        norm_yr = mpl.colors.Normalize(vmin=np.min(epochs/365.25), vmax=np.max(epochs[-1,:]/365.25))
+        norm_yr = mpl.colors.Normalize(
+            vmin=np.min(Time(epochs,format='mjd').decimalyear), 
+            vmax=np.max(Time(epochs,format='mjd').decimalyear)
+        )
 
         # Create figure for orbit plots
         fig, ax = plt.subplots()

diff --git a/orbitize/sampler.py b/orbitize/sampler.py
@@ -376,7 +376,7 @@ def _logl(self, params, include_logp=False):
                 we need model predictions for. Must be in the same order
                 documented in System() above. If M=1, this can be a 1d array.
 
-            include_logp (bool): if True, also includ elog prior in this function
+            include_logp (bool): if True, also include log prior in this function
 
         Returns:
             lnlikes (float): sum of all log likelihoods of the data given input model

diff --git a/tests/GJ504.csv b/tests/GJ504.csv
@@ -1,10 +1,10 @@
 # Table 12 of Blunt et al 2017
 # Previous Observations of GJ 504 b
 epoch,object,sep,sep_err,pa,pa_err
-734601.7575,1,2479,16,327.94,0.39
-734660.1975,1,2483,8,327.45,0.19
-734740.5525,1,2481,33,326.84,0.94
-734744.205,1,2448,24,325.82,0.66
-734941.44,1,2483,15,326.46,0.36
-734985.27,1,2487,8,326.54,0.18
-735025.4475,1,2499,26,326.14,0.61
+55645.95,1,2479,16,327.94,0.39
+55702.89,1,2483,8,327.45,0.19
+55785.015,1,2481,33,326.84,0.94
+55787.935,1,2448,24,325.82,0.66
+55985.19400184,1,2483,15,326.46,0.36
+56029.11400323,1,2487,8,326.54,0.18
+56072.30200459,1,2499,26,326.14,0.61
diff --git a/tests/GJ504_1epoch.csv b/tests/GJ504_1epoch.csv
@@ -1,4 +1,4 @@
 # Table 12 of Blunt et al 2017
 # One-epoch Observation of GJ 504 b
 epoch,object,sep,sep_err,pa,pa_err
-734601.7575,1,2479,16,327.94,0.39
+55645.95,1,2479,16,327.94,0.39