add upperlimit functionality to all of the builtin models

docs/tutorials/AT2022cmc_Andreoni2022.txt

@@ -1,18 +1,18 @@
 facility,date,dt,nu,F_nu,F_err,upperlimit
-NOEMA,2022-03-24 22:14,41.48,86,3914,30,false
+NOEMA,2022-03-24 22:14,41.48,86,3914,30,true
 NOEMA,2022-03-24 22:14,41.48,102,3609,34,false
 NOEMA,2022-03-25 00:45,41.58,136,3045,41,false
 NOEMA,2022-03-25 00:45,41.58,152,2750,51,false
 VLA,2022-03-31 04:08,47.73,31.4,2130,30,false
 VLA,2022-03-31 04:08,47.73,33.5,2260,30,false
 VLA,2022-03-31 04:08,47.73,35.6,2350,40,false
-VLA,2022-03-31 04:08,47.73,37.5,2360,40,false
+VLA,2022-03-31 04:08,47.73,37.5,2360,40,true
 VLA,2022-03-31 04:13,47.73,8.5,270,12,false
 VLA,2022-03-31 04:13,47.73,9.5,336,11,false
 VLA,2022-03-31 04:13,47.73,10.5,385,12,false
 VLA,2022-03-31 04:13,47.73,11.5,438,14,false
 VLA,2022-03-31 04:23,47.74,12.8,583,12,false
 VLA,2022-03-31 04:23,47.74,14.3,724,12,false
-VLA,2022-03-31 04:23,47.74,15.9,801,14,false
+VLA,2022-03-31 04:23,47.74,15.9,801,14,true
 VLA,2022-03-31 04:23,47.74,17.4,935,15,false
 

src/syncfit/mcmc.py

@@ -11,7 +11,8 @@
 
 def do_emcee(theta_init:list[float], nu:list[float], F_muJy:list[float],
              F_error:list[float], model:BaseModel=B5, niter:int=2000,
-             nwalkers:int=100, fix_p:float=None, day:str=None, plot:bool=False
+             nwalkers:int=100, fix_p:float=None, upperlimits:list[bool]=None,
+             day:str=None, plot:bool=False
              ) -> tuple[list[float],list[float]]:
     """
     Fit the data with the given model using the emcee package.
@@ -27,6 +28,7 @@ def do_emcee(theta_init:list[float], nu:list[float], F_muJy:list[float],
         nwalkers (int): The number of walkers to use for emcee
         fix_p (float): Will fix the p value to whatever you give, do not provide p in theta_init
                                if this is the case!
+        upperlimits (list[bool]): True if the point is an upperlimit, False otherwise.
         day (string): day of observation, used for labeling plots
         plot (bool): If True, generate the plots used for debugging. Default is False.
     
@@ -43,8 +45,15 @@ def do_emcee(theta_init:list[float], nu:list[float], F_muJy:list[float],
     ndim = len(theta_init)
 
     # get some values from the import
+    nu = np.array(nu)
+    F_muJy = np.array(F_muJy)
+    F_error = np.array(F_error)
+    if upperlimits is not None:
+        upperlimits = np.array(upperlimits)
+
     pos, labels, emcee_args = model.unpack_util(theta_init, nu, F_muJy, F_error,
-                                                nwalkers, p=fix_p)
+                                                nwalkers, p=fix_p,
+                                                upperlimit=upperlimits)
 
     # setup and run the MCMC
     nwalkers, ndim = pos.shape

src/syncfit/models/b1b2_b3b4_weighted_model.py

@@ -57,14 +57,18 @@ def SED(nu, p, log_F_nu, log_nu_a, log_nu_m):
 
         return F
 
-    def lnprior(theta, p=None, **kwargs):
+    def lnprior(theta, nu, F, upperlimit, p=None, **kwargs):
         ''' Priors: '''
+        uppertest = BaseModel._is_below_upperlimits(
+            nu, F, upperlimit, theta, B1B2_B3B4_Weighted.SED, p=p
+        )
+
         if p is None:
             p, log_F_nu, log_nu_a, log_nu_m= theta
         else:
             log_F_nu, log_nu_a, log_nu_m= theta
 
-        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 12 and 6 < log_nu_m < 12:
+        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 12 and 6 < log_nu_m < 12 and uppertest:
             return 0.0
 
         else:

src/syncfit/models/b1b2_model.py

@@ -34,14 +34,18 @@ def SED(nu, p, log_F_nu, log_nu_a, log_nu_m):
 
         return F_nu * term1 * term2
 
-    def lnprior(theta, p=None, **kwargs):
+    def lnprior(theta, nu, F, upperlimit, p=None, **kwargs):
         ''' Priors: '''
+        uppertest = BaseModel._is_below_upperlimits(
+            nu, F, upperlimit, theta, B1B2.SED, p=p
+        )
+
         if p is None:
             p, log_F_nu, log_nu_a, log_nu_m = theta
         else:
             log_F_nu, log_nu_a, log_nu_m = theta
 
-        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 12 and log_nu_m > log_nu_a:
+        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 12 and log_nu_m > log_nu_a and uppertest:
             return 0.0
 
         else:

src/syncfit/models/b4b5_model.py

@@ -35,14 +35,18 @@ def SED(nu, p, log_F_nu, log_nu_a, log_nu_m):
 
         return F_nu * term1 * term2
 
-    def lnprior(theta, p=None, **kwargs):
+    def lnprior(theta, nu, F, upperlimit, p=None, **kwargs):
         ''' Priors: '''
+        uppertest = BaseModel._is_below_upperlimits(
+            nu, F, upperlimit, theta, B4B5.SED, p=p
+        )
+
         if p is None:
             p, log_F_nu, log_nu_a, log_nu_m = theta
         else:
             log_F_nu, log_nu_a, log_nu_m = theta
 
-        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and log_nu_m < log_nu_a:
+        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and log_nu_m < log_nu_a and uppertest:
             return 0.0
 
         else:

src/syncfit/models/b4b5b3_model.py

@@ -39,13 +39,17 @@ def SED(nu, p, log_F_nu, log_nu_a, log_nu_m, log_nu_c):
 
         return F_nu * term1 * term2 * term3
 
-    def lnprior(theta, p=None, **kwargs):
+    def lnprior(theta, nu, F, upperlimit, p=None, **kwargs):
         ''' Priors: '''
+        uppertest = BaseModel._is_below_upperlimits(
+            nu, F, upperlimit, theta, B4B5B3.SED, p=p
+        )
+
         if p is None:
             p, log_F_nu, log_nu_a, log_nu_m, log_nu_c= theta
         else:
             log_F_nu, log_nu_a, log_nu_m, log_nu_c= theta
-        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and log_nu_m < log_nu_a and log_nu_a < log_nu_c:
+        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and log_nu_m < log_nu_a and log_nu_a < log_nu_c and uppertest:
             return 0.0
 
         else:

src/syncfit/models/b5_model.py

@@ -28,13 +28,18 @@ def SED(nu, p, log_F_nu, log_nu_a):
 
         return F_nu*term**(-1/s)
 
-    def lnprior(theta, p=None, **kwargs):
+    def lnprior(theta, nu, F, upperlimit, p=None, **kwargs):
         ''' Priors: '''
+        uppertest = BaseModel._is_below_upperlimits(
+            nu, F, upperlimit, theta, B5.SED, p=p
+        )
+
         if p is None:
-            p, log_F_nu, log_nu_a= theta
+            p, log_F_nu, log_nu_a = theta
         else:
             log_F_nu, log_nu_a = theta
-        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11:
+
+        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and uppertest:
             return 0.0
 
         else:

src/syncfit/models/b5b3_model.py

@@ -36,14 +36,18 @@ def SED(nu, p, log_F_nu, log_nu_a, log_nu_c):
 
         return F_nu * term1 * term2
 
-    def lnprior(theta, p=None, **kwargs):
+    def lnprior(theta, nu, F, upperlimit, p=None, **kwargs):
         ''' Priors: '''
+        uppertest = BaseModel._is_below_upperlimits(
+            nu, F, upperlimit, theta, B5B3.SED, p=p
+        )
+
         if p is None:
             p, log_F_nu, log_nu_a, log_nu_c = theta
         else:
             log_F_nu, log_nu_a, log_nu_c = theta
 
-        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and log_nu_c > log_nu_a:
+        if 2< p < 4 and -4 < log_F_nu < 2 and 6 < log_nu_a < 11 and log_nu_c > log_nu_a and uppertest:
             return 0.0
 
         else:

src/syncfit/models/base_model.py

@@ -49,7 +49,7 @@ def get_pos(theta_init:list, nwalkers:int) -> list[float]:
 
         return pos
 
-    def get_emcee_args(nu:list, F_muJy:list, F_error:list, p:float=None) -> dict:
+    def get_emcee_args(nu:list, F_muJy:list, F_error:list, p:float=None, upperlimit:list=None) -> dict:
         '''
         Packages up the args to be passed into the model based on the user input.
 
@@ -67,13 +67,13 @@ def get_emcee_args(nu:list, F_muJy:list, F_error:list, p:float=None) -> dict:
         F_error = np.array(F_error)*1e-3
 
         if p is None:
-            return {'nu':nu, 'F':F, 'F_error':F_error} 
+            return {'nu':nu, 'F':F, 'F_error':F_error, 'upperlimit':upperlimit} 
         else:
-            return {'nu':nu, 'F':F, 'F_error':F_error, 'p':p} 
+            return {'nu':nu, 'F':F, 'F_error':F_error, 'p':p, 'upperlimit':upperlimit} 
 
     # package those up for easy getting in do_emcee
     @classmethod
-    def unpack_util(cls, theta_init, nu, F_muJy, F_error, nwalkers, p=None):
+    def unpack_util(cls, theta_init, nu, F_muJy, F_error, nwalkers, p=None, upperlimit=None):
         '''
         A wrapper on the utility functions.
 
@@ -87,7 +87,7 @@ def unpack_util(cls, theta_init, nu, F_muJy, F_error, nwalkers, p=None):
         '''
         return (cls.get_pos(theta_init,nwalkers),
                 cls.get_labels(p=p),
-                cls.get_emcee_args(nu, F_muJy, F_error, p))
+                cls.get_emcee_args(nu, F_muJy, F_error, p, upperlimit))
 
     @classmethod
     def lnprob(cls, theta:list, **kwargs):
@@ -130,7 +130,26 @@ def loglik(cls, theta, nu, F, F_error, p=None, **kwargs):
         chi2 = np.sum((F - model_result)**2/sigma2)
         ll = -0.5*chi2
         return ll
-
+
+    @staticmethod
+    def _is_below_upperlimits(nu, F, upperlimits, theta, model, p=None):
+        '''
+        Checks that the location of theta is below any upperlimits
+        '''
+
+        if upperlimits is None:
+            return True
+
+        where_upperlimit = np.where(upperlimits)[0]
+        F_upperlimits = F[where_upperlimit]
+
+        if p is None:
+            test_fluxes = model(nu, *theta)[where_upperlimit]
+        else:
+            test_fluxes = model(nu, p, *theta)[where_upperlimit]
+
+        return np.all(F_upperlimits > test_fluxes)
+
     # Some *required* abstract methods
     @staticmethod
     @abstractmethod