Merge pull request #148 from lsst-uk/147-alderdata-fields

147 alderdata fields

src/adler/adler.py

@@ -100,8 +100,8 @@ def runAdler(cli_args):
 
             # initial simple phase curve filter model with fixed G12
             pc = PhaseCurve(
-                abs_mag=sso.H * u.mag,
-                phase_param=0.62,
+                H=sso.H * u.mag,
+                phase_parameter_1=0.62,
                 model_name="HG12_Pen16",
             )
 
@@ -168,15 +168,15 @@ def runAdler(cli_args):
 
             # if sum(obs_mask) < N_pc_fit:
             #     # use an assumed value of G12 until more data is available
-            #     pc_fit = PhaseCurve(abs_mag=sso.H * u.mag, phase_param=0.62, model_name="HG12_Pen16")
+            #     pc_fit = PhaseCurve(H=sso.H * u.mag, phase_parameter_1=0.62, model_name="HG12_Pen16")
             # else:
             #     # do a simple HG12_Pen16 fit to the past data
             #     pc_fit = pc.FitModel(alpha[obs_mask], red_mag[obs_mask], mag_err[obs_mask])
             #     pc_fit = pc.InitModelSbpy(pc_fit)
 
             # print(pc_fit)
-            # print(pc_fit.abs_mag)
-            # print(pc_fit.phase_param)
+            # print(pc_fit.H)
+            # print(pc_fit.phase_parameter_1)
 
             # # now check if the new observations are outlying
             # alpha_new = alpha[~obs_mask]

src/adler/science/PhaseCurve.py

@@ -1,7 +1,26 @@
 from sbpy.photometry import HG, HG1G2, HG12, HG12_Pen16, LinearPhaseFunc
 import astropy.units as u
+import numpy as np
 from astropy.modeling.fitting import LevMarLSQFitter
 
+# translation between sbpy and adler field names
+SBPY_ADLER_DICT = {
+    "H": "H",
+    "G": "phase_parameter_1",
+    "G12": "phase_parameter_1",
+    "G1": "phase_parameter_1",
+    "G2": "phase_parameter_2",
+    "S": "phase_parameter_1",
+}
+# translation between adler and sbpy fields, depends on model
+ADLER_SBPY_DICT = {
+    "HG": {"phase_parameter_1": "G"},
+    "HG1G2": {"phase_parameter_1": "G1", "phase_parameter_2": "G2"},
+    "HG12": {"phase_parameter_1": "G12"},
+    "HG12_Pen16": {"phase_parameter_1": "G12"},
+    "LinearPhaseFunc": {"phase_parameter_1": "S"},
+}
+
 
 class PhaseCurve:
     """A class to define the phasecurve model and associated functions.
@@ -10,14 +29,14 @@ class PhaseCurve:
 
     Attributes
     -----------
-    abs_mag : float
+    H : float
        Absolute magnitude, H, of the phasecurve model (often units of mag).
 
-    phase_param: float
+    phase_parameter_1: float
        The first phase parameter of the phasecurve model, e.g. G from HG
        (often dimensionless units unless S from LinearPhaseFunc, which has units mag/deg or mag/rad).
 
-    phase_param2: float
+    phase_parameter_2: float
        The second phase parameter, only used for the 3 parameter HG1G2 phasecurve model.
 
     model_name: str
@@ -27,32 +46,87 @@ class PhaseCurve:
     bounds: bool
        Flag for using the default bounds on the sbpy model (True) or ignoring bounds (False).
 
+    H_err : float
+       Uncertainty in absolute magnitude.
+
+    phase_parameter_1_err : float
+       Uncertainty in the (first) phase parameter.
+
+    phase_parameter_2_err : float
+       Uncertainty in the (second, if required) phase parameter.
+
     """
 
-    def __init__(self, abs_mag=18, phase_param=0.2, phase_param2=None, model_name="HG"):
-        self.abs_mag = abs_mag
-        self.phase_param = phase_param
-        self.phase_param2 = phase_param2
+    def __init__(
+        self,
+        H=18,
+        phase_parameter_1=0.2,
+        phase_parameter_2=None,
+        model_name="HG",
+        H_err=None,
+        phase_parameter_1_err=None,
+        phase_parameter_2_err=None,
+    ):
+        self.H = H
+        self.phase_parameter_1 = phase_parameter_1
+        self.phase_parameter_2 = phase_parameter_2
         self.model_name = model_name
+        self.H_err = H_err
+        self.phase_parameter_1_err = phase_parameter_1_err
+        self.phase_parameter_2_err = phase_parameter_2_err
 
         if model_name == "HG":
-            self.model_function = HG(H=abs_mag, G=self.phase_param)
+            self.model_function = HG(H=H, G=self.phase_parameter_1)
         elif model_name == "HG1G2":
-            self.model_function = HG1G2(H=abs_mag, G1=self.phase_param, G2=self.phase_param)
+            self.model_function = HG1G2(H=H, G1=self.phase_parameter_1, G2=self.phase_parameter_1)
         elif model_name == "HG12":
-            self.model_function = HG12(H=abs_mag, G12=self.phase_param)
+            self.model_function = HG12(H=H, G12=self.phase_parameter_1)
         elif model_name == "HG12_Pen16":
-            self.model_function = HG12_Pen16(H=abs_mag, G12=self.phase_param)
+            self.model_function = HG12_Pen16(H=H, G12=self.phase_parameter_1)
         elif model_name == "LinearPhaseFunc":
-            self.model_function = LinearPhaseFunc(H=abs_mag, S=self.phase_param)
+            self.model_function = LinearPhaseFunc(H=H, S=self.phase_parameter_1)
         else:
             print("no model selected")
 
     def SetModelBounds(self, param, bound_vals=(None, None)):
+        """Set the "bounds" attribute of an sbpy model parameter, i.e. the lower and upper constraints for the fitter.
+
+        Parameters
+        -----------
+        param : str
+           Parameter name bounds to be set fix, if not an sbpy parameter name (e.g. G) the corresponding adler name is looked up (e.g. phase_parameter_1)
+
+        bound_vals : tuple
+           Set the fitter constraints to (upper, lower) for the param.
+
+        """
         model_sbpy = self.model_function
-        param_names = model_sbpy.param_names
+        if param not in model_sbpy.__dict__:
+            param = ADLER_SBPY_DICT[self.model_name][param]
         x = getattr(model_sbpy, param)
         setattr(x, "bounds", bound_vals)
+        return
+
+    def FixParam(self, param, fix_flag=True):
+        """Set the "fixed" attribute of an sbpy model parameter.
+        E.g. use this to fit for absolute magnitude whilst keeping phase parameter fixed.
+
+        Parameters
+        -----------
+        param : str
+           Parameter name to fix, if not an sbpy parameter name (e.g. G) the corresponding adler name is looked up (e.g. phase_parameter_1)
+
+        fix_flag : bool
+           Set True to keep the param fixed when fitting
+
+        """
+
+        model_sbpy = self.model_function
+        if param not in model_sbpy.__dict__:
+            param = ADLER_SBPY_DICT[self.model_name][param]
+        x = getattr(model_sbpy, param)
+        setattr(x, "fixed", fix_flag)
+        return
 
     def ReturnModelDict(self):
         """Return the values for the PhaseCurve class as a dict
@@ -74,7 +148,7 @@ def InitModelDict(self, model_dict):
         Parameters
         -----------
         model_dict : dict
-           Dictionary containing the PhaseCurve parameters you wish to set, e.g. abs_mag, phase_param
+           Dictionary containing the PhaseCurve parameters you wish to set, e.g. H, phase_parameter_1
 
         Returns
         ----------
@@ -91,7 +165,7 @@ def InitModelDict(self, model_dict):
 
     def InitModelSbpy(self, model_sbpy):
         """Set up a new PhaseCurve model object from an existing sbpy model
-        ### or create dict from sbpy model and then use InitModelDict?
+        ### TODO or create dict from sbpy model and then use InitModelDict?
 
         Parameters
         -----------
@@ -109,19 +183,33 @@ def InitModelSbpy(self, model_sbpy):
         model_name = model_sbpy.__class__.name
 
         # get the sbpy model parameters
-        param_names = model_sbpy.param_names
-        parameters = []
+        param_names = list(model_sbpy.param_names)
+        # we will also check for any uncertainties we have stored in the sbpy object
+        param_names_err = ["{}_err".format(x) for x in param_names]
+        param_names = param_names + param_names_err
+
+        # create a dictionary of phase curve parameters from sbpy in a format accepted by PhaseCurve
+        parameters = {}
         for p in param_names:
-            # try get the quantity (value with units)
-            x = getattr(model_sbpy, p).quantity
-            # if there are no units get just the value
-            if x is None:
-                x = getattr(model_sbpy, p).value
-            parameters.append(x)
-        # print(param_names, parameters)
+            if p in model_sbpy.__dict__:  # check that the parameter is available in the sbpy object
+                x = getattr(model_sbpy, p)
+                # try get the quantity (value with units)
+                if hasattr(x, "unit"):
+                    if (x.unit is None) or (
+                        x.unit == ""
+                    ):  # if there are no units (or weird blank units?) get just the value
+                        x = x.value
+                    else:
+                        x = x.quantity
+                # look up the correct adler parameter name (accounting for additional uncertainty, "_err", parameters)
+                if p.endswith("_err"):  # assumes the uncertainty parameter always ends in "_err"
+                    _p = SBPY_ADLER_DICT[p.split("_err")[0]] + "_err"
+                    parameters[_p] = x
+                else:
+                    parameters[SBPY_ADLER_DICT[p]] = x
 
         # create a PhaseCurve object with the extracted parameters
-        model = PhaseCurve(*parameters, model_name=model_name)
+        model = PhaseCurve(**parameters, model_name=model_name)
 
         return model
 
@@ -167,7 +255,8 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None):
            Uncertainty on the reduced magnitude, used to weight the fit.
 
         fitter : object
-           Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter
+           Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter.
+           N.B. that LevMarLSQFitter cannot handle inequality constraints for the HG1G2 model, use something like SLSQPLSQFitter from astropy.modeling.fitting (does not return covariance matrix!).
 
         Returns
         ----------
@@ -180,14 +269,34 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None):
         # use the LevMarLSQFitter by default
         if fitter is None:
             fitter = LevMarLSQFitter()
-        # print(fitter)
 
         if mag_err is not None:  # fit weighted by photometric uncertainty
             model_fit = fitter(self.model_function, phase_angle, reduced_mag, weights=1.0 / mag_err)
         else:  # unweighted fit
             model_fit = fitter(self.model_function, phase_angle, reduced_mag)
 
+        # Add fitted uncertainties as an additional attribute within the sbpy object
+        if "param_cov" in fitter.fit_info:
+            # get the covariance matrix from the fit
+            covariance = fitter.fit_info["param_cov"]
+            if covariance is not None:
+                # get fit uncertainties as square of the diagonal of the covariance
+                fit_errs = np.sqrt(np.diag(covariance))
+                # update only the uncertainties for parameters used in the fit
+                param_names = np.array(model_fit.param_names)
+                fit_mask = ~np.array([getattr(model_fit, x).fixed for x in param_names])
+                for i, x in enumerate(param_names[fit_mask]):
+                    setattr(model_fit, "{}_err".format(x), fit_errs[i])
+            # else:
+            ### TODO log covariance is None error here
+
+        ### TODO
+        # else:
+        #     log lack of uncertainties for fitter
+        #     run an MCMC estimate of uncertainty?
+
         ### if overwrite_model: # add an overwrite option?
         # redo __init__ with the new fitted parameters
+        # this would then return an adler PhaseCurve object rather than an sbpy object
 
         return model_fit