add option to exclude hip iad from likelihood

sblunt · May 2, 2024 · c601501 · c601501
1 parent 0f6ae0f
commit c601501
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 21 deletions.
diff --git a/orbitize/hipparcos.py b/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:

diff --git a/orbitize/sampler.py b/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(