Option to exclude hip iad from absastromlike & minor bugfix

orbitize/hipparcos.py

@@ -164,6 +164,8 @@ class HipparcosLogProb(object):
             should be False, but it's helpful for testing. Check out
             `orbitize.hipparcos.nielsen_iad_refitting_test()` for an example
             using this renormalization.
+        include_hip_iad_in_likelihood (bool): if False, then don't add the Hipparcos
+            log(likelihood) to the overall log(likelihood computed in sampler.py)
 
     Written: Sarah Blunt & Rob de Rosa, 2021
     """
@@ -175,9 +177,11 @@ def __init__(
         num_secondary_bodies,
         alphadec0_epoch=1991.25,
         renormalize_errors=False,
+        include_hip_iad_in_likelihood=True
     ):
         self.path_to_iad_file = path_to_iad_file
         self.renormalize_errors = renormalize_errors
+        self.include_hip_iad_in_likelihood = include_hip_iad_in_likelihood
 
         # infer if the IAD file is an older DVD file or a new file
         with open(path_to_iad_file, "r") as f:
@@ -282,14 +286,14 @@ def __init__(
         else:
             iad = np.transpose(np.loadtxt(path_to_iad_file))
 
-        n_lines = len(iad)
-
         times = iad[1] + 1991.25
         self.cos_phi = iad[3]  # scan direction
         self.sin_phi = iad[4]
         self.R = iad[5]  # abscissa residual [mas]
         self.eps = iad[6]  # error on abscissa residual [mas]
 
+        n_lines = len(self.eps)
+
         # reject negative errors (scans that were rejected by Hipparcos team)
         good_scans = np.where(self.eps > 0)[0]
 

orbitize/sampler.py

@@ -100,30 +100,31 @@ def _logl(self, params):
             lnlikes_sum += self.custom_lnlike(params)
 
         if self.system.hipparcos_IAD is not None:
-            # compute Ra/Dec predictions at the Hipparcos IAD epochs
-            raoff_model, deoff_model, _ = self.system.compute_all_orbits(
-                params, epochs=self.system.hipparcos_IAD.epochs_mjd
-            )
+            if self.system.hipparcos_IAD.include_hip_iad_in_likelihood:
+                # compute Ra/Dec predictions at the Hipparcos IAD epochs
+                raoff_model, deoff_model, _ = self.system.compute_all_orbits(
+                    params, epochs=self.system.hipparcos_IAD.epochs_mjd
+                )
 
-            (
-                raoff_model_hip_epoch,
-                deoff_model_hip_epoch,
-                _,
-            ) = self.system.compute_all_orbits(
-                params, epochs=Time([1991.25], format="decimalyear").mjd
-            )
+                (
+                    raoff_model_hip_epoch,
+                    deoff_model_hip_epoch,
+                    _,
+                ) = self.system.compute_all_orbits(
+                    params, epochs=Time([1991.25], format="decimalyear").mjd
+                )
 
-            # subtract off position of star at reference Hipparcos epoch
-            raoff_model[:, 0, :] -= raoff_model_hip_epoch[:, 0, :]
-            deoff_model[:, 0, :] -= deoff_model_hip_epoch[:, 0, :]
+                # subtract off position of star at reference Hipparcos epoch
+                raoff_model[:, 0, :] -= raoff_model_hip_epoch[:, 0, :]
+                deoff_model[:, 0, :] -= deoff_model_hip_epoch[:, 0, :]
 
-            # select body 0 raoff/deoff predictions & feed into Hip IAD lnlike fn
-            lnlikes_sum += self.system.hipparcos_IAD.compute_lnlike(
-                raoff_model[:, 0, :],
-                deoff_model[:, 0, :],
-                params,
-                self.system.param_idx,
-            )
+                # select body 0 raoff/deoff predictions & feed into Hip IAD lnlike fn
+                lnlikes_sum += self.system.hipparcos_IAD.compute_lnlike(
+                    raoff_model[:, 0, :],
+                    deoff_model[:, 0, :],
+                    params,
+                    self.system.param_idx,
+                )
 
         if self.system.gaia is not None:
             gaiahip_epochs = Time(

orbitize/system.py

@@ -454,8 +454,11 @@ def compute_all_orbits(self, params_arr, epochs=None, comp_rebound=False):
                     within_orbit = np.where(all_smas <= sma)
                     outside_orbit = np.where(all_smas > sma)
                     all_pl_masses = params_arr[self.secondary_mass_indx]
-                    inside_masses = all_pl_masses[within_orbit]
-                    mtot = np.sum(inside_masses) + m0
+                    inside_masses = all_pl_masses[within_orbit].reshape((-1, n_orbits))
+                    if n_orbits == 1:
+                        mtot = np.sum(inside_masses) + m0
+                    else:
+                        mtot = np.sum(inside_masses, axis=0) + m0
 
                 else:
                     m_pl = np.zeros(self.num_secondary_bodies)