From c1b1c8837fca58af92eb54ada7e17dd3fb341c96 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Thu, 6 Jun 2024 14:35:02 +0100 Subject: [PATCH] add more tests --- src/adler/science/PhaseCurve.py | 57 +++++++++++++---- tests/adler/science/test_PhaseCurve.py | 87 ++++++++++++++++++++++---- 2 files changed, 121 insertions(+), 23 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index c4bfbec..99fc16d 100644 --- a/src/adler/science/PhaseCurve.py +++ b/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) diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index b3c0ceb..4cbf76b 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/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