add more tests

src/adler/science/PhaseCurve.py

@@ -12,6 +12,14 @@
     "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:
@@ -80,15 +88,45 @@ def __init__(
         else:
             print("no model selected")
 
-        ### set bounds here by default using SetModelBounds?
-        # for x in model phase parameters:
-        #   self.SetModelBounds(x)
-
     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
@@ -127,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
         -----------
@@ -155,17 +193,14 @@ def InitModelSbpy(self, model_sbpy):
         for p in param_names:
             if p in model_sbpy.__dict__:  # check that the parameter is available in the sbpy object
                 x = getattr(model_sbpy, p)
-                print(p, x)
                 # try get the quantity (value with units)
                 if hasattr(x, "unit"):
-                    print(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
-                    print(x)
                 # 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"
@@ -220,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.
 
         Returns
         ----------
@@ -233,7 +269,6 @@ 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)

tests/adler/science/test_PhaseCurve.py

@@ -1,11 +1,12 @@
+from adler.science.PhaseCurve import PhaseCurve
 from numpy.testing import assert_array_equal
 import pytest
+import numpy as np
+import astropy.units as u
 
 
 def test_PhaseCurve_init():
-    import numpy as np
-    import astropy.units as u
-    from adler.science.PhaseCurve import PhaseCurve
+    """Test intialising the model."""
 
     # test creating a model
     H = 18.9
@@ -19,9 +20,7 @@ def test_PhaseCurve_init():
 
 
 def test_PhaseCurve_ReducedMag():
-    import numpy as np
-    import astropy.units as u
-    from adler.science.PhaseCurve import PhaseCurve
+    """Test calculating the reduced magnitude from a PhaseCurve object."""
 
     # define the phase angles
     alpha = np.array([0, 10]) * u.deg
@@ -36,10 +35,8 @@ def test_PhaseCurve_ReducedMag():
     assert_array_equal(red_mag.value, np.array([18.0, 19.0]))
 
 
-def test_PhaseCurve_FitModel():
-    import numpy as np
-    import astropy.units as u
-    from adler.science.PhaseCurve import PhaseCurve
+def test_PhaseCurve_FitModel_linear():
+    """Test fitting a linear phase function to two data points."""
 
     # define the observations
     alpha = np.array([0, 10]) * u.deg
@@ -48,13 +45,79 @@ def test_PhaseCurve_FitModel():
     # empty linear phase curve model
     pc_lin = PhaseCurve(model_name="LinearPhaseFunc")
 
-    # fit the model to the data
+    # fit the model to the data - this returns an sbpy object
     pc_fit = pc_lin.FitModel(alpha, red_mag)
 
     assert pc_fit.H.unit == u.mag
     assert pc_fit.H.value == 18.0
     assert pc_fit.S.unit == u.mag / u.deg
     assert pc_fit.S.value == 0.1
+    assert (
+        hasattr(pc_fit, "H_err") == False
+    )  # with only two data points the covariance matrix is None and no uncertainties are stored
 
 
-### TODO - add more phase curve tests, test bounds, fixed values and fit uncertainties
+def test_PhaseCurve_FitModel_HG():
+    """Test fitting a HG model to generated data."""
+
+    # generate some model data
+    pc1 = PhaseCurve(H=18.0 * u.mag, phase_parameter_1=0.15, model_name="HG")
+    alpha = np.linspace(0, 30) * u.deg
+    red_mag = pc1.ReducedMag(alpha)
+
+    # fit the same phase curve model to the data
+    pc_fit = pc1.FitModel(alpha, red_mag)
+    # convert from sbpy to adler PhaseCurve object
+    pc2 = pc1.InitModelSbpy(pc_fit)
+
+    # the new fitted model should have the same parameters as the input model
+    assert pc2.H == pc1.H
+    assert pc2.phase_parameter_1 == pc1.phase_parameter_1
+    assert pc2.phase_parameter_2 is None
+    assert pc1.phase_parameter_1_err is None  # the first model had no uncertainties
+    assert pc2.phase_parameter_1_err is not None  # the fitted model has some uncertainties
+
+
+def test_PhaseCurve_FitModel_HG_fixed():
+    """Test fitting a just H whilst keeping G fixed."""
+
+    # generate some model data
+    pc1 = PhaseCurve(H=18.0 * u.mag, phase_parameter_1=0.15, model_name="HG")
+    alpha = np.linspace(0, 30) * u.deg
+    red_mag = pc1.ReducedMag(alpha)
+
+    # fix phase_parameter_1
+    pc1.FixParam("phase_parameter_1")
+
+    # fit the same phase curve model to the data
+    pc_fit = pc1.FitModel(alpha, red_mag)
+    # convert from sbpy to adler PhaseCurve object
+    pc2 = pc1.InitModelSbpy(pc_fit)
+
+    # the new fitted model should have the same parameters as the input model, but G is fixed
+    assert pc_fit._fixed["G"] is True
+    assert pc2.H == pc1.H
+    assert pc2.phase_parameter_1 == pc1.phase_parameter_1
+    assert pc2.phase_parameter_1_err is None  # the fitted model has no uncertainties when param is fixed
+
+
+def test_PhaseCurve_FitModel_HG_bounds():
+    """Test fitting a just H whilst keeping G fixed."""
+
+    # generate some model data
+    pc1 = PhaseCurve(H=18.0 * u.mag, phase_parameter_1=0.15, model_name="HG")
+    alpha = np.linspace(0, 30) * u.deg
+    red_mag = pc1.ReducedMag(alpha)
+
+    # set bounds on phase parameter
+    pc1.SetModelBounds("phase_parameter_1", (0.0, 0.1))
+
+    # fit the same phase curve model to the data
+    pc_fit = pc1.FitModel(alpha, red_mag)
+    # convert from sbpy to adler PhaseCurve object
+    pc2 = pc1.InitModelSbpy(pc_fit)
+
+    # the new fitted model should have the same parameters as the input model, but G is fixed
+    assert pc_fit.G.bounds == (0.0, 0.1)
+    assert pc2.phase_parameter_1 == 0.1
+    assert pc2.phase_parameter_1_err is not None