diff --git a/src/adler/adler_run.py b/src/adler/adler_run.py index 181fd6a..3be24fa 100644 --- a/src/adler/adler_run.py +++ b/src/adler/adler_run.py @@ -1,6 +1,7 @@ import logging import argparse import astropy.units as u +from astropy.time import Time import numpy as np import pandas as pd import matplotlib.pyplot as plt @@ -103,7 +104,8 @@ def runAdler(cli_args): # initial simple phase curve filter model with fixed G12 pc = PhaseCurve( - H=sso.H * u.mag, + # H=sso.H * u.mag, + H=sso.H, phase_parameter_1=0.62, model_name="HG12_Pen16", ) @@ -118,14 +120,24 @@ def runAdler(cli_args): # do a HG12_Pen16 fit to the past data pc_fit = pc.FitModel( - np.array(df_obs["phaseAngle"]) * u.deg, - np.array(df_obs["reduced_mag"]) * u.mag, - np.array(df_obs["magErr"]) * u.mag, + # np.array(df_obs["phaseAngle"]) * u.deg, + # np.array(df_obs["reduced_mag"]) * u.mag, + # np.array(df_obs["magErr"]) * u.mag, + np.radians(np.array(df_obs["phaseAngle"])), + np.array(df_obs["reduced_mag"]), + np.array(df_obs["magErr"]), ) pc_fit = pc.InitModelSbpy(pc_fit) - # Store the fitted values in an AdlerData object - adler_data.populate_phase_parameters(filt, **pc_fit.__dict__) + # Store the fitted values, and metadata, in an AdlerData object + ad_params = pc_fit.__dict__ + ad_params["phaseAngle_min"] = np.amin(df_obs["phaseAngle"]) # * u.deg + ad_params["phaseAngle_range"] = np.ptp(df_obs["phaseAngle"]) # * u.deg + ad_params["arc"] = np.ptp(df_obs["midPointMjdTai"]) # * u.d + ad_params["nobs"] = len(df_obs) + ad_params["modelFitMjd"] = Time.now().mjd + # adler_data.populate_phase_parameters(filt, **pc_fit.__dict__) + adler_data.populate_phase_parameters(filt, **ad_params) # add to plot ax1 = fig.axes[0] @@ -134,13 +146,18 @@ def runAdler(cli_args): alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg ax1.plot( alpha.value, - pc_fit.ReducedMag(alpha).value, - label="{}, H={:.2f}, G12={:.2f}".format(filt, pc_fit.H.value, pc_fit.phase_parameter_1), + # pc_fit.ReducedMag(alpha).value, + # label="{}, H={:.2f}, G12={:.2f}".format(filt, pc_fit.H.value, pc_fit.phase_parameter_1), + pc_fit.ReducedMag(alpha), + label="{}, H={:.2f}, G12={:.2f}".format(filt, pc_fit.H, pc_fit.phase_parameter_1), ) - - # TODO: save the figures if an outpath is provided ax1.legend() - if cli_args.outpath: + + # TODO: Use a CLI arg flag to open figure interactively instead of saving? + if cli_args.plot_show: + plt.show() + # Save figures at the outpath location + else: fig_file = "{}/phase_curve_{}_{}.png".format( cli_args.outpath, cli_args.ssObjectId, int(np.amax(df_obs["midPointMjdTai"])) ) @@ -149,11 +166,15 @@ def runAdler(cli_args): logger.info(msg) fig = plot_errorbar(planetoid, fig=fig, filename=fig_file) # TODO: add titles with filter name? plt.close() - else: - plt.show() - # TODO: output adler values to a database + # Output adler values to a database if a db_name is provided print(adler_data.__dict__) + if cli_args.db_name: + adler_db = "{}/{}".format(cli_args.outpath, cli_args.db_name) + msg = "write to {}".format(adler_db) + print(msg) + logger.info(msg) + adler_data.write_row_to_database(adler_db) # analyse colours for the filters provided logger.info("Calculate colours: {}".format(cli_args.colour_list)) @@ -183,10 +204,10 @@ def runAdler(cli_args): # determine the filt_obs - filt_ref colour # generate a plot - if cli_args.outpath: - plot_dir = cli_args.outpath - else: + if cli_args.plot_show: plot_dir = None + else: + plot_dir = cli_args.outpath col_dict = col_obs_ref( planetoid, @@ -198,6 +219,7 @@ def runAdler(cli_args): # x1 = x1, # x2 = x2, plot_dir=plot_dir, + plot_show=cli_args.plot_show, ) print(col_dict) @@ -253,17 +275,27 @@ def main(): optional_group.add_argument( "-n", "--db_name", - help="Stem filename of output database. If this doesn't exist, it will be created. Default: adler_out.", + # help="Stem filename of output database. If this doesn't exist, it will be created. Default: adler_out.", + # type=str, + # default="adler_out", + help="Optional filename of output database, used to store Adler results in a db if provided.", type=str, - default="adler_out", + default=None, ) optional_group.add_argument( "-i", "--sql_filename", - help="Optional input path location of a sql database file containing observations", + help="Optional input path location of a sql database file containing observations.", type=str, default=None, ) + # TODO: add flag argument to display plots instead of saving them + optional_group.add_argument( + "-p", + "--plot_show", + help="Optional flag to display plots interactively instead of saving to file.", + action="store_true", + ) args = parser.parse_args() diff --git a/src/adler/objectdata/AdlerData.py b/src/adler/objectdata/AdlerData.py index b8db020..d1f94f1 100644 --- a/src/adler/objectdata/AdlerData.py +++ b/src/adler/objectdata/AdlerData.py @@ -7,7 +7,7 @@ from datetime import datetime, timezone -FILTER_DEPENDENT_KEYS = ["phaseAngle_min", "phaseAngle_range", "nobs", "arc", "modelFitMjd"] +FILTER_DEPENDENT_KEYS = ["phaseAngle_min", "phaseAngle_range", "nobs", "arc"] MODEL_DEPENDENT_KEYS = [ "H", "H_err", @@ -15,6 +15,7 @@ "phase_parameter_1_err", "phase_parameter_2", "phase_parameter_2_err", + "modelFitMjd", ] ALL_FILTER_LIST = ["u", "g", "r", "i", "z", "y"] @@ -522,6 +523,7 @@ class PhaseModelDependentAdler: phase_parameter_1_err: float = np.nan phase_parameter_2: float = np.nan phase_parameter_2_err: float = np.nan + modelFitMjd: float = np.nan class PhaseParameterOutput: diff --git a/src/adler/science/Colour.py b/src/adler/science/Colour.py index 37c2d58..fb4e0fd 100644 --- a/src/adler/science/Colour.py +++ b/src/adler/science/Colour.py @@ -19,6 +19,7 @@ def col_obs_ref( x1=None, x2=None, plot_dir=None, + plot_show=False, ): """A function to calculate the colour of an Adler planetoid object. An observation in a given filter (filt_obs) is compared to previous observation in a reference filter (filt_ref). @@ -99,7 +100,8 @@ def col_obs_ref( x_obs = df_obs.iloc[i][x_col] y_obs = df_obs.iloc[i][y_col] yerr_obs = df_obs.iloc[i][yerr_col] - obsId = df_obs.iloc[i][obsId_col] + # obsId = df_obs.iloc[i][obsId_col] # NB: for some reason iloc here doesn't preserve the int64 dtype of obsId_col...? + obsId = df_obs[obsId_col].iloc[i] # select observations in the reference filter from before the new obs ref_mask = df_obs_ref[x_col] < x_obs @@ -112,11 +114,11 @@ def col_obs_ref( # select only the N_ref ref obs for comparison _df_obs_ref = df_obs_ref[ref_mask].iloc[-_N_ref:] - print(len(_df_obs_ref)) - print(np.array(_df_obs_ref[x_col])) + # print(len(_df_obs_ref)) + # print(np.array(_df_obs_ref[x_col])) if len(_df_obs_ref) == 0: - print("no reference observations") # TODO: add proper error handling and logging here - return df_obs + # print("insufficient reference observations") # TODO: add proper error handling and logging here + return # determine reference observation values y_ref = np.mean(_df_obs_ref[y_col]) # TODO: add option to choose statistic, e.g. mean or median? @@ -127,7 +129,9 @@ def col_obs_ref( # Create the colour dict col_dict = {} - col_dict[obsId_col] = obsId + col_dict[obsId_col] = np.int64( + obsId + ) # store id as an int to make sure it doesn't get stored as float e notation! col_dict[x_col] = x_obs col_dict[colour] = y_obs - y_ref col_dict[delta_t_col] = x_obs - x2_ref @@ -141,7 +145,7 @@ def col_obs_ref( # need to test error case where there are no r filter obs yet # TODO: add a plotting option? - if plot_dir: + if plot_dir or plot_show: fig = plt.figure() gs = gridspec.GridSpec(1, 1) ax1 = plt.subplot(gs[0, 0]) @@ -160,13 +164,16 @@ def col_obs_ref( ax1.legend() ax1.invert_yaxis() - fname = "{}/colour_plot_{}_{}-{}_{}.png".format( - plot_dir, planetoid.ssObjectId, filt_obs, filt_ref, int(x_obs) - ) - print("Save figure: {}".format(fname)) - plt.savefig(fname, facecolor="w", transparent=True, bbox_inches="tight") - - # plt.show() # TODO: add option to display figure, or to return the fig object? - plt.close() + if plot_dir: + fname = "{}/colour_plot_{}_{}-{}_{}.png".format( + plot_dir, planetoid.ssObjectId, filt_obs, filt_ref, int(x_obs) + ) + print("Save figure: {}".format(fname)) + plt.savefig(fname, facecolor="w", transparent=True, bbox_inches="tight") + + if plot_show: + plt.show() # TODO: add option to display figure, or to return the fig object? + else: + plt.close() return col_dict diff --git a/src/adler/utilities/AdlerCLIArguments.py b/src/adler/utilities/AdlerCLIArguments.py index b3e2980..3a717b3 100644 --- a/src/adler/utilities/AdlerCLIArguments.py +++ b/src/adler/utilities/AdlerCLIArguments.py @@ -25,6 +25,8 @@ def __init__(self, args): self.outpath = args.outpath self.db_name = args.db_name self.sql_filename = args.sql_filename + self.plot_show = args.plot_show + self.phase_model = args.phase_model self.validate_arguments() @@ -47,6 +49,9 @@ def validate_arguments(self): if self.colour_list: self._validate_colour_list() + if self.phase_model: + self._validate_phase_model() + def _validate_filter_list(self): """Validation checks for the filter_list command-line argument.""" expected_filters = ["u", "g", "r", "i", "z", "y"] @@ -149,3 +154,16 @@ def _validate_sql_filename(self): raise ValueError( "The file supplied for the command-line argument --sql_filename cannot be found." ) + + def _validate_phase_model(self): + """Validation checks for the phase_model command-line argument.""" + expected_models = ["HG", "HG1G2", "HG12", "HG12_Pen16", "LinearPhaseFunc"] + err_msg_model = ( + "Unexpected model in --phase_model command-line arguments. Please select from {}".format( + expected_models + ) + ) + + if self.phase_model not in expected_models: + logging.error(err_msg_model) + raise ValueError(err_msg_model) diff --git a/src/adler/utilities/science_utilities.py b/src/adler/utilities/science_utilities.py index 90ff1df..146f9c0 100644 --- a/src/adler/utilities/science_utilities.py +++ b/src/adler/utilities/science_utilities.py @@ -78,7 +78,8 @@ def zero_func(x, axis=None): return 0 -def sigma_clip(data_res, kwargs={"maxiters": 1, "cenfunc": zero_func}): +# def sigma_clip(data_res, kwargs={"sigma":3, "maxiters": 1, "cenfunc": zero_func}): +def sigma_clip(data_res, **kwargs): """Wrapper function for astropy.stats.sigma_clip, here we define the default centre of the data (the data - model residuals) to be zero Parameters @@ -206,7 +207,9 @@ def get_df_obs_filt(planetoid, filt, x_col="midPointMjdTai", x1=None, x2=None, c if "AbsMag" in col_list: # calculate the model absolute magnitude # TODO: add robustness to the units, phaseAngle and reduced_mag must match pc_model - df_obs["AbsMag"] = pc_model.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value + # For now we must assume that there are no units, and that degrees have been passed... + # df_obs["AbsMag"] = pc_model.AbsMag(obs.phaseAngle * u.deg, obs.reduced_mag * u.mag).value + df_obs["AbsMag"] = pc_model.AbsMag(np.radians(obs.phaseAngle), obs.reduced_mag) # select only the required columns df_obs = df_obs[col_list] diff --git a/tests/adler/objectdata/test_AdlerData.py b/tests/adler/objectdata/test_AdlerData.py index dd4af55..4751521 100644 --- a/tests/adler/objectdata/test_AdlerData.py +++ b/tests/adler/objectdata/test_AdlerData.py @@ -4,219 +4,394 @@ import pandas as pd import sqlite3 -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_almost_equal from adler.objectdata.AdlerData import AdlerData from adler.utilities.tests_utilities import get_test_data_filepath # setting up the AdlerData object to be used for testing -test_object = AdlerData("8268570668335894776", ["u", "g", "r"]) - -u_model_1 = { - "model_name": "model_1", - "phaseAngle_min": 11.0, - "phaseAngle_range": 12.0, - "nobs": 13, - "arc": 14.0, - "H": 15.0, - "H_err": 16.0, - "phase_parameter_1": 17.0, - "phase_parameter_1_err": 18.0, -} -u_model_2 = { - "model_name": "model_2", - "H": 25.0, - "H_err": 26.0, - "phase_parameter_1": 27.0, - "phase_parameter_1_err": 28.0, - "phase_parameter_2": 29.0, - "phase_parameter_2_err": 30.0, -} +test_object = AdlerData("8268570668335894776", ["g", "r", "i"]) +model_1 = "HG12_Pen16" +model_2 = "HG" +# test data was generated in the tests/data dir using: +# python adler_run_test_data.py +# the output is copied here g_model_1 = { - "model_name": "model_1", - "phaseAngle_min": 31.0, - "phaseAngle_range": 32.0, - "nobs": 33, - "arc": 34.0, - "H": 35.0, - "H_err": 36.0, - "phase_parameter_1": 37.0, - "phase_parameter_1_err": 38.0, + "filter_name": "g", + "phaseAngle_min": 27.426668167114258, + "phaseAngle_range": 36.17388343811035, + "nobs": 9, + "arc": 3306.0079999999944, + "model_name": "HG12_Pen16", + "H": 20.719465178426468, + "H_err": 0.018444134023977106, + "phase_parameter_1": 0.62, + "phase_parameter_1_err": np.nan, + "phase_parameter_2": np.nan, + "phase_parameter_2_err": np.nan, +} +r_model_1 = { + "filter_name": "r", + "phaseAngle_min": 2.553332567214966, + "phaseAngle_range": 124.23803400993347, + "nobs": 38, + "arc": 3338.0655999999944, + "model_name": "HG12_Pen16", + "H": 19.92863542616601, + "H_err": 0.018525355171274356, + "phase_parameter_1": 1.0, + "phase_parameter_1_err": 0.05300829494059732, + "phase_parameter_2": np.nan, + "phase_parameter_2_err": np.nan, +} +i_model_1 = { + "filter_name": "i", + "phaseAngle_min": 10.0595064163208, + "phaseAngle_range": 101.74150371551514, + "nobs": 32, + "arc": 3342.0585599999977, + "model_name": "HG12_Pen16", + "H": 19.653155995892117, + "H_err": 0.01415178691419005, + "phase_parameter_1": 1.0, + "phase_parameter_1_err": 0.01516569963597136, + "phase_parameter_2": np.nan, + "phase_parameter_2_err": np.nan, +} +g_model_2 = { + "filter_name": "g", + "phaseAngle_min": 27.426668167114258, + "phaseAngle_range": 36.17388343811035, + "nobs": 9, + "arc": 3306.0079999999944, + "model_name": "HG", + "H": 20.72954332871528, + "H_err": 0.018444134023977106, + "phase_parameter_1": 0.15, + "phase_parameter_1_err": np.nan, + "phase_parameter_2": np.nan, + "phase_parameter_2_err": np.nan, } - r_model_2 = { - "model_name": "model_2", - "phaseAngle_min": 41.0, - "phaseAngle_range": 42.0, - "nobs": 43, - "arc": 44.0, - "H": 45.0, - "H_err": 46.0, - "phase_parameter_1": 47.0, - "phase_parameter_1_err": 48.0, - "phase_parameter_2": 49.0, - "phase_parameter_2_err": 50.0, + "filter_name": "r", + "phaseAngle_min": 2.553332567214966, + "phaseAngle_range": 124.23803400993347, + "nobs": 38, + "arc": 3338.0655999999944, + "model_name": "HG", + "H": 18.879873513575124, + "H_err": 0.007456882346021157, + "phase_parameter_1": -0.253, + "phase_parameter_1_err": 1.6701987135262256e-05, + "phase_parameter_2": np.nan, + "phase_parameter_2_err": np.nan, +} +i_model_2 = { + "filter_name": "i", + "phaseAngle_min": 10.0595064163208, + "phaseAngle_range": 101.74150371551514, + "nobs": 32, + "arc": 3342.0585599999977, + "model_name": "HG", + "H": 18.628894992991583, + "H_err": 0.01716803531553185, + "phase_parameter_1": -0.253, + "phase_parameter_1_err": 0.00014324314956736168, + "phase_parameter_2": np.nan, + "phase_parameter_2_err": np.nan, } -test_object.populate_phase_parameters("u", **u_model_1) -test_object.populate_phase_parameters("u", **u_model_2) -test_object.populate_phase_parameters("g", **g_model_1) -test_object.populate_phase_parameters("r", **r_model_2) +# u_model_1 = { +# "model_name": "model_1", +# "phaseAngle_min": 11.0, +# "phaseAngle_range": 12.0, +# "nobs": 13, +# "arc": 14.0, +# "H": 15.0, +# "H_err": 16.0, +# "phase_parameter_1": 17.0, +# "phase_parameter_1_err": 18.0, +# } + +# u_model_2 = { +# "model_name": "model_2", +# "H": 25.0, +# "H_err": 26.0, +# "phase_parameter_1": 27.0, +# "phase_parameter_1_err": 28.0, +# "phase_parameter_2": 29.0, +# "phase_parameter_2_err": 30.0, +# } + +# g_model_1 = { +# "model_name": "model_1", +# "phaseAngle_min": 31.0, +# "phaseAngle_range": 32.0, +# "nobs": 33, +# "arc": 34.0, +# "H": 35.0, +# "H_err": 36.0, +# "phase_parameter_1": 37.0, +# "phase_parameter_1_err": 38.0, +# } + +# r_model_2 = { +# "model_name": "model_2", +# "phaseAngle_min": 41.0, +# "phaseAngle_range": 42.0, +# "nobs": 43, +# "arc": 44.0, +# "H": 45.0, +# "H_err": 46.0, +# "phase_parameter_1": 47.0, +# "phase_parameter_1_err": 48.0, +# "phase_parameter_2": 49.0, +# "phase_parameter_2_err": 50.0, +# } + +_g_model_1 = g_model_1.copy() +del _g_model_1["filter_name"] +_g_model_2 = g_model_2.copy() +del _g_model_2["filter_name"] +_r_model_1 = r_model_1.copy() +del _r_model_1["filter_name"] +_i_model_2 = i_model_2.copy() +del _i_model_2["filter_name"] + +test_object.populate_phase_parameters("g", **_g_model_1) +test_object.populate_phase_parameters("g", **_g_model_2) +test_object.populate_phase_parameters("r", **_r_model_1) +test_object.populate_phase_parameters("i", **_i_model_2) def test_populate_phase_parameters(): # test to make sure the object is correctly populated - assert test_object.filter_list == ["u", "g", "r"] - - assert test_object.filter_dependent_values[0].model_list == ["model_1", "model_2"] - assert test_object.filter_dependent_values[1].model_list == ["model_1"] - assert test_object.filter_dependent_values[2].model_list == ["model_2"] - - assert_array_equal([a.phaseAngle_min for a in test_object.filter_dependent_values], [11.0, 31.0, 41.0]) - assert_array_equal([a.phaseAngle_range for a in test_object.filter_dependent_values], [12.0, 32.0, 42.0]) - assert_array_equal([a.nobs for a in test_object.filter_dependent_values], [13, 33, 43]) - assert_array_equal([a.arc for a in test_object.filter_dependent_values], [14.0, 34.0, 44.0]) - - assert test_object.filter_dependent_values[0].model_dependent_values[0].__dict__ == { - "filter_name": "u", - "model_name": "model_1", - "H": 15.0, - "H_err": 16.0, - "phase_parameter_1": 17.0, - "phase_parameter_1_err": 18.0, - "phase_parameter_2": np.nan, - "phase_parameter_2_err": np.nan, - } - - assert test_object.filter_dependent_values[0].model_dependent_values[1].__dict__ == { - "filter_name": "u", - "model_name": "model_2", - "H": 25.0, - "H_err": 26.0, - "phase_parameter_1": 27.0, - "phase_parameter_1_err": 28.0, - "phase_parameter_2": 29.0, - "phase_parameter_2_err": 30.0, - } - - assert test_object.filter_dependent_values[1].model_dependent_values[0].__dict__ == { - "filter_name": "g", - "model_name": "model_1", - "H": 35.0, - "H_err": 36.0, - "phase_parameter_1": 37.0, - "phase_parameter_1_err": 38.0, - "phase_parameter_2": np.nan, - "phase_parameter_2_err": np.nan, - } - - assert test_object.filter_dependent_values[2].model_dependent_values[0].__dict__ == { - "filter_name": "r", - "model_name": "model_2", - "H": 45.0, - "H_err": 46.0, - "phase_parameter_1": 47.0, - "phase_parameter_1_err": 48.0, - "phase_parameter_2": 49.0, - "phase_parameter_2_err": 50.0, - } + assert test_object.filter_list == ["g", "r", "i"] + + assert test_object.filter_dependent_values[0].model_list == [model_1, model_2] + assert test_object.filter_dependent_values[1].model_list == [model_1] + assert test_object.filter_dependent_values[2].model_list == [model_2] + + assert_almost_equal( + [a.phaseAngle_min for a in test_object.filter_dependent_values], + [g_model_1["phaseAngle_min"], r_model_1["phaseAngle_min"], i_model_1["phaseAngle_min"]], + ) + assert_array_equal( + [a.phaseAngle_range for a in test_object.filter_dependent_values], + [g_model_1["phaseAngle_range"], r_model_1["phaseAngle_range"], i_model_1["phaseAngle_range"]], + ) + assert_array_equal( + [a.nobs for a in test_object.filter_dependent_values], + [g_model_1["nobs"], r_model_1["nobs"], i_model_1["nobs"]], + ) + assert_array_equal( + [a.arc for a in test_object.filter_dependent_values], + [g_model_1["arc"], r_model_1["arc"], i_model_1["arc"]], + ) + + # assert test_object.filter_dependent_values[0].model_dependent_values[0].__dict__ == g_model_1 + # { + # "filter_name": "u", + # "model_name": "model_1", + # "H": 15.0, + # "H_err": 16.0, + # "phase_parameter_1": 17.0, + # "phase_parameter_1_err": 18.0, + # "phase_parameter_2": np.nan, + # "phase_parameter_2_err": np.nan, + # } + + # assert test_object.filter_dependent_values[0].model_dependent_values[1].__dict__ == g_model_2 + # { + # "filter_name": "u", + # "model_name": "model_2", + # "H": 25.0, + # "H_err": 26.0, + # "phase_parameter_1": 27.0, + # "phase_parameter_1_err": 28.0, + # "phase_parameter_2": 29.0, + # "phase_parameter_2_err": 30.0, + # } + + # assert test_object.filter_dependent_values[1].model_dependent_values[0].__dict__ == r_model_1 + # { + # "filter_name": "g", + # "model_name": "model_1", + # "H": 35.0, + # "H_err": 36.0, + # "phase_parameter_1": 37.0, + # "phase_parameter_1_err": 38.0, + # "phase_parameter_2": np.nan, + # "phase_parameter_2_err": np.nan, + # } + + # assert test_object.filter_dependent_values[2].model_dependent_values[0].__dict__ == i_model_2 + # { + # "filter_name": "r", + # "model_name": "model_2", + # "H": 45.0, + # "H_err": 46.0, + # "phase_parameter_1": 47.0, + # "phase_parameter_1_err": 48.0, + # "phase_parameter_2": 49.0, + # "phase_parameter_2_err": 50.0, + # } + + test_dict_list = [ + test_object.get_phase_parameters_in_filter("g", model_1).__dict__, + test_object.get_phase_parameters_in_filter("g", model_2).__dict__, + test_object.get_phase_parameters_in_filter("r", model_1).__dict__, + test_object.get_phase_parameters_in_filter("i", model_2).__dict__, + ] + expect_dict_list = [g_model_1, g_model_2, r_model_1, i_model_2] + + for test_dict, expect_dict in zip(test_dict_list, expect_dict_list): + for x in test_dict.keys(): + # print(x) + test_val = test_dict[x] + expect_val = expect_dict[x] + if type(expect_val) == str: + assert test_val == expect_val + else: + assert_almost_equal(test_val, expect_val) # check to make sure model-dependent parameter is correctly updated (then return it to previous) - test_object.populate_phase_parameters("u", model_name="model_1", H=99.0) - assert test_object.filter_dependent_values[0].model_dependent_values[0].H == 99.0 - test_object.populate_phase_parameters("u", model_name="model_1", H=15.0) + test_object.populate_phase_parameters("g", model_name=model_1, H=99.0) + assert test_object.get_phase_parameters_in_filter("g", model_1).H == 99.0 + test_object.populate_phase_parameters("g", model_name=model_1, H=g_model_1["H"]) # check to make sure filter-dependent parameter is correctly updated (then return it to previous) - test_object.populate_phase_parameters("u", nobs=99) - assert test_object.filter_dependent_values[0].nobs == 99 - test_object.populate_phase_parameters("u", nobs=13) + # print(test_object.__dict__) + test_object.populate_phase_parameters("r", model_name=model_1, nobs=99) + # print(test_object.__dict__) + # print(r_model_1["nobs"]) + # print(test_object.get_phase_parameters_in_filter("r",model_1).__dict__) + # print(test_object.get_phase_parameters_in_filter("r",model_1).nobs) + # assert test_object.filter_dependent_values[0].nobs == 99 # TODO: the indexing in this line does not access the correct filter + assert test_object.get_phase_parameters_in_filter("r", model_1).nobs == 99 + # print(r_model_1["nobs"]) + test_object.populate_phase_parameters("r", model_name=model_1, nobs=r_model_1["nobs"]) + # print(test_object.get_phase_parameters_in_filter("r",model_1).nobs) # testing to make sure the correct error messages trigger + test_filt = "y" with pytest.raises(ValueError) as error_info_1: - test_object.populate_phase_parameters("y") + test_object.populate_phase_parameters(test_filt) - assert error_info_1.value.args[0] == "Filter y does not exist in AdlerData.filter_list." + assert error_info_1.value.args[0] == "Filter {} does not exist in AdlerData.filter_list.".format( + test_filt + ) with pytest.raises(NameError) as error_info_2: - test_object.populate_phase_parameters("u", H=4.0) + test_object.populate_phase_parameters("r", H=4.0) assert error_info_2.value.args[0] == "No model name given. Cannot update model-specific phase parameters." def test_get_phase_parameters_in_filter(): - assert test_object.get_phase_parameters_in_filter("u", "model_1").__dict__ == { - "filter_name": "u", - "phaseAngle_min": 11.0, - "phaseAngle_range": 12.0, - "nobs": 13, - "arc": 14.0, - "model_name": "model_1", - "H": 15.0, - "H_err": 16.0, - "phase_parameter_1": 17.0, - "phase_parameter_1_err": 18.0, - "phase_parameter_2": np.nan, - "phase_parameter_2_err": np.nan, - } - - assert test_object.get_phase_parameters_in_filter("u", "model_2").__dict__ == { - "filter_name": "u", - "phaseAngle_min": 11.0, - "phaseAngle_range": 12.0, - "nobs": 13, - "arc": 14.0, - "model_name": "model_2", - "H": 25.0, - "H_err": 26.0, - "phase_parameter_1": 27.0, - "phase_parameter_1_err": 28.0, - "phase_parameter_2": 29.0, - "phase_parameter_2_err": 30.0, - } - - assert test_object.get_phase_parameters_in_filter("g", "model_1").__dict__ == { - "filter_name": "g", - "phaseAngle_min": 31.0, - "phaseAngle_range": 32.0, - "nobs": 33, - "arc": 34.0, - "model_name": "model_1", - "H": 35.0, - "H_err": 36.0, - "phase_parameter_1": 37.0, - "phase_parameter_1_err": 38.0, - "phase_parameter_2": np.nan, - "phase_parameter_2_err": np.nan, - } - - assert test_object.get_phase_parameters_in_filter("r", "model_2").__dict__ == { - "filter_name": "r", - "phaseAngle_min": 41.0, - "phaseAngle_range": 42.0, - "nobs": 43, - "arc": 44.0, - "model_name": "model_2", - "H": 45.0, - "H_err": 46.0, - "phase_parameter_1": 47.0, - "phase_parameter_1_err": 48.0, - "phase_parameter_2": 49.0, - "phase_parameter_2_err": 50.0, - } + # assert test_object.get_phase_parameters_in_filter("g", model_1).__dict__ == g_model_1 + # { + # "filter_name": "u", + # "phaseAngle_min": 11.0, + # "phaseAngle_range": 12.0, + # "nobs": 13, + # "arc": 14.0, + # "model_name": "model_1", + # "H": 15.0, + # "H_err": 16.0, + # "phase_parameter_1": 17.0, + # "phase_parameter_1_err": 18.0, + # "phase_parameter_2": np.nan, + # "phase_parameter_2_err": np.nan, + # } + + # assert test_object.get_phase_parameters_in_filter("g", model_2).__dict__ == g_model_2 + # { + # "filter_name": "u", + # "phaseAngle_min": 11.0, + # "phaseAngle_range": 12.0, + # "nobs": 13, + # "arc": 14.0, + # "model_name": "model_2", + # "H": 25.0, + # "H_err": 26.0, + # "phase_parameter_1": 27.0, + # "phase_parameter_1_err": 28.0, + # "phase_parameter_2": 29.0, + # "phase_parameter_2_err": 30.0, + # } + + # assert test_object.get_phase_parameters_in_filter("r", model_1).__dict__ == r_model_1 + # { + # "filter_name": "g", + # "phaseAngle_min": 31.0, + # "phaseAngle_range": 32.0, + # "nobs": 33, + # "arc": 34.0, + # "model_name": "model_1", + # "H": 35.0, + # "H_err": 36.0, + # "phase_parameter_1": 37.0, + # "phase_parameter_1_err": 38.0, + # "phase_parameter_2": np.nan, + # "phase_parameter_2_err": np.nan, + # } + + # assert test_object.get_phase_parameters_in_filter("i", model_2).__dict__ == i_model_2 + # { + # "filter_name": "r", + # "phaseAngle_min": 41.0, + # "phaseAngle_range": 42.0, + # "nobs": 43, + # "arc": 44.0, + # "model_name": "model_2", + # "H": 45.0, + # "H_err": 46.0, + # "phase_parameter_1": 47.0, + # "phase_parameter_1_err": 48.0, + # "phase_parameter_2": 49.0, + # "phase_parameter_2_err": 50.0, + # } + + test_dict_list = [ + test_object.get_phase_parameters_in_filter("g", model_1).__dict__, + test_object.get_phase_parameters_in_filter("g", model_2).__dict__, + test_object.get_phase_parameters_in_filter("r", model_1).__dict__, + test_object.get_phase_parameters_in_filter("i", model_2).__dict__, + ] + expect_dict_list = [g_model_1, g_model_2, r_model_1, i_model_2] + + for test_dict, expect_dict in zip(test_dict_list, expect_dict_list): + for x in test_dict.keys(): + # print(x) + test_val = test_dict[x] + expect_val = expect_dict[x] + if type(expect_val) == str: + assert test_val == expect_val + else: + assert_almost_equal(test_val, expect_val) # checking the error messages + test_filt = "f" with pytest.raises(ValueError) as error_info_1: - error_dict = test_object.get_phase_parameters_in_filter("f", model_name="model_2") + error_dict = test_object.get_phase_parameters_in_filter(test_filt, model_name=model_2) - assert error_info_1.value.args[0] == "Filter f does not exist in AdlerData.filter_list." + assert error_info_1.value.args[0] == "Filter {} does not exist in AdlerData.filter_list.".format( + test_filt + ) + test_filt = "r" with pytest.raises(ValueError) as error_info_2: - error_dict_2 = test_object.get_phase_parameters_in_filter("r", model_name="model_1") + error_dict_2 = test_object.get_phase_parameters_in_filter(test_filt, model_name=model_2) - assert error_info_2.value.args[0] == "Model model_1 does not exist for filter r in AdlerData.model_lists." + print(error_info_2.value.args[0]) + print("Model {} does not exist for filter {} in AdlerData.model_lists.".format(model_2, test_filt)) + assert error_info_2.value.args[ + 0 + ] == "Model {} does not exist for filter {} in AdlerData.model_lists.".format(model_2, test_filt) # here the capsys fixture captures any output to the terminal @@ -226,13 +401,17 @@ def test_print_data(capsys): # get what was printed to the terminal captured = capsys.readouterr() - expected = "Filter: u\nPhase angle minimum: 11.0\nPhase angle range: 12.0\nNumber of observations: 13\nArc: 14.0\nModel: model_1.\n\tH: 15.0\n\tH error: 16.0\n\tPhase parameter 1: 17.0\n\tPhase parameter 1 error: 18.0\n\tPhase parameter 2: nan\n\tPhase parameter 2 error: nan\nModel: model_2.\n\tH: 25.0\n\tH error: 26.0\n\tPhase parameter 1: 27.0\n\tPhase parameter 1 error: 28.0\n\tPhase parameter 2: 29.0\n\tPhase parameter 2 error: 30.0\n\n\nFilter: g\nPhase angle minimum: 31.0\nPhase angle range: 32.0\nNumber of observations: 33\nArc: 34.0\nModel: model_1.\n\tH: 35.0\n\tH error: 36.0\n\tPhase parameter 1: 37.0\n\tPhase parameter 1 error: 38.0\n\tPhase parameter 2: nan\n\tPhase parameter 2 error: nan\n\n\nFilter: r\nPhase angle minimum: 41.0\nPhase angle range: 42.0\nNumber of observations: 43\nArc: 44.0\nModel: model_2.\n\tH: 45.0\n\tH error: 46.0\n\tPhase parameter 1: 47.0\n\tPhase parameter 1 error: 48.0\n\tPhase parameter 2: 49.0\n\tPhase parameter 2 error: 50.0\n\n\n" + expected = "Filter: g\nPhase angle minimum: 27.426668167114258\nPhase angle range: 36.17388343811035\nNumber of observations: 9\nArc: 3306.0079999999944\nModel: HG12_Pen16.\n\tH: 20.719465178426468\n\tH error: 0.018444134023977106\n\tPhase parameter 1: 0.62\n\tPhase parameter 1 error: nan\n\tPhase parameter 2: nan\n\tPhase parameter 2 error: nan\nModel: HG.\n\tH: 20.72954332871528\n\tH error: 0.018444134023977106\n\tPhase parameter 1: 0.15\n\tPhase parameter 1 error: nan\n\tPhase parameter 2: nan\n\tPhase parameter 2 error: nan\n\n\nFilter: r\nPhase angle minimum: 2.553332567214966\nPhase angle range: 124.23803400993347\nNumber of observations: 38\nArc: 3338.0655999999944\nModel: HG12_Pen16.\n\tH: 19.92863542616601\n\tH error: 0.018525355171274356\n\tPhase parameter 1: 1.0\n\tPhase parameter 1 error: 0.05300829494059732\n\tPhase parameter 2: nan\n\tPhase parameter 2 error: nan\n\n\nFilter: i\nPhase angle minimum: 10.0595064163208\nPhase angle range: 101.74150371551514\nNumber of observations: 32\nArc: 3342.0585599999977\nModel: HG.\n\tH: 18.628894992991583\n\tH error: 0.01716803531553185\n\tPhase parameter 1: -0.253\n\tPhase parameter 1 error: 0.00014324314956736168\n\tPhase parameter 2: nan\n\tPhase parameter 2 error: nan\n\n\n" + # print(captured) + # print(expected) assert captured.out == expected def test_write_row_to_database(tmp_path): db_location = os.path.join(tmp_path, "test_AdlerData_database.db") + # print(db_location) + # print(test_object.__dict__) test_object.write_row_to_database(db_location) con = sqlite3.connect(db_location) @@ -240,17 +419,33 @@ def test_write_row_to_database(tmp_path): con.close() expected_data_filepath = get_test_data_filepath("test_SQL_database_table.csv") - expected_data = pd.read_csv(expected_data_filepath) + expected_data = pd.read_csv(expected_data_filepath, index_col=0) # we don't expect the timestamp column to be the same, obviously - expected_data = expected_data.drop(columns="timestamp") - written_data = written_data.drop(columns="timestamp") + drop_cols = [x for x in written_data if (x == "timestamp") or ("modelFitMjd" in x)] + expected_data = expected_data.drop(columns=drop_cols) + written_data = written_data.drop(columns=drop_cols) # note that because I'm using Pandas there's some small dtype and np.nan/None stuff to clear up # but this makes for a quick streamlined test anyway - expected_data = expected_data.replace({np.nan: None}) + # expected_data = expected_data.replace({np.nan: None}) # TODO: fix weirdness between nan and None? + written_data = written_data.replace({None: np.nan}) # TODO: fix weirdness between nan and None? expected_data = expected_data.astype({"ssObjectId": str}) - pd.testing.assert_frame_equal(expected_data, written_data, check_dtype=False) + # print(expected_data.dtypes) + # print(written_data.dtypes) + # print(expected_data.iloc[0].to_dict()) + # print(written_data.iloc[0].to_dict()) + for x in written_data.columns: + # print(x) + # print(expected_data.iloc[0][x],written_data.iloc[0][x]) + # print(type(expected_data.iloc[0][x]),type(written_data.iloc[0][x])) + # print(expected_data.iloc[0][x]==written_data.iloc[0][x]) + if type(written_data.iloc[0][x]) == str: + assert expected_data.iloc[0][x] == written_data.iloc[0][x] + else: + assert_almost_equal(expected_data.iloc[0][x], written_data.iloc[0][x]) + # pd.testing.assert_frame_equal(expected_data, written_data, check_dtype=False) + # assert expected_data.iloc[0].to_dict() == pytest.approx(written_data.iloc[0].to_dict()) def test_read_row_from_database(): @@ -259,28 +454,50 @@ def test_read_row_from_database(): db_location = get_test_data_filepath("test_AdlerData_database.db") - new_object = AdlerData("8268570668335894776", ["u", "g", "r"]) + new_object = AdlerData("8268570668335894776", ["g", "r", "i"]) new_object.populate_from_database(db_location) - assert new_object.__dict__ == test_object.__dict__ + # print(new_object.__dict__) + # print(test_object.__dict__) + # assert new_object.__dict__ == test_object.__dict__ + + for filt, model in zip(["g", "g", "r", "i"], [model_1, model_2, model_1, model_2]): + # print(filt,model) + new = new_object.get_phase_parameters_in_filter(filt, model) + # print(new.__dict__) + test = test_object.get_phase_parameters_in_filter(filt, model) + # print(test.__dict__) + # assert new.__dict__ == test.__dict__ + # assert new.__dict__ == pytest.approx(test.__dict__) + + for x in new.__dict__.keys(): + print(x) + new_val = new.__dict__[x] + test_val = test.__dict__[x] + print(new_val, test_val) + if type(new_val) == str: + assert new_val == test_val + else: + assert_almost_equal(new_val, test_val) with pytest.raises(ValueError) as error_info_1: - empty_data = AdlerData("pretend_object", ["u", "g", "r"]) + empty_data = AdlerData("pretend_object", ["g", "r", "i"]) empty_data.populate_from_database(db_location) assert error_info_1.value.args[0] == "No data found in this database for the supplied ssObjectId." with pytest.raises(ValueError) as error_info_2: - bad_filter = AdlerData("8268570668335894776", ["u", "g", "h"]) + bad_filter = AdlerData("8268570668335894776", ["g", "r", "h"]) bad_filter.populate_from_database(db_location) - assert ( - error_info_2.value.args[0] - == "Data does not exist for some of the requested filters in this database. Filters in database for this object: ['u', 'g', 'r']" + assert error_info_2.value.args[ + 0 + ] == "Data does not exist for some of the requested filters in this database. Filters in database for this object: {}".format( + test_object.filter_list ) with pytest.raises(ValueError) as error_info_3: - bad_filter = AdlerData("8268570668335894776", ["u", "g", "h"]) + bad_filter = AdlerData("8268570668335894776", ["g", "r", "h"]) bad_filter.populate_from_database("./dummy_location.db") assert error_info_3.value.args[0] == "Database cannot be found at given filepath." diff --git a/tests/adler/objectdata/test_AdlerPlanetoid.py b/tests/adler/objectdata/test_AdlerPlanetoid.py index 4716b64..1bdbaa5 100644 --- a/tests/adler/objectdata/test_AdlerPlanetoid.py +++ b/tests/adler/objectdata/test_AdlerPlanetoid.py @@ -6,7 +6,7 @@ from adler.objectdata.AdlerPlanetoid import AdlerPlanetoid -ssoid = 8268570668335894776 +ssoid = "8268570668335894776" test_db_path = get_test_data_filepath("testing_database.db") @@ -160,25 +160,43 @@ def test_failed_SQL_queries(): def test_attach_previous_adlerdata(): test_planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, test_db_path, filter_list=["g", "r"]) + # TODO: this setup is currently a bit dodgy. Because AdlerData write_row_to_database appends a new line to the db the most recent model for a given object may be nan in that row + # as such this test depends on the order/number of times the adler commands have been run to make it + + # the test database can be recreated by running the Adler commands in the tests/data dir: + # adler -s 8268570668335894776 -i testing_database.db -n test_AdlerData_database.db + # adler -s 8268570668335894776 -i testing_database.db -n test_AdlerData_database.db -m HG db_location = get_test_data_filepath("test_AdlerData_database.db") + print(db_location) test_planetoid.attach_previous_adler_data(db_location) - test_output = test_planetoid.PreviousAdlerData.get_phase_parameters_in_filter("g", "model_1") + test_output = test_planetoid.PreviousAdlerData.get_phase_parameters_in_filter("r", "HG12_Pen16") + print(test_output.__dict__) expected_output = { - "filter_name": "g", - "phaseAngle_min": 31.0, - "phaseAngle_range": 32.0, - "nobs": 33, - "arc": 34.0, - "model_name": "model_1", - "H": 35.0, - "H_err": 36.0, - "phase_parameter_1": 37.0, - "phase_parameter_1_err": 38.0, + "filter_name": "r", + "phaseAngle_min": 2.553332567214966, + "phaseAngle_range": 124.23803400993347, + "nobs": 38, + "arc": 3338.0655999999944, + "model_name": "HG12_Pen16", + "H": 19.92863542616601, + "H_err": 0.018525355171274356, + "phase_parameter_1": 1.0, + "phase_parameter_1_err": 0.05300829494059732, "phase_parameter_2": np.nan, "phase_parameter_2_err": np.nan, } - - assert test_output.__dict__ == expected_output + print(expected_output) + + # assert test_output.__dict__ == pytest.approx(expected_output) # TODO: why isn't this working now? + + for x in test_output.__dict__.keys(): + print(x) + test_val = test_output.__dict__[x] + expect_val = expected_output[x] + if type(expect_val) == str: + assert test_val == expect_val + else: + assert_almost_equal(test_val, expect_val) diff --git a/tests/adler/science/test_Colour.py b/tests/adler/science/test_Colour.py index 307eb68..2240762 100644 --- a/tests/adler/science/test_Colour.py +++ b/tests/adler/science/test_Colour.py @@ -53,11 +53,18 @@ H = sso.H G12 = sso.G12 - pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name="HG12_Pen16") + pc = PhaseCurve( + H=H, + # H=H * u.mag, + phase_parameter_1=G12, + model_name="HG12_Pen16", + ) pc_fit = pc.FitModel( - np.array(getattr(obs, "phaseAngle")) * u.deg, - np.array(getattr(obs, "reduced_mag")) * u.mag, + np.radians(np.array(getattr(obs, "phaseAngle"))), + np.array(getattr(obs, "reduced_mag")), + # np.array(getattr(obs, "phaseAngle")) * u.deg, + # np.array(getattr(obs, "reduced_mag")) * u.mag, ) pc = pc.InitModelSbpy(pc_fit) @@ -110,7 +117,17 @@ def test_col_obs_ref( ) # merge with observation data (avoiding duplicating x_col) # compare results df to stored df + # NB: see also pytest.approx for comparing dictionaries for x in list(df_col): print(x) assert_array_almost_equal(np.array(df_col[x]), np.array(df_ref[x])) # TODO: diasourceId failing on ubuntu tests, due to float? need to save as str/int? + + +# test_col_obs_ref( +# planetoid=planetoid, +# adler_data=adler_data, +# column_dict=column_dict, +# N_ref_list=[1, 3, 5], +# df_ref_list=[df_N_ref_1, df_N_ref_3, df_N_ref_5], +# ) diff --git a/tests/adler/utilities/test_AdlerCLIArguments.py b/tests/adler/utilities/test_AdlerCLIArguments.py index 1a49be3..ecf4368 100644 --- a/tests/adler/utilities/test_AdlerCLIArguments.py +++ b/tests/adler/utilities/test_AdlerCLIArguments.py @@ -16,6 +16,8 @@ def __init__( outpath, db_name, sql_filename, + plot_show, + phase_model, ): self.ssObjectId = ssObjectId self.ssObjectId_list = ssObjectId_list @@ -25,6 +27,8 @@ def __init__( self.outpath = outpath self.db_name = db_name self.sql_filename = sql_filename + self.plot_show = plot_show + self.phase_model = phase_model def test_AdlerCLIArguments_population(): @@ -38,6 +42,8 @@ def test_AdlerCLIArguments_population(): "outpath": "./", "db_name": "output", "sql_filename": None, + "plot_show": False, + "phase_model": "HG12_Pen16", } good_arguments = args(**good_input_dict) good_arguments_object = AdlerCLIArguments(good_arguments) @@ -64,7 +70,16 @@ def test_AdlerCLIArguments_population(): def test_AdlerCLIArguments_badSSOID(): # test that a bad ssObjectId triggers the right error bad_ssoid_arguments = args( - "hello!", None, ["g", "r", "i"], ["g-r", "r-i"], [60000.0, 67300.0], "./", "output", "None" + "hello!", + None, + ["g", "r", "i"], + ["g-r", "r-i"], + [60000.0, 67300.0], + "./", + "output", + "None", + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_ssoid_error: @@ -79,7 +94,16 @@ def test_AdlerCLIArguments_badSSOID(): def test_AdlerCLIArguments_badfilters(): # test that non-LSST or unexpected filters trigger the right error bad_filter_arguments = args( - "666", None, ["g", "r", "i", "m"], ["g-r", "r-i"], [60000.0, 67300.0], "./", "output", None + "666", + None, + ["g", "r", "i", "m"], + ["g-r", "r-i"], + [60000.0, 67300.0], + "./", + "output", + None, + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_filter_error: @@ -91,7 +115,16 @@ def test_AdlerCLIArguments_badfilters(): ) bad_filter_arguments_2 = args( - "666", None, ["pony"], ["g-r", "r-i"], [60000.0, 67300.0], "./", "output", None + "666", + None, + ["pony"], + ["g-r", "r-i"], + [60000.0, 67300.0], + "./", + "output", + None, + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_filter_error_2: @@ -108,7 +141,16 @@ def test_AdlerCLIArguments_badfilters(): # the colours are not in the right format bad_colour_arguments_1 = args( - "666", None, ["g", "r", "i"], ["g - r", "r x i"], [60000.0, 67300.0], "./", "output", None + "666", + None, + ["g", "r", "i"], + ["g - r", "r x i"], + [60000.0, 67300.0], + "./", + "output", + None, + False, + "HG12_Pen16", ) err_msg1 = "Unexpected filters found in --colour_list command-line argument. --colour_list must contain LSST filters in the format 'filter2-filter1'." @@ -119,7 +161,16 @@ def test_AdlerCLIArguments_badfilters(): # colours are requested in filters that are not available bad_colour_arguments_2 = args( - "666", None, ["g", "r", "i"], ["g-r", "r-j"], [60000.0, 67300.0], "./", "output", None + "666", + None, + ["g", "r", "i"], + ["g-r", "r-j"], + [60000.0, 67300.0], + "./", + "output", + None, + False, + "HG12_Pen16", ) err_msg2 = err_msg1 @@ -130,7 +181,16 @@ def test_AdlerCLIArguments_badfilters(): # colours are requested in filters that are not available bad_colour_arguments_3 = args( - "666", None, ["g", "r"], ["g-r", "r-i"], [60000.0, 67300.0], "./", "output", None + "666", + None, + ["g", "r"], + ["g-r", "r-i"], + [60000.0, 67300.0], + "./", + "output", + None, + False, + "HG12_Pen16", ) err_msg3 = "The filters required to calculate the colours have not been requested in --filter-list" @@ -143,7 +203,16 @@ def test_AdlerCLIArguments_badfilters(): def test_AdlerCLIArguments_baddates(): # test that overly-large dates trigger the right error big_date_arguments = args( - "666", None, ["g", "r", "i"], ["g-r", "r-i"], [260000.0, 267300.0], "./", "output", None + "666", + None, + ["g", "r", "i"], + ["g-r", "r-i"], + [260000.0, 267300.0], + "./", + "output", + None, + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as big_date_error: @@ -156,7 +225,16 @@ def test_AdlerCLIArguments_baddates(): # test that unexpected date values trigger the right error bad_date_arguments = args( - "666", None, ["g", "r", "i"], ["g-r", "r-i"], [60000.0, "cheese"], "./", "output", None + "666", + None, + ["g", "r", "i"], + ["g-r", "r-i"], + [60000.0, "cheese"], + "./", + "output", + None, + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_date_error: @@ -178,6 +256,8 @@ def test_AdlerCLIArguments_badoutput(): "./definitely_fake_folder/", "output", None, + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_output_error: @@ -199,6 +279,8 @@ def test_AdlerCLIArguments_badlist(): "./", "output", "./", + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_list_error: @@ -220,6 +302,8 @@ def test_AdlerCLIArguments_badsql(): "./", "output", "./dummy_database.db", + False, + "HG12_Pen16", ) with pytest.raises(ValueError) as bad_sql_error: @@ -229,3 +313,6 @@ def test_AdlerCLIArguments_badsql(): bad_sql_error.value.args[0] == "The file supplied for the command-line argument --sql_filename cannot be found." ) + + +# TODO: test plot_show somehow, and phase_model diff --git a/tests/adler/utilities/test_science_utilities.py b/tests/adler/utilities/test_science_utilities.py index b8a0588..19eef3c 100644 --- a/tests/adler/utilities/test_science_utilities.py +++ b/tests/adler/utilities/test_science_utilities.py @@ -213,7 +213,7 @@ def test_zero_func(): def test_sigma_clip(): - outliers = sci_utils.sigma_clip(y_res, {"maxiters": 1, "cenfunc": sci_utils.zero_func}) + outliers = sci_utils.sigma_clip(y_res, sigma=3, maxiters=1, cenfunc=sci_utils.zero_func) assert_array_equal(outliers, outliers_y) diff --git a/tests/data/adler_run_test_data.py b/tests/data/adler_run_test_data.py new file mode 100644 index 0000000..4ddbf1e --- /dev/null +++ b/tests/data/adler_run_test_data.py @@ -0,0 +1,247 @@ +import logging +import argparse +import astropy.units as u +from astropy.time import Time +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import os +import sqlite3 + +from adler.objectdata.AdlerPlanetoid import AdlerPlanetoid +from adler.objectdata.AdlerData import AdlerData +from adler.science.PhaseCurve import PhaseCurve, ADLER_SBPY_DICT +from adler.science.Colour import col_obs_ref +from adler.utilities.AdlerCLIArguments import AdlerCLIArguments +from adler.utilities.adler_logging import setup_adler_logging +from adler.utilities.readin_utilities import read_in_SSObjectID_file +from adler.utilities.plotting_utilities import plot_errorbar +import adler.utilities.science_utilities as sci_utils + +""" +Copy of adler_run.oy which can be used to make the AdlerData testing database +""" + +ssObjectId = "8268570668335894776" +filter_list = ["g", "r", "i"] +colour_list = None +date_range = [60000.0, 67300.0] +outpath = "." +db_name = "test_AdlerData_database.db" +sql_filename = "testing_database.db" +plot_show = False +# phase_model="HG12_Pen16" +phase_model_list = ["HG12_Pen16", "HG"] + +# adler parameters +N_pc_fit = 10 # minimum number of data points to fit phase curve +diff_cut = 1.0 # magnitude difference used to identify outliers +obs_cols = ["diaSourceId", "midPointMjdTai", "outlier"] # observation columns to use + + +# Define colour parameters +# set number of reference observations to use for colour estimate +N_ref = 5 + +ssObjectId_list = [ssObjectId] + +# consider each ssObjectId in the list separately +for i, ssObjectId in enumerate(ssObjectId_list): + print( + "Query data in the range: {} <= date <= {}".format(date_range[0], date_range[1]) + ) # the adler planetoid date_range is used in an SQL BETWEEN statement which is inclusive + + # load ssObjectId data + if sql_filename: # load from a local database, if provided + msg = "query sql database {}".format(sql_filename) + print(msg) + planetoid = AdlerPlanetoid.construct_from_SQL( + ssObjectId, + filter_list=filter_list, + date_range=date_range, + sql_filename=sql_filename, + ) + else: # otherwise load from the Rubin Science Platform + msg = "query RSP" + print(msg) + planetoid = AdlerPlanetoid.construct_from_RSP(ssObjectId, filter_list, date_range) + + # TODO: Here we would load the AdlerData object from our data tables + adler_data = AdlerData(ssObjectId, planetoid.filter_list) + print(adler_data.__dict__) + + # now let's do some phase curves! + + # make an empty figure + fig = plot_errorbar(planetoid, filt_list=[]) + + # do each model + for phase_model in phase_model_list: + + print(phase_model) + # get the name of the phase parameter + # TODO: make an option for two parameter HG1G2 + phase_parameter_1 = ADLER_SBPY_DICT[phase_model]["phase_parameter_1"] + print(phase_model, phase_parameter_1) + + # set default phase parameter values + # TODO: make an option for two parameter HG1G2 + if phase_model == "HG": + phase_param_1_default = 0.15 + else: + phase_param_1_default = 0.62 + + # operate on each filter in turn + for filt in filter_list: + + # get the filter SSObject metadata + try: + sso = planetoid.SSObject_in_filter(filt) + except: + continue + + # get the observations + obs = planetoid.observations_in_filter(filt) + df_obs = pd.DataFrame(obs.__dict__) + df_obs["outlier"] = [False] * len(df_obs) + + # Determine the reference phase curve model + # TODO: We would load the best phase curve model available in AdlerData here + + # initial simple phase curve filter model with fixed phase_parameter + # use the ssObject value for H as initial guess, this is in HG12_Pen16 + # TODO: use the ssObject value for phase parameter as initial guess? + pc = PhaseCurve( + # H=sso.H * u.mag, + H=sso.H, + phase_parameter_1=phase_param_1_default, + model_name=phase_model, + ) + + # only fit phase_parameter when sufficient data is available + if len(df_obs) < N_pc_fit: + msg = "Do not fit {}, use {}={:.2f}".format( + phase_parameter_1, phase_parameter_1, pc.phase_parameter_1 + ) + # pc.model_function.G12.fixed = True + getattr(pc.model_function, phase_parameter_1).fixed = True + else: + # pc.model_function.G12.fixed = False + getattr(pc.model_function, phase_parameter_1).fixed = False + + # do a phase model fit to the past data + pc_fit = pc.FitModel( + # np.array(df_obs["phaseAngle"]) * u.deg, + # np.array(df_obs["reduced_mag"]) * u.mag, + # np.array(df_obs["magErr"]) * u.mag, + np.radians(np.array(df_obs["phaseAngle"])), + np.array(df_obs["reduced_mag"]), + np.array(df_obs["magErr"]), + ) + pc_fit = pc.InitModelSbpy(pc_fit) + + # Store the fitted values, and metadata, in an AdlerData object + ad_params = pc_fit.__dict__ + ad_params["phaseAngle_min"] = np.amin(df_obs["phaseAngle"]) # * u.deg + ad_params["phaseAngle_range"] = np.ptp(df_obs["phaseAngle"]) # * u.deg + ad_params["arc"] = np.ptp(df_obs["midPointMjdTai"]) # * u.d + ad_params["nobs"] = len(df_obs) + ad_params["modelFitMjd"] = Time.now().mjd + # adler_data.populate_phase_parameters(filt, **pc_fit.__dict__) + # TODO: replace any None with np.nan? e.g. phase_parameter_2? + adler_data.populate_phase_parameters(filt, **ad_params) + + # print the test data set + print(adler_data.get_phase_parameters_in_filter(filt, phase_model).__dict__) + + # add to plot + ax1 = fig.axes[0] + # TODO: set colours based on filter + ax1.scatter(df_obs["phaseAngle"], df_obs["reduced_mag"]) + alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg + ax1.plot( + alpha.value, + # pc_fit.ReducedMag(alpha).value, + # label="{}, H={:.2f}, G12={:.2f}".format(filt, pc_fit.H.value, pc_fit.phase_parameter_1), + pc_fit.ReducedMag(alpha), + label="{}, H={:.2f}, {}={:.2f}".format( + filt, pc_fit.H, phase_parameter_1, pc_fit.phase_parameter_1 + ), + ) + ax1.legend() + + # TODO: Use a CLI arg flag to open figure interactively instead of saving? + if plot_show: + plt.show() + # Save figures at the outpath location + else: + fig_file = "{}/phase_curve_{}_{}_{}.png".format( + outpath, ssObjectId, phase_model, int(np.amax(df_obs["midPointMjdTai"])) + ) + msg = "Save figure: {}".format(fig_file) + print(msg) + fig = plot_errorbar(planetoid, fig=fig, filename=fig_file) # TODO: add titles with filter name? + plt.close() + + # analyse colours for the filters provided + + # if requested cycle through the filters, calculating a colour relative to the next filter + if not colour_list: + colour_list = [] + else: + colour_list = colour_list + + # note that the order in which cli_args.filter_list is passed will determine which colours are calculated + for colour in colour_list: + col_filts = colour.split("-") + filt_obs = col_filts[0] + filt_ref = col_filts[1] + + if ~np.isin([filt_obs, filt_ref], planetoid.filter_list).all(): + missing_filts = np.array([filt_obs, filt_ref])[ + ~np.isin([filt_obs, filt_ref], planetoid.filter_list) + ] + continue + + # determine the filt_obs - filt_ref colour + # generate a plot + if plot_show: + plot_dir = None + else: + plot_dir = outpath + + col_dict = col_obs_ref( + planetoid, + adler_data, + phase_model=phase_model, + filt_obs=filt_obs, + filt_ref=filt_ref, + N_ref=N_ref, + # x1 = x1, + # x2 = x2, + plot_dir=plot_dir, + plot_show=plot_show, + ) + + # print(col_dict) + + # Output adler values to a database if a db_name is provided + print(adler_data.__dict__) + if db_name: + adler_db = "{}/{}".format(outpath, db_name) + + if os.path.isfile(adler_db): + print("clear old db file") + os.remove(adler_db) + + msg = "write to {}".format(adler_db) + print(msg) + adler_data.write_row_to_database(adler_db) + +# create the test csv file +fname_csv = "{}/test_SQL_database_table.csv".format(outpath) +conn = sqlite3.Connection(adler_db) +df = pd.read_sql("SELECT * from AdlerData", conn) +print("write to {}".format(fname_csv)) +df.to_csv(fname_csv) diff --git a/tests/data/test_AdlerData_database.db b/tests/data/test_AdlerData_database.db index 31bd8ef..b6138a4 100644 Binary files a/tests/data/test_AdlerData_database.db and b/tests/data/test_AdlerData_database.db differ diff --git a/tests/data/test_SQL_database_table.csv b/tests/data/test_SQL_database_table.csv index 6c34fe1..242eb3d 100644 --- a/tests/data/test_SQL_database_table.csv +++ b/tests/data/test_SQL_database_table.csv @@ -1,2 +1,2 @@ -ssObjectId,timestamp,u_phaseAngle_min,u_phaseAngle_range,u_nobs,u_arc,u_model_1_H,u_model_1_H_err,u_model_1_phase_parameter_1,u_model_1_phase_parameter_1_err,u_model_1_phase_parameter_2,u_model_1_phase_parameter_2_err,u_model_2_H,u_model_2_H_err,u_model_2_phase_parameter_1,u_model_2_phase_parameter_1_err,u_model_2_phase_parameter_2,u_model_2_phase_parameter_2_err,g_phaseAngle_min,g_phaseAngle_range,g_nobs,g_arc,g_model_1_H,g_model_1_H_err,g_model_1_phase_parameter_1,g_model_1_phase_parameter_1_err,g_model_1_phase_parameter_2,g_model_1_phase_parameter_2_err,r_phaseAngle_min,r_phaseAngle_range,r_nobs,r_arc,r_model_2_H,r_model_2_H_err,r_model_2_phase_parameter_1,r_model_2_phase_parameter_1_err,r_model_2_phase_parameter_2,r_model_2_phase_parameter_2_err -8268570668335894776,2024-04-18 13:32:07.096776+00:00,11.0,12.0,13,14.0,15.0,16.0,17.0,18.0,,,25.0,26.0,27.0,28.0,29.0,30.0,31.0,32.0,33,34.0,35.0,36.0,37.0,38.0,,,41.0,42.0,43,44.0,45.0,46.0,47.0,48.0,49.0,50.0 +,ssObjectId,timestamp,g_phaseAngle_min,g_phaseAngle_range,g_nobs,g_arc,g_HG12_Pen16_H,g_HG12_Pen16_H_err,g_HG12_Pen16_phase_parameter_1,g_HG12_Pen16_phase_parameter_1_err,g_HG12_Pen16_phase_parameter_2,g_HG12_Pen16_phase_parameter_2_err,g_HG12_Pen16_modelFitMjd,g_HG_H,g_HG_H_err,g_HG_phase_parameter_1,g_HG_phase_parameter_1_err,g_HG_phase_parameter_2,g_HG_phase_parameter_2_err,g_HG_modelFitMjd,r_phaseAngle_min,r_phaseAngle_range,r_nobs,r_arc,r_HG12_Pen16_H,r_HG12_Pen16_H_err,r_HG12_Pen16_phase_parameter_1,r_HG12_Pen16_phase_parameter_1_err,r_HG12_Pen16_phase_parameter_2,r_HG12_Pen16_phase_parameter_2_err,r_HG12_Pen16_modelFitMjd,r_HG_H,r_HG_H_err,r_HG_phase_parameter_1,r_HG_phase_parameter_1_err,r_HG_phase_parameter_2,r_HG_phase_parameter_2_err,r_HG_modelFitMjd,i_phaseAngle_min,i_phaseAngle_range,i_nobs,i_arc,i_HG12_Pen16_H,i_HG12_Pen16_H_err,i_HG12_Pen16_phase_parameter_1,i_HG12_Pen16_phase_parameter_1_err,i_HG12_Pen16_phase_parameter_2,i_HG12_Pen16_phase_parameter_2_err,i_HG12_Pen16_modelFitMjd,i_HG_H,i_HG_H_err,i_HG_phase_parameter_1,i_HG_phase_parameter_1_err,i_HG_phase_parameter_2,i_HG_phase_parameter_2_err,i_HG_modelFitMjd +0,8268570668335894776,2024-11-06 17:23:21.213435+00:00,27.426668167114258,36.17388343811035,9,3306.0079999999944,20.719465178426468,0.018444134023977106,0.62,,,,60620.724548566934,20.72954332871528,0.018444134023977106,0.15,,,,60620.72454979901,2.553332567214966,124.23803400993347,38,3338.0655999999944,19.92863542616601,0.018525355171274356,1.0,0.05300829494059732,,,60620.72454870314,18.879873513575124,0.007456882346021157,-0.253,1.6701987135262256e-05,,,60620.72454985305,10.0595064163208,101.74150371551514,32,3342.0585599999977,19.653155995892117,0.01415178691419005,1.0,0.01516569963597136,,,60620.72454882437,18.628894992991583,0.01716803531553185,-0.253,0.00014324314956736168,,,60620.724549914594