diff --git a/notebooks/colour_functions_demo.ipynb b/notebooks/colour_functions_demo.ipynb new file mode 100644 index 0000000..902f481 --- /dev/null +++ b/notebooks/colour_functions_demo.ipynb @@ -0,0 +1,391 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "d1ef57c2", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61f50e0d", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82537a61", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.dataclasses.AdlerData import AdlerData\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "from adler.science.Colour import col_obs_ref\n", + "from adler.utilities.plotting_utilities import plot_errorbar\n", + "from adler.utilities.science_utilities import apparition_gap_finder, get_df_obs_filt\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import astropy.units as u\n", + "from scipy.interpolate import splrep, BSpline, CubicSpline\n", + "from astropy.modeling.fitting import SLSQPLSQFitter\n", + "from astropy.time import Time\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9f6dcda", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"6098332225018\" # good MBA test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6f705c1", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the object data\n", + "# fname = \"../tests/data/testing_database.db\"\n", + "fname = \"/Users/jrobinson/lsst-adler/notebooks/gen_test_data/adler_demo_testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21e0b93e", + "metadata": {}, + "outputs": [], + "source": [ + "# check orbit parameters\n", + "e = planetoid.MPCORB.e\n", + "incl = planetoid.MPCORB.incl\n", + "q = planetoid.MPCORB.q\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee4064e7", + "metadata": {}, + "outputs": [], + "source": [ + "# plot and fit the phase curves in g and r filters\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "adler_cols = AdlerData(ssoid, planetoid.filter_list)\n", + "\n", + "for filt in [\"g\", \"r\"]:\n", + " # get observations and define a phase curve model to fit\n", + "\n", + " sso = planetoid.SSObject_in_filter(filt)\n", + " obs = planetoid.observations_in_filter(filt)\n", + "\n", + " H = sso.H\n", + " G12 = sso.G12\n", + " pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name=\"HG12_Pen16\")\n", + "\n", + " alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg\n", + "\n", + " pc_fit = pc.FitModel(\n", + " np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n", + " np.array(getattr(obs, \"reduced_mag\")) * u.mag,\n", + " )\n", + " pc = pc.InitModelSbpy(pc_fit)\n", + " red_mag = pc.ReducedMag(alpha)\n", + "\n", + " adler_cols.populate_phase_parameters(filt, **pc.ReturnModelDict())\n", + "\n", + " # add this phase curve to the figure\n", + " fig = plot_errorbar(planetoid, filt_list=[filt], fig=fig)\n", + " ax1 = fig.axes[0]\n", + " ax1.plot(\n", + " alpha.value,\n", + " pc.ReducedMag(alpha).value,\n", + " label=\"{}: H={:.2f}, G12={:.2f}\".format(filt, pc.H, pc.phase_parameter_1),\n", + " )\n", + " ax1.legend()\n", + "\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0df3c51", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the r filter phase curve model\n", + "adler_cols.get_phase_parameters_in_filter(\"r\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c46fb9eb", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the g filter phase curve model\n", + "adler_cols.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73871e56", + "metadata": {}, + "outputs": [], + "source": [ + "# get the boundary times for each apparation of the object in the survey\n", + "df_obs_all = pd.DataFrame()\n", + "for filt in [\"r\", \"g\"]:\n", + " obs = planetoid.observations_in_filter(filt)\n", + " _df_obs = pd.DataFrame(obs.__dict__)\n", + " df_obs_all = pd.concat([df_obs_all, _df_obs])\n", + "df_obs_all = df_obs_all.sort_values(\"midPointMjdTai\")\n", + "\n", + "t_app = apparition_gap_finder(np.array(df_obs_all[\"midPointMjdTai\"]))\n", + "print(t_app)" + ] + }, + { + "cell_type": "markdown", + "id": "63980010", + "metadata": {}, + "source": [ + "Now we can inpsect how the colour of the object varies (or not) as a function of time. The adler function `col_obs_ref` will compare the latest observation in a given filter with observations in another filter. By setting parameter `N_ref` one can set how many past obsevrations to use when calculating the latest colour." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd0e0185", + "metadata": {}, + "outputs": [], + "source": [ + "# Here we simulate observations coming night-by-night and the calculation of a g-r colour for the object\n", + "\n", + "# define colour function parameters\n", + "\n", + "# set number of reference observations to use for colour estimate\n", + "# N_ref = 5\n", + "N_ref = 3\n", + "# N_ref = 1\n", + "# N_ref = None\n", + "\n", + "# observation and filter field names\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = \"AbsMag\"\n", + "y_plot = \"reduced_mag\"\n", + "yerr_plot = \"magErr\"\n", + "filt_obs = \"g\" # observation filter\n", + "filt_ref = \"r\" # reference filter (we are calculating a filt_obs - filt_ref colour)\n", + "\n", + "# define colour field names\n", + "colour = \"{}-{}\".format(filt_obs, filt_ref)\n", + "colErr = \"{}-{}Err\".format(filt_obs, filt_ref)\n", + "delta_t_col = \"delta_t_{}\".format(colour)\n", + "y_ref_col = \"{}_{}\".format(y_plot, filt_ref)\n", + "x1_ref_col = \"{}1_{}\".format(x_plot, filt_ref)\n", + "x2_ref_col = \"{}2_{}\".format(x_plot, filt_ref)\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "col_dict_list = []\n", + "for app_i in range(len(t_app) - 1):\n", + " # consider only one apparition\n", + " if app_i != 3:\n", + " continue\n", + "\n", + " time_min = t_app[app_i]\n", + " time_max = t_app[app_i + 1]\n", + " print(time_min, time_max)\n", + "\n", + " _df_obs_all = df_obs_all[\n", + " (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\n", + " print(app_i, len(_df_obs_all))\n", + " _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n", + "\n", + " # get the phase curve model and observations for each filter\n", + "\n", + " # get the stored AdlerData parameters for the observation filter\n", + " ad_g = adler_cols.get_phase_parameters_in_filter(filt_obs, \"HG12_Pen16\")\n", + " pc_g = PhaseCurve().InitModelDict(ad_g.__dict__) # make the PhaseCurve object from AdlerData\n", + " # get the phase curve model for the reference filter\n", + " ad_r = adler_cols.get_phase_parameters_in_filter(filt_ref, \"HG12_Pen16\")\n", + " pc_r = PhaseCurve().InitModelDict(ad_r.__dict__)\n", + " # get the observations in both filters\n", + " df_obs = get_df_obs_filt(\n", + " planetoid, filt_obs, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_g\n", + " )\n", + " df_obs_ref = get_df_obs_filt(\n", + " planetoid, filt_ref, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_r\n", + " )\n", + "\n", + " ax1.errorbar(df_obs[x_plot], df_obs[y_plot], df_obs[yerr_plot], fmt=\"o\", label=filt_obs)\n", + " ax1.errorbar(df_obs_ref[x_plot], df_obs_ref[y_plot], df_obs_ref[yerr_plot], fmt=\"o\", label=filt_ref)\n", + "\n", + " # simulate stepping through each filt_obs observation\n", + " x1 = time_min\n", + " for xi in range(len(df_obs)):\n", + " x2 = df_obs.iloc[xi][x_plot]\n", + " print(xi, x1, x2)\n", + "\n", + " # do the colour finding function here\n", + " col_dict = col_obs_ref(\n", + " planetoid,\n", + " adler_cols,\n", + " filt_obs=filt_obs,\n", + " filt_ref=filt_ref,\n", + " N_ref=N_ref,\n", + " x_col=x_plot,\n", + " y_col=y_plot,\n", + " yerr_col=yerr_plot,\n", + " x1=x1,\n", + " x2=x2,\n", + " )\n", + " col_dict_list.append(col_dict)\n", + "\n", + " # plot some lines to show the colour and mean reference\n", + " ax1.vlines(df_obs.iloc[xi][x_plot], df_obs.iloc[xi][y_plot], col_dict[y_ref_col], color=\"k\", ls=\":\")\n", + " ax1.hlines(col_dict[y_ref_col], col_dict[x1_ref_col], col_dict[x2_ref_col], color=\"k\", ls=\"--\")\n", + "\n", + "# store running colour parameters as a dataframe\n", + "df_col = pd.DataFrame(col_dict_list)\n", + "df_col = pd.concat([df_obs, df_col], axis=1)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7734122", + "metadata": {}, + "outputs": [], + "source": [ + "# recorded colour parameters\n", + "df_col" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cd867d9", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot how the colour changes as a function of time\n", + "# find filt_obs - filt_ref of newest filt_obs observation to the mean of the previous N_ref filt_ref observations\n", + "# colour code by time diff between obs and most recent obs_ref\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot = colour\n", + "y_plot_err = colErr\n", + "c_plot = delta_t_col\n", + "df_plot = df_col\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "\n", + "s1 = ax1.scatter(df_plot[x_plot], df_plot[y_plot], c=df_plot[c_plot], zorder=3)\n", + "cbar1 = plt.colorbar(s1)\n", + "ax1.errorbar(df_plot[x_plot], df_plot[y_plot], df_plot[yerr_plot], fmt=\".\", zorder=1)\n", + "\n", + "# TODO: leave the most recent obs out of this mean?\n", + "# loop through new obs and calculate current mean and std?\n", + "obs_ref_mean = np.mean(df_plot[y_plot])\n", + "obs_ref_median = np.median(df_plot[y_plot])\n", + "obs_ref_std = np.std(df_plot[y_plot])\n", + "print(obs_ref_median, obs_ref_std)\n", + "\n", + "ax1.axhline(obs_ref_mean, c=\"k\")\n", + "ax1.axhspan(obs_ref_mean - obs_ref_std, obs_ref_mean + obs_ref_std, zorder=0, color=\"k\", alpha=0.2)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "cbar1.set_label(c_plot)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05abfee3", + "metadata": {}, + "outputs": [], + "source": [ + "# These colours can then be run through the previously written outlier detection functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "605a043a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "adler-dev", + "language": "python", + "name": "adler-dev" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/gen_test_data/adler_demo_testing_database.db b/notebooks/gen_test_data/adler_demo_testing_database.db index 1c3cdb5..455b28b 100644 Binary files a/notebooks/gen_test_data/adler_demo_testing_database.db and b/notebooks/gen_test_data/adler_demo_testing_database.db differ diff --git a/src/adler/dataclasses/AdlerData.py b/src/adler/dataclasses/AdlerData.py index 86c2e2a..e3231f6 100644 --- a/src/adler/dataclasses/AdlerData.py +++ b/src/adler/dataclasses/AdlerData.py @@ -101,7 +101,7 @@ def populate_phase_parameters(self, filter_name, **kwargs): # update the value if it's in **kwargs for model_key in MODEL_DEPENDENT_KEYS: - if kwargs.get(model_key): + if model_key in kwargs: setattr( self.filter_dependent_values[filter_index].model_dependent_values[model_index], model_key, diff --git a/src/adler/dataclasses/AdlerPlanetoid.py b/src/adler/dataclasses/AdlerPlanetoid.py index b9d31b8..4592e18 100644 --- a/src/adler/dataclasses/AdlerPlanetoid.py +++ b/src/adler/dataclasses/AdlerPlanetoid.py @@ -240,7 +240,9 @@ def populate_observations( observations_sql_query = f""" SELECT ssObject.ssObjectId, mag, magErr, band, midPointMjdTai, ra, dec, phaseAngle, - topocentricDist, heliocentricDist + topocentricDist, heliocentricDist, heliocentricX, heliocentricY, heliocentricZ, + topocentricX, topocentricY, topocentricZ, + eclipticLambda, eclipticBeta FROM {schema}ssObject JOIN {schema}diaSource ON {schema}ssObject.ssObjectId = {schema}diaSource.ssObjectId @@ -291,7 +293,7 @@ def populate_MPCORB(self, ssObjectId, service=None, sql_filename=None, schema="d MPCORB_sql_query = f""" SELECT - ssObjectId, mpcDesignation, mpcNumber, mpcH, mpcG, epoch, peri, node, incl, e, n, q, + ssObjectId, mpcDesignation, fullDesignation, mpcNumber, mpcH, mpcG, epoch, tperi, peri, node, incl, e, n, q, uncertaintyParameter, flags FROM {schema}MPCORB diff --git a/src/adler/dataclasses/MPCORB.py b/src/adler/dataclasses/MPCORB.py index 289084d..844cef9 100644 --- a/src/adler/dataclasses/MPCORB.py +++ b/src/adler/dataclasses/MPCORB.py @@ -4,10 +4,12 @@ MPCORB_KEYS = { "mpcDesignation": str, + "fullDesignation": str, "mpcNumber": int, "mpcH": float, "mpcG": float, "epoch": float, + "tperi": float, "peri": float, "node": float, "incl": float, @@ -31,6 +33,9 @@ class MPCORB: mpcDesignation: str Number or provisional designation (in packed form) + fullDesignation: str + Number or provisional designation (in readable form) + mpcNumber: int MPC number (if the asteroid has been numbered; NULL otherwise). Provided for convenience. @@ -43,6 +48,9 @@ class MPCORB: epoch: float Epoch (in MJD, .0 TT) + tperi: float + MJD of pericentric passage + peri: float Argument of perihelion, J2000.0 (degrees) @@ -71,10 +79,12 @@ class MPCORB: ssObjectId: str = "" mpcDesignation: str = "" + fullDesignation: str = "" mpcNumber: int = 0 mpcH: float = 0.0 mpcG: float = 0.0 epoch: float = 0.0 + tperi: float = 0.0 peri: float = 0.0 node: float = 0.0 incl: float = 0.0 diff --git a/src/adler/dataclasses/Observations.py b/src/adler/dataclasses/Observations.py index e99a758..eb92af2 100644 --- a/src/adler/dataclasses/Observations.py +++ b/src/adler/dataclasses/Observations.py @@ -12,6 +12,14 @@ "phaseAngle": np.ndarray, "topocentricDist": np.ndarray, "heliocentricDist": np.ndarray, + "heliocentricX": np.ndarray, + "heliocentricY": np.ndarray, + "heliocentricZ": np.ndarray, + "topocentricX": np.ndarray, + "topocentricY": np.ndarray, + "topocentricZ": np.ndarray, + "eclipticLambda": np.ndarray, + "eclipticBeta": np.ndarray, } @@ -53,6 +61,30 @@ class Observations: heliocentricDist: array_like of floats Heliocentric distance. + heliocentricX: array_like of floats + x-axis component of the heliocentric distance. + + heliocentricY: array_like of floats + y-axis component of the heliocentric distance. + + heliocentricZ: array_like of floats + z-axis component of the heliocentric distance. + + topocentricX: array_like of floats + x-axis component of the topocentric distance. + + topocentricY: array_like of floats + y-axis component of the topocentric distance. + + topocentricZ: array_like of floats + z-axis component of the topocentric distance. + + eclipticLambda: array_like of floats + The ecliptic longitude. + + eclipticBeta: array_like of floats + The ecliptic latitude. + reduced_mag: array_like of floats The reduced magnitude. @@ -71,6 +103,14 @@ class Observations: phaseAngle: np.ndarray = field(default_factory=lambda: np.zeros(0)) topocentricDist: np.ndarray = field(default_factory=lambda: np.zeros(0)) heliocentricDist: np.ndarray = field(default_factory=lambda: np.zeros(0)) + heliocentricX: np.ndarray = field(default_factory=lambda: np.zeros(0)) + heliocentricY: np.ndarray = field(default_factory=lambda: np.zeros(0)) + heliocentricZ: np.ndarray = field(default_factory=lambda: np.zeros(0)) + topocentricX: np.ndarray = field(default_factory=lambda: np.zeros(0)) + topocentricY: np.ndarray = field(default_factory=lambda: np.zeros(0)) + topocentricZ: np.ndarray = field(default_factory=lambda: np.zeros(0)) + eclipticLambda: np.ndarray = field(default_factory=lambda: np.zeros(0)) + eclipticBeta: np.ndarray = field(default_factory=lambda: np.zeros(0)) reduced_mag: np.ndarray = field(default_factory=lambda: np.zeros(0)) num_obs: int = 0 diff --git a/src/adler/science/Colour.py b/src/adler/science/Colour.py new file mode 100644 index 0000000..0144684 --- /dev/null +++ b/src/adler/science/Colour.py @@ -0,0 +1,126 @@ +import numpy as np +from adler.science.PhaseCurve import PhaseCurve +from adler.utilities.science_utilities import get_df_obs_filt + +# TODO: add colour estimate function + + +def col_obs_ref( + planetoid, + adler_cols, + filt_obs="g", + filt_ref="r", + N_ref=3, + x_col="midPointMjdTai", + y_col="AbsMag", + yerr_col="magErr", + x1=None, + x2=None, +): + """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). + By setting N_ref one can control how many reference observations to include in the colour calculation. + Note that the observations are considered to be in reduced magnitudes, hence why an adlerData object with phase curve models is passed. + + Parameters + ----------- + planetoid : object + Adler planetoid object + + adler_cols : object + AdlerData object + + filt_obs: str + Filter name of the new observation for which to calculate a filt_obs - filt_ref colour + + filt_ref:str + Filter name of the reference observations when calculating the filt_obs - filt_ref colour + + N_ref: int + Number of reference observations to use when calculated the reference absolute magnitude. + Set to 1 to use only the most recent of ref observations. + Set to None to use all past ref obs. + + x_col: str + Time (or phase angle) column name + + y_col:str + Magnitude column name (e.g. AbsMag or reduced_mag) + + yerr_col: str + Magnitude uncertainty column name + + x1: float + Lower limit on x_col values (>=) + + x2: float + Upper limit on x_col values (<=) + + Returns + ---------- + + col_dict : dict + Dictionary containing the colour information for the most recent obs in filt_obs + + """ + # define fields to be recorded as a colour dictionary + colour = "{}-{}".format(filt_obs, filt_ref) + colErr = "{}-{}Err".format(filt_obs, filt_ref) + delta_t_col = "delta_t_{}".format(colour) + y_ref_col = "{}_{}".format(y_col, filt_ref) + x1_ref_col = "{}1_{}".format(x_col, filt_ref) + x2_ref_col = "{}2_{}".format(x_col, filt_ref) + + # get the stored AdlerData parameters for this filter # TODO: do this bit better + ad_obs = adler_cols.get_phase_parameters_in_filter(filt_obs, "HG12_Pen16") + # make the PhaseCurve object from AdlerData + pc_obs = PhaseCurve().InitModelDict(ad_obs.__dict__) + ad_ref = adler_cols.get_phase_parameters_in_filter(filt_ref, "HG12_Pen16") + pc_ref = PhaseCurve().InitModelDict(ad_ref.__dict__) + df_obs = get_df_obs_filt(planetoid, filt_obs, x1=x1, x2=x2, col_list=[y_col, yerr_col], pc_model=pc_obs) + df_obs_ref = get_df_obs_filt( + planetoid, filt_ref, x1=x1, x2=x2, col_list=[y_col, yerr_col], pc_model=pc_ref + ) + + # select the values of the new observation in the selected filter + i = -1 + x_obs = df_obs.iloc[i][x_col] + y_obs = df_obs.iloc[i][y_col] + yerr_obs = df_obs.iloc[i][yerr_col] + + # select observations in the reference filter from before the new obs + ref_mask = df_obs_ref[x_col] < x_obs + + # set the number of ref obs to use + if N_ref is None: + _N_ref = len(df_obs_ref[ref_mask]) # use all available ref obs + else: + _N_ref = N_ref + + # select only the N_ref ref obs for comparison + _df_obs_ref = df_obs_ref[ref_mask].iloc[-_N_ref:] + if len(_df_obs_ref) == 0: + print("no reference observations") # TODO: add proper error handling and logging here + return df_obs + + # determine reference observation values + y_ref = np.mean(_df_obs_ref[y_col]) # TODO: add option to choose statistic, e.g. mean or median? + yerr_ref = np.std(_df_obs_ref[y_col]) # TODO: propagate ref uncertainty properly + # determine the ref obs time range + x1_ref = np.array(_df_obs_ref[x_col])[0] + x2_ref = np.array(_df_obs_ref[x_col])[-1] + + # Create the colour dict + col_dict = {} + col_dict[colour] = y_obs - y_ref + col_dict[delta_t_col] = x_obs - x2_ref + col_dict[colErr] = np.sqrt((yerr_obs**2.0) + (yerr_ref**2.0)) + col_dict[y_ref_col] = y_ref + col_dict[x1_ref_col] = x1_ref + col_dict[x2_ref_col] = x2_ref + + # TODO: + # could also record phase angle diff and phase curve residual? + # need to test error case where there are no r filter obs yet + + return col_dict diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index a40d1cb..f0b0463 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -158,9 +158,18 @@ def InitModelDict(self, model_dict): """ - model = PhaseCurve() + # clean the input dictionary + del_keys = [] for key, value in model_dict.items(): - setattr(model, key, value) + if not hasattr(self, key): + del_keys.append(key) + model_dict = model_dict.copy() # make a copy to avoid changing the original dict + for key in del_keys: + model_dict.pop(key, None) + + # initialise a new model + model = PhaseCurve(**model_dict) + return model def InitModelSbpy(self, model_sbpy): @@ -221,7 +230,8 @@ def ReducedMag(self, phase_angle): Parameters ----------- phase_angle : float or array - value or array of phase angles at which to evaluate the phasecurve model, must have astropy units of degrees. + value or array of phase angles at which to evaluate the phasecurve model. + Must have astropy units of degrees if sbpy model uses units, otherwise pass radian values. Returns ---------- @@ -233,6 +243,59 @@ def ReducedMag(self, phase_angle): return self.model_function(phase_angle) + def ModelResiduals(self, phase_angle, reduced_mag): + """For a set of phase curve observations, return the residuals to the PhaseCurve model. + NB that units must match the sbpy model. E.g. phase_angle should be passed with units of degrees, or be in values of radians + + Parameters + ----------- + phase_angle : float or array + value or array of phase angles at which to evaluate the phasecurve model. + + reduced_mag : float or array + value or array of reduced magnitudes at which to evaluate the phasecurve model. + + Returns + ---------- + + residuals : float or array + The residuals of the observations minus PhaseCurve model values + + """ + + residuals = reduced_mag - self.ReducedMag(phase_angle) + + return residuals + + def AbsMag(self, phase_angle, reduced_mag): + """For a set of phase curve observations, return the absolute magnitude from the fitted phase curve model. + I.e. this is the model residuals, shifted by the fitted absolute magnitude + NB that units for phase_angle and reduced_mag must match the sbpy model. E.g. phase_angle should be passed with units of degrees, or be in values of radians + + Parameters + ----------- + phase_angle : float or array + value or array of phase angles at which to evaluate the phasecurve model. + + reduced_mag : float or array + value or array of reduced magnitudes at which to evaluate the phasecurve model. + + Returns + ---------- + + abs_mag : float or array + The residuals of the observations minus PhaseCurve model values shifted by the model absolute magnitude + + """ + + # TODO: add option to pass model & filt instead of reduced_mag & phase angle - will calculate all absmags automatically + # probably not possible as observations are contained in a separate object + # if (phase_angle is None) and (reduced_mag is None) and (filt is not None)... + + abs_mag = reduced_mag - self.ReducedMag(phase_angle) + self.H + + return abs_mag + def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): """Fit the phasecurve model parameters to observations. starts with a phase curve model as an initial guess for parameters. diff --git a/src/adler/utilities/science_utilities.py b/src/adler/utilities/science_utilities.py index 68beaf7..90ff1df 100644 --- a/src/adler/utilities/science_utilities.py +++ b/src/adler/utilities/science_utilities.py @@ -1,10 +1,15 @@ import numpy as np +import pandas as pd from astropy.stats import sigma_clip as astropy_sigma_clip +import astropy.units as u def outlier_diff(new_res, diff_cut=1.0): """Test whether new data point(s) is an outlier compared to the model by considering the residual difference. + Parameters + ----------- + new_res: array The residuals of the new data points compared to the model diff_cut: float @@ -29,6 +34,9 @@ def outlier_diff(new_res, diff_cut=1.0): def outlier_std(new_res, data_res, std_cut=3.0): """Test whether new data point(s) is an outlier compared to the model by considering the standard deviation of the residuals. + Parameters + ----------- + new_res: array The residuals of the new data point(s) compared to the model data_res: array @@ -59,6 +67,9 @@ def zero_func(x, axis=None): """Dummy function to return a zero. Can be used as the centre function in astropy.stats.sigma_clip to get std relative to zero rather than median/mean value. + Parameters + ----------- + x: Dummy variable axis: @@ -70,6 +81,9 @@ def zero_func(x, axis=None): def sigma_clip(data_res, kwargs={"maxiters": 1, "cenfunc": zero_func}): """Wrapper function for astropy.stats.sigma_clip, here we define the default centre of the data (the data - model residuals) to be zero + Parameters + ----------- + data_res: array The residuals of the data compared to the model. kwargs: dict @@ -90,6 +104,9 @@ def sigma_clip(data_res, kwargs={"maxiters": 1, "cenfunc": zero_func}): def outlier_sigma_diff(data_res, data_sigma, std_sigma=1): """Function to identify outliers by comparing the uncertainty of measurements to their residuals + Parameters + ----------- + data_res: array The residuals of the data compared to the model. data_sigma: array @@ -106,3 +123,101 @@ def outlier_sigma_diff(data_res, data_sigma, std_sigma=1): outlier_flag = np.abs(data_res) > (std_sigma * data_sigma) return outlier_flag + + +def apparition_gap_finder(x, dx=100.0): + """Function to find gaps in a data series. E.g. given an array of observation times, find the different apparitions of an asteroid from gaps between observations larger than a given value. + + Parameters + ----------- + + x: array + The SORTED data array to search for gaps + dx: float + The size of gap to identify in data series + + Returns + ----------- + x_gaps: array + Values of x which define the groups in the data, where each group is x_gaps[i] <= x < x_gaps[i+1] + """ + + ## TODO: check that x is sorted? + + # find the difference between each data point + x_diff = np.diff(x) + + # select only differences greater than value dx + x_diff_mask = x_diff > dx + + # find the values in x that correspond to the start of each group + x_gaps = x[1:][x_diff_mask] + + # add the very first observation to the array + x_gaps = np.insert(x_gaps, 0, x[0]) + + # add the very last observation to the array, including a shift of dx to ensure the boundary inequalities work + x_gaps = np.insert(x_gaps, len(x_gaps), x[-1] + dx) + + return x_gaps + + +def get_df_obs_filt(planetoid, filt, x_col="midPointMjdTai", x1=None, x2=None, col_list=None, pc_model=None): + """Retrieve a dataframe of observations in a given filter. Has the option to limit the observations to a range of values, e.g. times/phase angles, if required. + + Parameters + ----------- + + planetoid: object + Adler planetoid object containging the observations + filt: str + The filter to query + x_col: str + Column name to use for ordering values and limiting observations to a range: x1 <= df_obs[x_col] <= x2 + x1: float + Lower limit value for x_col + x2: float + Upper limit value for x_col + col_list: list + List of column names to retrieve in addition to x_col, otherwise all columns are retrieved. + N.B. if AbsMag is included in col_list then the absolute magnitude is calculated for the phase curve model pc_model + pc_model: object + Adler PhaseCurve model used to calculate AbsMag if required + + Returns + ----------- + df_obs: DataFrame + DataFrame of observations in the requested filter, ordered by x_col and with any x_col limits applied + """ + + # get observations in filter as a dataframe + obs = planetoid.observations_in_filter(filt) + # TODO: add option to get all filters? calculating AbsMag gets complicated... + # TODO: split the dataframe functions into a separate function? + df_obs = pd.DataFrame(obs.__dict__) + + # if no list of columns is provided, retrieve all columns + if col_list is None: + col_list = list(df_obs) + else: + col_list = [x_col] + col_list + + # calculate the phase curve absolute magnitudes (i.e. remove phase curve variation from reduced_mag) + 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 + + # select only the required columns + df_obs = df_obs[col_list] + + # sort into x_col order + df_obs = df_obs.sort_values(x_col) + # limit to range in x_col if limits x1, x2 are supplied + if (x1 is not None) & (x2 is not None): # TODO: add functionality to set only upper/lower limit + x_mask = (df_obs[x_col] >= x1) & (df_obs[x_col] <= x2) + df_obs = df_obs[x_mask] + # reset the dataframe index + df_obs = df_obs.reset_index(drop=True) + + return df_obs diff --git a/tests/adler/dataclasses/test_MPCORB.py b/tests/adler/dataclasses/test_MPCORB.py index d139d00..ee43b3d 100644 --- a/tests/adler/dataclasses/test_MPCORB.py +++ b/tests/adler/dataclasses/test_MPCORB.py @@ -12,7 +12,7 @@ def test_construct_MPCORB_from_data_table(): test_query = f""" SELECT - ssObjectId, mpcDesignation, mpcNumber, mpcH, mpcG, epoch, peri, node, incl, e, n, q, + ssObjectId, mpcDesignation, fullDesignation, mpcNumber, mpcH, mpcG, epoch, tperi, peri, node, incl, e, n, q, uncertaintyParameter, flags FROM MPCORB @@ -25,6 +25,9 @@ def test_construct_MPCORB_from_data_table(): assert test_MPCORB.ssObjectId == 8268570668335894776 assert test_MPCORB.mpcDesignation == "2014 QL4" + assert ( + test_MPCORB.fullDesignation == "2011 2014 QL433" + ) # N.B. that there are known DP0.3 issues with mpcDesignation and fullDesignation, https://dp0-3.lsst.io/data-products-dp0-3/data-simulation-dp0-3.html#known-issues assert test_MPCORB.mpcNumber == 0 assert_almost_equal(test_MPCORB.mpcH, 19.8799991607666, decimal=6) assert_almost_equal(test_MPCORB.mpcG, 0.15000000596046448, decimal=6) diff --git a/tests/adler/dataclasses/test_Observations.py b/tests/adler/dataclasses/test_Observations.py index 44a3717..39e1a57 100644 --- a/tests/adler/dataclasses/test_Observations.py +++ b/tests/adler/dataclasses/test_Observations.py @@ -15,7 +15,9 @@ def test_construct_observations_from_data_table(): test_query = f""" SELECT ssObject.ssObjectId, mag, magErr, band, midPointMjdTai, ra, dec, phaseAngle, - topocentricDist, heliocentricDist + topocentricDist, heliocentricDist, heliocentricX, heliocentricY, heliocentricZ, + topocentricX, topocentricY, topocentricZ, + eclipticLambda, eclipticBeta FROM ssObject JOIN diaSource ON ssObject.ssObjectId = diaSource.ssObjectId diff --git a/tests/adler/science/test_Colour.py b/tests/adler/science/test_Colour.py new file mode 100644 index 0000000..08e64c6 --- /dev/null +++ b/tests/adler/science/test_Colour.py @@ -0,0 +1,112 @@ +import numpy as np +import pandas as pd +from numpy.testing import assert_almost_equal, assert_array_equal, assert_array_almost_equal +import pytest +import os +import astropy.units as u + +from adler.science.Colour import col_obs_ref +from adler.utilities.tests_utilities import get_test_data_filepath +from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid +from adler.dataclasses.AdlerData import AdlerData +from adler.science.PhaseCurve import PhaseCurve +from adler.utilities.science_utilities import get_df_obs_filt + +# set up a good test object (main belt asteroid) +ssoid = 6098332225018 +test_db_path = get_test_data_filepath("testing_database.db") +planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=test_db_path) +adler_cols = AdlerData(ssoid, planetoid.filter_list) + +# Set time boundaries for a given apparition (apparition index 3 for this object) +test_app_i = 3 +t_app1 = 61469 +t_app2 = 61678 + +# define columns +column_dict = { + "x_col": "midPointMjdTai", + "y_col": "AbsMag", + "y_col": "reduced_mag", + "yerr_col": "magErr", + "filt_obs": "g", # observation filter + "filt_ref": "r", # reference filter (we are calculating a filt_obs - filt_ref colour) +} + +# Load the stored colour results for this apparition from a df csv in tests/data, see also lsst-adler/notebooks/colour_functions_testing.ipynb +# test_data_path = "/Users/jrobinson/lsst-adler/tests/data" # pytest does not seem to like relative file paths, use absolute +test_data_path = "/".join(os.path.abspath(__file__).split("/")[:-3]) + "/data" +df_col_fname = "{}/df_{}_{}_{}_app_{}".format( + test_data_path, ssoid, column_dict["filt_obs"], column_dict["filt_ref"], test_app_i +) +df_N_ref_1 = pd.read_csv("{}_N_ref_{}.csv".format(df_col_fname, 1), index_col=0) +df_N_ref_3 = pd.read_csv("{}_N_ref_{}.csv".format(df_col_fname, 3), index_col=0) +df_N_ref_5 = pd.read_csv("{}_N_ref_{}.csv".format(df_col_fname, 5), index_col=0) + +# fit adler phase curve models +# TODO: replace this with the stored alder_cols when they are implemented in the database +adler_cols = AdlerData(ssoid, planetoid.filter_list) +for filt in [column_dict["filt_obs"], column_dict["filt_ref"]]: + sso = planetoid.SSObject_in_filter(filt) + obs = planetoid.observations_in_filter(filt) + + H = sso.H + G12 = sso.G12 + pc = PhaseCurve(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, + ) + pc = pc.InitModelSbpy(pc_fit) + + adler_cols.populate_phase_parameters(filt, **pc.ReturnModelDict()) + + +def test_col_obs_ref( + planetoid=planetoid, + adler_cols=adler_cols, + 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], +): + # gather the observations for the apparition + ad_g = adler_cols.get_phase_parameters_in_filter(column_dict["filt_obs"], "HG12_Pen16") + pc_g = PhaseCurve().InitModelDict(ad_g.__dict__) + df_obs = get_df_obs_filt( + planetoid, + column_dict["filt_obs"], + x1=t_app1, + x2=t_app2, + col_list=[column_dict["y_col"], column_dict["yerr_col"]], + pc_model=pc_g, + ) + + # test multiple N_ref + for N_ref, df_ref in zip(N_ref_list, df_ref_list): + col_dict_list = [] + + # simulate stepping through each filt_obs observation + x1 = t_app1 + for xi in range(len(df_obs)): + x2 = df_obs.iloc[xi][column_dict["x_col"]] + + # do the colour finding function here + col_dict = col_obs_ref( + planetoid, + adler_cols, + N_ref=N_ref, + x1=x1, + x2=x2, + **column_dict, + ) + col_dict_list.append(col_dict) + + # store results as a dataframe + df_col = pd.DataFrame(col_dict_list) + df_col = pd.concat([df_obs, df_col], axis=1) + + # compare results df to stored df + for x in list(df_col): + print(x) + assert_array_almost_equal(np.array(df_col[x]), np.array(df_ref[x])) diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index 77e0385..d67fc3c 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/tests/adler/science/test_PhaseCurve.py @@ -143,3 +143,6 @@ def test_PhaseCurve_FitModel_HG_bounds(): 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 + + +# TODO: test absmag diff --git a/tests/adler/utilities/test_science_utilities.py b/tests/adler/utilities/test_science_utilities.py index 3a2d753..b8a0588 100644 --- a/tests/adler/utilities/test_science_utilities.py +++ b/tests/adler/utilities/test_science_utilities.py @@ -220,3 +220,11 @@ def test_sigma_clip(): def test_outlier_sigma_diff(): outliers = sci_utils.outlier_sigma_diff(y_res, y_err, std_sigma=5.0) assert_array_equal(outliers, outliers_y) + + +# TODO: test apparition_gap_finder, correct number of groups identified? correct boundaries? +# check boundaires are in ascending order + +# TODO: test get_df_obs_filt, check filter and sorting/time range constraints, +# check the column list behaviour +# check AbsMag calculation diff --git a/tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv b/tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv new file mode 100644 index 0000000..9d0ab62 --- /dev/null +++ b/tests/data/df_6098332225018_g_r_app_3_N_ref_1.csv @@ -0,0 +1,16 @@ +,midPointMjdTai,reduced_mag,magErr,g-r,delta_t_g-r,g-rErr,reduced_mag_r,midPointMjdTai1_r,midPointMjdTai2_r +0,61485.3588,17.82546893645328,0.13199999928474426,0.6294446491054693,7.961480000005395,0.13199999928474426,17.19602428734781,61477.39732,61477.39732 +1,61500.36601,17.548929910716545,0.052000001072883606,0.4263135941635241,3.041850000001432,0.052000001072883606,17.12261631655302,61497.32416,61497.32416 +2,61504.33885,17.545916408017177,0.05700000002980232,0.42122628677132923,0.024510000002919696,0.05700000002980232,17.124690121245848,61504.31434,61504.31434 +3,61524.2819,17.60109262744087,0.06199999898672104,0.4764025061950221,19.967560000004596,0.06199999898672104,17.124690121245848,61504.31434,61504.31434 +4,61525.31814,17.530157869823057,0.04699999839067459,0.4042173365957282,0.024610000000393484,0.04699999839067459,17.12594053322733,61525.29353,61525.29353 +5,61536.27685,17.474481419144603,0.039000000804662704,0.3952621411742143,4.973070000000007,0.039000000804662704,17.079219277970388,61531.30378,61531.30378 +6,61556.37692,17.262436542715754,0.03200000151991844,0.39311468820022455,0.024060000003373716,0.03200000151991844,16.86932185451553,61556.35286,61556.35286 +7,61558.20527,17.330018817677058,0.03200000151991844,0.46069696316152786,1.8524099999995087,0.03200000151991844,16.86932185451553,61556.35286,61556.35286 +8,61558.20674,17.302025752576462,0.029999999329447746,0.43270389806093235,1.8538800000023912,0.029999999329447746,16.86932185451553,61556.35286,61556.35286 +9,61589.25295,17.086183244116974,0.028999999165534973,0.5580614970225177,1.968010000004142,0.028999999165534973,16.528121747094456,61587.28494,61587.28494 +10,61589.25519,17.035176049202192,0.032999999821186066,0.5070543021077363,1.9702500000057626,0.032999999821186066,16.528121747094456,61587.28494,61587.28494 +11,61590.002,17.013050114521022,0.04399999976158142,0.4279341627624582,0.023699999997916166,0.04399999976158142,16.585115951758564,61589.9783,61589.9783 +12,61590.24357,17.091348238945375,0.02199999988079071,0.5062322871868119,0.26526999999623513,0.02199999988079071,16.585115951758564,61589.9783,61589.9783 +13,61616.20243,17.44818378972448,0.032999999821186066,0.5004836109665014,0.19327999999950407,0.032999999821186066,16.94770017875798,61616.00915,61616.00915 +14,61648.07405,17.627385699986323,0.054999999701976776,0.4917527634000258,0.024980000001960434,0.054999999701976776,17.135632936586298,61648.04907,61648.04907 diff --git a/tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv b/tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv new file mode 100644 index 0000000..d55b95f --- /dev/null +++ b/tests/data/df_6098332225018_g_r_app_3_N_ref_3.csv @@ -0,0 +1,16 @@ +,midPointMjdTai,reduced_mag,magErr,g-r,delta_t_g-r,g-rErr,reduced_mag_r,midPointMjdTai1_r,midPointMjdTai2_r +0,61485.3588,17.82546893645328,0.13199999928474426,0.6646031133749091,7.961480000005395,0.1387787184136571,17.16086582307837,61469.38552,61477.39732 +1,61500.36601,17.548929910716545,0.052000001072883606,0.3998204721333032,3.041850000001432,0.061730424176583504,17.149109438583242,61477.39732,61497.32416 +2,61504.33885,17.545916408017177,0.05700000002980232,0.41306100220004893,0.024510000002919696,0.06719475757019573,17.13285540581713,61503.35033,61504.31434 +3,61524.2819,17.60109262744087,0.06199999898672104,0.4682372216237418,19.967560000004596,0.07148521046953224,17.13285540581713,61503.35033,61504.31434 +4,61525.31814,17.530157869823057,0.04699999839067459,0.41017516413599964,0.024610000000393484,0.047603929903510755,17.119982705687057,61504.31434,61525.29353 +5,61536.27685,17.474481419144603,0.039000000804662704,0.3696556612160329,4.973070000000007,0.0435304535193636,17.10482575792857,61524.3074,61531.30378 +6,61556.37692,17.262436542715754,0.03200000151991844,0.23760932081133745,0.024060000003373716,0.11609811521344783,17.024827221904417,61525.29353,61556.35286 +7,61558.20527,17.330018817677058,0.03200000151991844,0.30519159577264077,1.8524099999995087,0.11609811521344783,17.024827221904417,61525.29353,61556.35286 +8,61558.20674,17.302025752576462,0.029999999329447746,0.27719853067204525,1.8538800000023912,0.11556284964731135,17.024827221904417,61525.29353,61556.35286 +9,61589.25295,17.086183244116974,0.028999999165534973,0.483080272053364,1.968010000004142,0.09278157251016264,16.60310297206361,61562.31723,61587.28494 +10,61589.25519,17.035176049202192,0.032999999821186066,0.43207307713858256,1.9702500000057626,0.09410855558372905,16.60310297206361,61562.31723,61587.28494 +11,61590.002,17.013050114521022,0.04399999976158142,0.4132712746425149,0.023699999997916166,0.05221424730174361,16.599778839878507,61589.277,61589.9783 +12,61590.24357,17.091348238945375,0.02199999988079071,0.4915693990668686,0.26526999999623513,0.03569772593630112,16.599778839878507,61589.277,61589.9783 +13,61616.20243,17.44818378972448,0.032999999821186066,0.5299196193560611,0.19327999999950407,0.041117958024899405,16.91826417036842,61610.20336,61616.00915 +14,61648.07405,17.627385699986323,0.054999999701976776,0.5006338385771478,0.024980000001960434,0.08026750494410353,17.126751861409176,61620.16575,61648.04907 diff --git a/tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv b/tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv new file mode 100644 index 0000000..f083e76 --- /dev/null +++ b/tests/data/df_6098332225018_g_r_app_3_N_ref_5.csv @@ -0,0 +1,16 @@ +,midPointMjdTai,reduced_mag,magErr,g-r,delta_t_g-r,g-rErr,reduced_mag_r,midPointMjdTai1_r,midPointMjdTai2_r +0,61485.3588,17.82546893645328,0.13199999928474426,0.6646031133749091,7.961480000005395,0.1387787184136571,17.16086582307837,61469.38552,61477.39732 +1,61500.36601,17.548929910716545,0.052000001072883606,0.40214961118913806,3.041850000001432,0.06408282905927162,17.146780299527407,61469.38552,61497.32416 +2,61504.33885,17.545916408017177,0.05700000002980232,0.411249516593152,0.024510000002919696,0.06403779113904624,17.134666891424025,61497.32416,61504.31434 +3,61524.2819,17.60109262744087,0.06199999898672104,0.4664257360168449,19.967560000004596,0.0685261888982895,17.134666891424025,61497.32416,61504.31434 +4,61525.31814,17.530157869823057,0.04699999839067459,0.4033930271697166,0.024610000000393484,0.05524451285749118,17.12676484265334,61503.35033,61525.29353 +5,61536.27685,17.474481419144603,0.039000000804662704,0.3506598671390968,4.973070000000007,0.05091030826493725,17.123821552005506,61503.35078,61531.30378 +6,61556.37692,17.262436542715754,0.03200000151991844,0.20073869280633616,0.024060000003373716,0.10276020693938347,17.061697849909418,61504.31434,61556.35286 +7,61558.20527,17.330018817677058,0.03200000151991844,0.2683209677676395,1.8524099999995087,0.10276020693938347,17.061697849909418,61504.31434,61556.35286 +8,61558.20674,17.302025752576462,0.029999999329447746,0.24032790266704396,1.8538800000023912,0.10215507815432866,17.061697849909418,61504.31434,61556.35286 +9,61589.25295,17.086183244116974,0.028999999165534973,0.40181661190363016,1.968010000004142,0.12432445391335133,16.684366632213344,61558.23271,61587.28494 +10,61589.25519,17.035176049202192,0.032999999821186066,0.3508094169888487,1.9702500000057626,0.12531787533089708,16.684366632213344,61558.23271,61587.28494 +11,61590.002,17.013050114521022,0.04399999976158142,0.4366844446633422,0.023699999997916166,0.057457452706149006,16.57636566985768,61586.24199,61589.9783 +12,61590.24357,17.091348238945375,0.02199999988079071,0.5149825690876959,0.26526999999623513,0.04300417290467258,16.57636566985768,61586.24199,61589.9783 +13,61616.20243,17.44818378972448,0.032999999821186066,0.6651807942177115,0.19327999999950407,0.17001002587811542,16.78300299550677,61589.27924,61616.00915 +14,61648.07405,17.627385699986323,0.054999999701976776,0.5616368620560408,0.024980000001960434,0.10406682109366068,17.065748837930283,61616.22684,61648.04907 diff --git a/tests/data/testing_database.db b/tests/data/testing_database.db index 7a220e3..455b28b 100644 Binary files a/tests/data/testing_database.db and b/tests/data/testing_database.db differ