From 165b544b85a9760c4d11d73ee45d7ef8afd6dd23 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 16 Jul 2024 10:43:27 +0100 Subject: [PATCH 01/10] montecarlo phase fit nb --- notebooks/mc_phase_fit.ipynb | 279 +++++++++++++++++++++++++++++++++++ 1 file changed, 279 insertions(+) create mode 100644 notebooks/mc_phase_fit.ipynb diff --git a/notebooks/mc_phase_fit.ipynb b/notebooks/mc_phase_fit.ipynb new file mode 100644 index 0000000..8a04d8b --- /dev/null +++ b/notebooks/mc_phase_fit.ipynb @@ -0,0 +1,279 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "3f265be4", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4eb75bb", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2154f51", + "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.utilities.plotting_utilities import plot_errorbar\n", + "\n", + "# from adler.utilities.science_utilities import 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 astropy.modeling.fitting import SLSQPLSQFitter" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdd1f0f2", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "# ssoid = \"8268570668335894776\" # NEO\n", + "ssoid = \"6098332225018\" # good MBA test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "472609a0", + "metadata": {}, + "outputs": [], + "source": [ + "# 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": "e2ae7c0b", + "metadata": {}, + "outputs": [], + "source": [ + "resample = 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b0c1412", + "metadata": {}, + "outputs": [], + "source": [ + "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.array(getattr(obs, \"phaseAngle\")) * u.deg,)\n", + " mag_err = np.array(getattr(obs, \"magErr\")) * u.mag\n", + " red_mag = np.array(getattr(obs, \"reduced_mag\")) * u.mag\n", + "\n", + " pc_fit = pc.FitModel(\n", + " np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n", + " red_mag,\n", + " mag_err,\n", + " )\n", + " pc = pc.InitModelSbpy(pc_fit)\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", + " alpha_fit = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg\n", + " ax1.plot(\n", + " alpha_fit.value,\n", + " pc.ReducedMag(alpha_fit).value,\n", + " label=\"{}: H={:.2f}, G12={:.2f}\".format(filt, pc.H, pc.phase_parameter_1),\n", + " )\n", + " ax1.legend()\n", + "\n", + " mc_models = []\n", + " for i in range(resample):\n", + " _red_mag = np.random.normal(loc=np.array(red_mag), scale=np.array(mag_err)) * u.mag\n", + " ax1.scatter(alpha[0].value, _red_mag.value, s=1, c=\"k\", alpha=0.1, zorder=5)\n", + "\n", + " _pc_fit = pc.FitModel(\n", + " np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n", + " _red_mag,\n", + " )\n", + " _pc = pc.InitModelSbpy(_pc_fit)\n", + " ax1.plot(\n", + " alpha_fit.value,\n", + " _pc.ReducedMag(alpha_fit).value,\n", + " c=\"k\",\n", + " alpha=0.1,\n", + " zorder=0,\n", + " label=\"{}: H={:.2f}, G12={:.2f}\".format(filt, _pc.H, _pc.phase_parameter_1),\n", + " )\n", + "\n", + " mc_models.append(_pc)\n", + "\n", + " break\n", + "\n", + "ax1.invert_yaxis()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "919f0e21", + "metadata": {}, + "outputs": [], + "source": [ + "mag_err" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea39b9b0", + "metadata": {}, + "outputs": [], + "source": [ + "mc_models[0].H, mc_models[0].phase_parameter_1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc96d3d6", + "metadata": {}, + "outputs": [], + "source": [ + "abs_mag = np.array([x.H.value for x in mc_models])\n", + "phase_param = np.array([x.phase_parameter_1 for x in mc_models])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "040867f2", + "metadata": {}, + "outputs": [], + "source": [ + "phase_param" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ecaff7e", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 2)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "ax2 = plt.subplot(gs[0, 1])\n", + "\n", + "for x, ax, ax_lab in zip([abs_mag, phase_param], [ax1, ax2], [\"H\", \"phase_parameter_1\"]):\n", + " ax.hist(x, bins=\"auto\", histtype=\"step\")\n", + " mean = np.mean(x)\n", + " std = np.std(x)\n", + " print(mean, std)\n", + " ax.axvline(mean)\n", + " ax.axvspan(mean - std, mean + std, alpha=0.3)\n", + " ax.set_xlabel(ax_lab)\n", + "\n", + "for x, x_err, ax in zip([\"H\", \"phase_parameter_1\"], [\"H_err\", \"phase_parameter_1_err\"], [ax1, ax2]):\n", + " x = getattr(pc, x)\n", + " x_err = getattr(pc, x_err)\n", + "\n", + " if hasattr(x, \"unit\"):\n", + " x = x.value\n", + "\n", + " ax.axvline(x, color=\"k\")\n", + " ax.axvspan(x - x_err, x + x_err, alpha=0.3, color=\"k\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "219b3e89", + "metadata": {}, + "outputs": [], + "source": [ + "pc.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa2bb33a", + "metadata": {}, + "outputs": [], + "source": [ + "_pc.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6a0cbce", + "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 +} From 7e04166ad96bdd0d257a06b032abcab8971523d4 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 16 Jul 2024 09:47:49 +0000 Subject: [PATCH 02/10] update nb --- notebooks/mc_phase_fit.ipynb | 166 +++++++++++++++++++++++++++++++---- 1 file changed, 149 insertions(+), 17 deletions(-) diff --git a/notebooks/mc_phase_fit.ipynb b/notebooks/mc_phase_fit.ipynb index 8a04d8b..1eda8a7 100644 --- a/notebooks/mc_phase_fit.ipynb +++ b/notebooks/mc_phase_fit.ipynb @@ -4,7 +4,15 @@ "cell_type": "code", "execution_count": null, "id": "3f265be4", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:06.759338Z", + "iopub.status.busy": "2024-07-16T09:46:06.759037Z", + "iopub.status.idle": "2024-07-16T09:46:07.686648Z", + "shell.execute_reply": "2024-07-16T09:46:07.685945Z", + "shell.execute_reply.started": "2024-07-16T09:46:06.759319Z" + } + }, "outputs": [], "source": [ "%matplotlib notebook" @@ -14,7 +22,15 @@ "cell_type": "code", "execution_count": null, "id": "e4eb75bb", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:07.688067Z", + "iopub.status.busy": "2024-07-16T09:46:07.687672Z", + "iopub.status.idle": "2024-07-16T09:46:07.692332Z", + "shell.execute_reply": "2024-07-16T09:46:07.691639Z", + "shell.execute_reply.started": "2024-07-16T09:46:07.688043Z" + } + }, "outputs": [], "source": [ "%matplotlib inline" @@ -24,7 +40,15 @@ "cell_type": "code", "execution_count": null, "id": "e2154f51", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:07.693397Z", + "iopub.status.busy": "2024-07-16T09:46:07.693191Z", + "iopub.status.idle": "2024-07-16T09:46:10.794076Z", + "shell.execute_reply": "2024-07-16T09:46:10.793249Z", + "shell.execute_reply.started": "2024-07-16T09:46:07.693382Z" + } + }, "outputs": [], "source": [ "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", @@ -46,31 +70,57 @@ "cell_type": "code", "execution_count": null, "id": "fdd1f0f2", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:10.797096Z", + "iopub.status.busy": "2024-07-16T09:46:10.796369Z", + "iopub.status.idle": "2024-07-16T09:46:10.800291Z", + "shell.execute_reply": "2024-07-16T09:46:10.799642Z", + "shell.execute_reply.started": "2024-07-16T09:46:10.797072Z" + } + }, "outputs": [], "source": [ "# ssObjectId of object to analyse\n", "# ssoid = \"8268570668335894776\" # NEO\n", - "ssoid = \"6098332225018\" # good MBA test object" + "# ssoid = \"6098332225018\" # good MBA test object\n", + "ssoid = \"5334524274243260416\"" ] }, { "cell_type": "code", "execution_count": null, "id": "472609a0", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:10.801360Z", + "iopub.status.busy": "2024-07-16T09:46:10.801151Z", + "iopub.status.idle": "2024-07-16T09:46:13.171505Z", + "shell.execute_reply": "2024-07-16T09:46:13.170567Z", + "shell.execute_reply.started": "2024-07-16T09:46:10.801344Z" + } + }, "outputs": [], "source": [ "# 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)" + "# fname = \"/Users/jrobinson/lsst-adler/notebooks/gen_test_data/adler_demo_testing_database.db\"\n", + "# planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)\n", + "planetoid = AdlerPlanetoid.construct_from_RSP(ssoid)" ] }, { "cell_type": "code", "execution_count": null, "id": "e2ae7c0b", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:13.173102Z", + "iopub.status.busy": "2024-07-16T09:46:13.172834Z", + "iopub.status.idle": "2024-07-16T09:46:13.177158Z", + "shell.execute_reply": "2024-07-16T09:46:13.176249Z", + "shell.execute_reply.started": "2024-07-16T09:46:13.173083Z" + } + }, "outputs": [], "source": [ "resample = 100" @@ -80,7 +130,15 @@ "cell_type": "code", "execution_count": null, "id": "5b0c1412", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:13.178331Z", + "iopub.status.busy": "2024-07-16T09:46:13.178106Z", + "iopub.status.idle": "2024-07-16T09:46:15.403027Z", + "shell.execute_reply": "2024-07-16T09:46:15.401842Z", + "shell.execute_reply.started": "2024-07-16T09:46:13.178314Z" + } + }, "outputs": [], "source": [ "fig = plt.figure()\n", @@ -155,17 +213,51 @@ "cell_type": "code", "execution_count": null, "id": "919f0e21", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.405030Z", + "iopub.status.busy": "2024-07-16T09:46:15.404641Z", + "iopub.status.idle": "2024-07-16T09:46:15.413186Z", + "shell.execute_reply": "2024-07-16T09:46:15.412407Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.405000Z" + } + }, "outputs": [], "source": [ "mag_err" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "321547a5-28ff-416a-b370-93e4f55d4370", + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:50.718468Z", + "iopub.status.busy": "2024-07-16T09:46:50.718060Z", + "iopub.status.idle": "2024-07-16T09:46:50.725967Z", + "shell.execute_reply": "2024-07-16T09:46:50.725006Z", + "shell.execute_reply.started": "2024-07-16T09:46:50.718441Z" + } + }, + "outputs": [], + "source": [ + "np.median(mag_err)" + ] + }, { "cell_type": "code", "execution_count": null, "id": "ea39b9b0", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.414457Z", + "iopub.status.busy": "2024-07-16T09:46:15.414220Z", + "iopub.status.idle": "2024-07-16T09:46:15.429426Z", + "shell.execute_reply": "2024-07-16T09:46:15.428699Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.414441Z" + } + }, "outputs": [], "source": [ "mc_models[0].H, mc_models[0].phase_parameter_1" @@ -175,7 +267,15 @@ "cell_type": "code", "execution_count": null, "id": "dc96d3d6", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.430747Z", + "iopub.status.busy": "2024-07-16T09:46:15.430475Z", + "iopub.status.idle": "2024-07-16T09:46:15.445087Z", + "shell.execute_reply": "2024-07-16T09:46:15.444131Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.430728Z" + } + }, "outputs": [], "source": [ "abs_mag = np.array([x.H.value for x in mc_models])\n", @@ -186,7 +286,15 @@ "cell_type": "code", "execution_count": null, "id": "040867f2", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.446446Z", + "iopub.status.busy": "2024-07-16T09:46:15.446193Z", + "iopub.status.idle": "2024-07-16T09:46:15.471748Z", + "shell.execute_reply": "2024-07-16T09:46:15.470912Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.446429Z" + } + }, "outputs": [], "source": [ "phase_param" @@ -196,7 +304,15 @@ "cell_type": "code", "execution_count": null, "id": "0ecaff7e", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.473021Z", + "iopub.status.busy": "2024-07-16T09:46:15.472795Z", + "iopub.status.idle": "2024-07-16T09:46:15.821399Z", + "shell.execute_reply": "2024-07-16T09:46:15.820487Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.473003Z" + } + }, "outputs": [], "source": [ "fig = plt.figure()\n", @@ -230,7 +346,15 @@ "cell_type": "code", "execution_count": null, "id": "219b3e89", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.822907Z", + "iopub.status.busy": "2024-07-16T09:46:15.822565Z", + "iopub.status.idle": "2024-07-16T09:46:15.829040Z", + "shell.execute_reply": "2024-07-16T09:46:15.828114Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.822887Z" + } + }, "outputs": [], "source": [ "pc.__dict__" @@ -240,7 +364,15 @@ "cell_type": "code", "execution_count": null, "id": "aa2bb33a", - "metadata": {}, + "metadata": { + "execution": { + "iopub.execute_input": "2024-07-16T09:46:15.830321Z", + "iopub.status.busy": "2024-07-16T09:46:15.830089Z", + "iopub.status.idle": "2024-07-16T09:46:15.884244Z", + "shell.execute_reply": "2024-07-16T09:46:15.883201Z", + "shell.execute_reply.started": "2024-07-16T09:46:15.830304Z" + } + }, "outputs": [], "source": [ "_pc.__dict__" From f43a741acaa0b920aa5625b1058d1dffbaad4995 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 16 Jul 2024 10:49:45 +0100 Subject: [PATCH 03/10] intitial mc phase func --- src/adler/science/PhaseCurve.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index a40d1cb..a3a14ce 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -233,7 +233,7 @@ def ReducedMag(self, phase_angle): return self.model_function(phase_angle) - def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): + def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None, resample=None): """Fit the phasecurve model parameters to observations. starts with a phase curve model as an initial guess for parameters. fits model to phase angle and reduced magnitude. @@ -258,6 +258,10 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): Select a fitting function from astropy.modeling.fitting, defaults to astropy.modeling.fitting.LevMarLSQFitter. N.B. that LevMarLSQFitter cannot handle inequality constraints for the HG1G2 model, use something like SLSQPLSQFitter from astropy.modeling.fitting (does not return covariance matrix!). + resample : int + Optional - if passed this forces a monte carlo resampling of data points within their uncertainties. + This the number of times to resample and fit the phase curve. + The phase curve parameter uncertainties are determined from the mean and std of the fitted values. Returns ---------- @@ -276,7 +280,7 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): model_fit = fitter(self.model_function, phase_angle, reduced_mag) # Add fitted uncertainties as an additional attribute within the sbpy object - if "param_cov" in fitter.fit_info: + if ("param_cov" in fitter.fit_info) and (resample is None): # get the covariance matrix from the fit covariance = fitter.fit_info["param_cov"] if covariance is not None: @@ -287,13 +291,19 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None): fit_mask = ~np.array([getattr(model_fit, x).fixed for x in param_names]) for i, x in enumerate(param_names[fit_mask]): setattr(model_fit, "{}_err".format(x), fit_errs[i]) + # TODO: return uncertainties with units if units are passed # else: ### TODO log covariance is None error here - ### TODO - # else: - # log lack of uncertainties for fitter - # run an MCMC estimate of uncertainty? + elif (resample is not None) and (mag_err is not None): + # run an MC estimate of uncertainty? + + _reduced_mag = np.random.normal(loc=np.array(reduced_mag), scale=np.array(mag_err)) * u.mag + _model_fit = fitter(self.model_function, phase_angle, _reduced_mag) + + else: + # log lack of uncertainties for fitter + print("no phase curve parameter uncertainties") ### if overwrite_model: # add an overwrite option? # redo __init__ with the new fitted parameters From a6f1c723e12059deb696edc611bddc0a4bd3f6a9 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 16 Jul 2024 12:26:11 +0100 Subject: [PATCH 04/10] add MC phase fitter --- src/adler/science/PhaseCurve.py | 34 ++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index a3a14ce..1d84197 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -261,7 +261,7 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None, resample resample : int Optional - if passed this forces a monte carlo resampling of data points within their uncertainties. This the number of times to resample and fit the phase curve. - The phase curve parameter uncertainties are determined from the mean and std of the fitted values. + The phase curve parameter value and uncertainties are determined from the mean and std of the fitted values respectively. Returns ---------- @@ -291,19 +291,39 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None, resample fit_mask = ~np.array([getattr(model_fit, x).fixed for x in param_names]) for i, x in enumerate(param_names[fit_mask]): setattr(model_fit, "{}_err".format(x), fit_errs[i]) - # TODO: return uncertainties with units if units are passed + # TODO: return uncertainties with units if units are passed - see MC resample code below # else: - ### TODO log covariance is None error here + ### TODO log covariance is None error here - no uncertainties + # run an MC reasmple fit to estimate parameter value and uncertainty elif (resample is not None) and (mag_err is not None): - # run an MC estimate of uncertainty? + mc_models = [] # list to store MC model fits + for i in range(resample): # TODO: try optimise/parallelise this loop? + _reduced_mag = np.random.normal(loc=np.array(reduced_mag), scale=np.array(mag_err)) * u.mag + _model_fit = fitter(self.model_function, phase_angle, _reduced_mag) + mc_models.append(_model_fit) + + # Update the model_fit parameters with the MC values + param_names = np.array(model_fit.param_names) + # update only the uncertainties for parameters used in the fit + fit_mask = ~np.array([getattr(model_fit, x).fixed for x in param_names]) + for i, x in enumerate(param_names[fit_mask]): + # check if the parameter has units and then get array of the MC model values + m = mc_models[0] + p = getattr(m, x) + if hasattr(p, "unit"): + fit_vals = np.array([getattr(m, x).value for m in mc_models]) * p.unit + else: + fit_vals = np.array([getattr(m, x) for m in mc_models]) - _reduced_mag = np.random.normal(loc=np.array(reduced_mag), scale=np.array(mag_err)) * u.mag - _model_fit = fitter(self.model_function, phase_angle, _reduced_mag) + # set the parameter value as the mean of the MC values + setattr(model_fit, "{}".format(x), np.mean(fit_vals)) + # set the parameter uncertainty as the std of the MC values + setattr(model_fit, "{}_err".format(x), np.std(fit_vals)) else: # log lack of uncertainties for fitter - print("no phase curve parameter uncertainties") + print("no phase curve parameter uncertainties calculated") ### if overwrite_model: # add an overwrite option? # redo __init__ with the new fitted parameters From 6e7091373a6ed81e97fdaebdd57ed2451a49ab46 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 16 Jul 2024 12:41:59 +0100 Subject: [PATCH 05/10] keep units with uncertainties if available --- src/adler/science/PhaseCurve.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index 1d84197..75efeea 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -199,7 +199,9 @@ def InitModelSbpy(self, model_sbpy): x.unit == "" ): # if there are no units (or weird blank units?) get just the value x = x.value - else: + if hasattr( + x, "quantity" + ): # catch any sbpy parameters returned as astropy.modeling.parameters.Parameter x = x.quantity # look up the correct adler parameter name (accounting for additional uncertainty, "_err", parameters) if p.endswith("_err"): # assumes the uncertainty parameter always ends in "_err" @@ -290,7 +292,11 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None, resample param_names = np.array(model_fit.param_names) fit_mask = ~np.array([getattr(model_fit, x).fixed for x in param_names]) for i, x in enumerate(param_names[fit_mask]): - setattr(model_fit, "{}_err".format(x), fit_errs[i]) + p = getattr(model_fit, x) + if hasattr(p, "unit") and (p.unit is not None): + setattr(model_fit, "{}_err".format(x), fit_errs[i] * p.unit) + else: + setattr(model_fit, "{}_err".format(x), fit_errs[i]) # TODO: return uncertainties with units if units are passed - see MC resample code below # else: ### TODO log covariance is None error here - no uncertainties From f6ede626ebd90f6c78129c803c00aa989ccb6983 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 16 Jul 2024 12:42:50 +0100 Subject: [PATCH 06/10] keep units with uncertainties if available --- src/adler/science/PhaseCurve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index 75efeea..9756090 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -317,7 +317,7 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None, resample # check if the parameter has units and then get array of the MC model values m = mc_models[0] p = getattr(m, x) - if hasattr(p, "unit"): + if hasattr(p, "unit") and (p.unit is not None): fit_vals = np.array([getattr(m, x).value for m in mc_models]) * p.unit else: fit_vals = np.array([getattr(m, x) for m in mc_models]) From 8519635d41c40aef780d3f108b57eb22d64584be Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 30 Jul 2024 10:11:12 +0100 Subject: [PATCH 07/10] update nb --- notebooks/mc_phase_fit.ipynb | 221 +++++++++++++---------------------- 1 file changed, 81 insertions(+), 140 deletions(-) diff --git a/notebooks/mc_phase_fit.ipynb b/notebooks/mc_phase_fit.ipynb index 1eda8a7..5f9d806 100644 --- a/notebooks/mc_phase_fit.ipynb +++ b/notebooks/mc_phase_fit.ipynb @@ -4,15 +4,7 @@ "cell_type": "code", "execution_count": null, "id": "3f265be4", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:06.759338Z", - "iopub.status.busy": "2024-07-16T09:46:06.759037Z", - "iopub.status.idle": "2024-07-16T09:46:07.686648Z", - "shell.execute_reply": "2024-07-16T09:46:07.685945Z", - "shell.execute_reply.started": "2024-07-16T09:46:06.759319Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook" @@ -22,15 +14,7 @@ "cell_type": "code", "execution_count": null, "id": "e4eb75bb", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:07.688067Z", - "iopub.status.busy": "2024-07-16T09:46:07.687672Z", - "iopub.status.idle": "2024-07-16T09:46:07.692332Z", - "shell.execute_reply": "2024-07-16T09:46:07.691639Z", - "shell.execute_reply.started": "2024-07-16T09:46:07.688043Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" @@ -40,15 +24,7 @@ "cell_type": "code", "execution_count": null, "id": "e2154f51", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:07.693397Z", - "iopub.status.busy": "2024-07-16T09:46:07.693191Z", - "iopub.status.idle": "2024-07-16T09:46:10.794076Z", - "shell.execute_reply": "2024-07-16T09:46:10.793249Z", - "shell.execute_reply.started": "2024-07-16T09:46:07.693382Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", @@ -70,57 +46,33 @@ "cell_type": "code", "execution_count": null, "id": "fdd1f0f2", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:10.797096Z", - "iopub.status.busy": "2024-07-16T09:46:10.796369Z", - "iopub.status.idle": "2024-07-16T09:46:10.800291Z", - "shell.execute_reply": "2024-07-16T09:46:10.799642Z", - "shell.execute_reply.started": "2024-07-16T09:46:10.797072Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "# ssObjectId of object to analyse\n", "# ssoid = \"8268570668335894776\" # NEO\n", - "# ssoid = \"6098332225018\" # good MBA test object\n", - "ssoid = \"5334524274243260416\"" + "ssoid = \"6098332225018\" # good MBA test object\n", + "# ssoid = \"5334524274243260416\"" ] }, { "cell_type": "code", "execution_count": null, "id": "472609a0", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:10.801360Z", - "iopub.status.busy": "2024-07-16T09:46:10.801151Z", - "iopub.status.idle": "2024-07-16T09:46:13.171505Z", - "shell.execute_reply": "2024-07-16T09:46:13.170567Z", - "shell.execute_reply.started": "2024-07-16T09:46:10.801344Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "# 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)\n", - "planetoid = AdlerPlanetoid.construct_from_RSP(ssoid)" + "fname = \"/Users/jrobinson/lsst-adler/notebooks/gen_test_data/adler_demo_testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)\n", + "# planetoid = AdlerPlanetoid.construct_from_RSP(ssoid)" ] }, { "cell_type": "code", "execution_count": null, "id": "e2ae7c0b", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:13.173102Z", - "iopub.status.busy": "2024-07-16T09:46:13.172834Z", - "iopub.status.idle": "2024-07-16T09:46:13.177158Z", - "shell.execute_reply": "2024-07-16T09:46:13.176249Z", - "shell.execute_reply.started": "2024-07-16T09:46:13.173083Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "resample = 100" @@ -130,15 +82,7 @@ "cell_type": "code", "execution_count": null, "id": "5b0c1412", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:13.178331Z", - "iopub.status.busy": "2024-07-16T09:46:13.178106Z", - "iopub.status.idle": "2024-07-16T09:46:15.403027Z", - "shell.execute_reply": "2024-07-16T09:46:15.401842Z", - "shell.execute_reply.started": "2024-07-16T09:46:13.178314Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", @@ -209,19 +153,71 @@ "plt.show()" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d0ef23a-ea32-4e32-95ed-a344c4f2760b", + "metadata": {}, + "outputs": [], + "source": [ + "pc_fit.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8fa4202e-e7f4-4438-a41e-dc7582312bca", + "metadata": {}, + "outputs": [], + "source": [ + "_pc_fit.param_names" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18e04044-8007-4b41-85d5-a9b1ad720223", + "metadata": {}, + "outputs": [], + "source": [ + "mc_models[0].__dict__.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f43eb646-cf4e-4103-a9d4-7f8b6c8d2738", + "metadata": {}, + "outputs": [], + "source": [ + "_pc_fit = pc.FitModel(np.array(getattr(obs, \"phaseAngle\")) * u.deg, red_mag, mag_err, resample=resample)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "871dc090-7c0e-41d7-bbc9-9a0ec0232d47", + "metadata": {}, + "outputs": [], + "source": [ + "_pc_fit.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71cce253-96f2-44a9-923b-9d4af3c50db4", + "metadata": {}, + "outputs": [], + "source": [ + "pc_fit.__dict__" + ] + }, { "cell_type": "code", "execution_count": null, "id": "919f0e21", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.405030Z", - "iopub.status.busy": "2024-07-16T09:46:15.404641Z", - "iopub.status.idle": "2024-07-16T09:46:15.413186Z", - "shell.execute_reply": "2024-07-16T09:46:15.412407Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.405000Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "mag_err" @@ -231,15 +227,7 @@ "cell_type": "code", "execution_count": null, "id": "321547a5-28ff-416a-b370-93e4f55d4370", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:50.718468Z", - "iopub.status.busy": "2024-07-16T09:46:50.718060Z", - "iopub.status.idle": "2024-07-16T09:46:50.725967Z", - "shell.execute_reply": "2024-07-16T09:46:50.725006Z", - "shell.execute_reply.started": "2024-07-16T09:46:50.718441Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "np.median(mag_err)" @@ -249,15 +237,7 @@ "cell_type": "code", "execution_count": null, "id": "ea39b9b0", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.414457Z", - "iopub.status.busy": "2024-07-16T09:46:15.414220Z", - "iopub.status.idle": "2024-07-16T09:46:15.429426Z", - "shell.execute_reply": "2024-07-16T09:46:15.428699Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.414441Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "mc_models[0].H, mc_models[0].phase_parameter_1" @@ -267,15 +247,7 @@ "cell_type": "code", "execution_count": null, "id": "dc96d3d6", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.430747Z", - "iopub.status.busy": "2024-07-16T09:46:15.430475Z", - "iopub.status.idle": "2024-07-16T09:46:15.445087Z", - "shell.execute_reply": "2024-07-16T09:46:15.444131Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.430728Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "abs_mag = np.array([x.H.value for x in mc_models])\n", @@ -286,15 +258,7 @@ "cell_type": "code", "execution_count": null, "id": "040867f2", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.446446Z", - "iopub.status.busy": "2024-07-16T09:46:15.446193Z", - "iopub.status.idle": "2024-07-16T09:46:15.471748Z", - "shell.execute_reply": "2024-07-16T09:46:15.470912Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.446429Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "phase_param" @@ -304,15 +268,7 @@ "cell_type": "code", "execution_count": null, "id": "0ecaff7e", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.473021Z", - "iopub.status.busy": "2024-07-16T09:46:15.472795Z", - "iopub.status.idle": "2024-07-16T09:46:15.821399Z", - "shell.execute_reply": "2024-07-16T09:46:15.820487Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.473003Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", @@ -335,6 +291,7 @@ "\n", " if hasattr(x, \"unit\"):\n", " x = x.value\n", + " x_err = x_err.value\n", "\n", " ax.axvline(x, color=\"k\")\n", " ax.axvspan(x - x_err, x + x_err, alpha=0.3, color=\"k\")\n", @@ -346,15 +303,7 @@ "cell_type": "code", "execution_count": null, "id": "219b3e89", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.822907Z", - "iopub.status.busy": "2024-07-16T09:46:15.822565Z", - "iopub.status.idle": "2024-07-16T09:46:15.829040Z", - "shell.execute_reply": "2024-07-16T09:46:15.828114Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.822887Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "pc.__dict__" @@ -364,15 +313,7 @@ "cell_type": "code", "execution_count": null, "id": "aa2bb33a", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-16T09:46:15.830321Z", - "iopub.status.busy": "2024-07-16T09:46:15.830089Z", - "iopub.status.idle": "2024-07-16T09:46:15.884244Z", - "shell.execute_reply": "2024-07-16T09:46:15.883201Z", - "shell.execute_reply.started": "2024-07-16T09:46:15.830304Z" - } - }, + "metadata": {}, "outputs": [], "source": [ "_pc.__dict__" From b8b9964a6f659a4cef457b57802680a80aa3b545 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 3 Dec 2024 18:22:24 +0000 Subject: [PATCH 08/10] testing nb --- ...e_fit.ipynb => mc_phase_fit_testing.ipynb} | 79 +++++++++++++++---- 1 file changed, 65 insertions(+), 14 deletions(-) rename notebooks/{mc_phase_fit.ipynb => mc_phase_fit_testing.ipynb} (81%) diff --git a/notebooks/mc_phase_fit.ipynb b/notebooks/mc_phase_fit_testing.ipynb similarity index 81% rename from notebooks/mc_phase_fit.ipynb rename to notebooks/mc_phase_fit_testing.ipynb index 5f9d806..78584e3 100644 --- a/notebooks/mc_phase_fit.ipynb +++ b/notebooks/mc_phase_fit_testing.ipynb @@ -27,8 +27,8 @@ "metadata": {}, "outputs": [], "source": [ - "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", - "from adler.dataclasses.AdlerData import AdlerData\n", + "from adler.objectdata.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.objectdata.AdlerData import AdlerData\n", "from adler.science.PhaseCurve import PhaseCurve\n", "from adler.utilities.plotting_utilities import plot_errorbar\n", "\n", @@ -190,7 +190,8 @@ "metadata": {}, "outputs": [], "source": [ - "_pc_fit = pc.FitModel(np.array(getattr(obs, \"phaseAngle\")) * u.deg, red_mag, mag_err, resample=resample)" + "pc_fit_resamp = pc.FitModel(np.array(getattr(obs, \"phaseAngle\")) * u.deg, red_mag, mag_err, resample=resample)\n", + "pc_fit_resamp = PhaseCurve().InitModelSbpy(pc_fit_resamp)" ] }, { @@ -200,7 +201,17 @@ "metadata": {}, "outputs": [], "source": [ - "_pc_fit.__dict__" + "pc_fit_resamp.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3867d1de", + "metadata": {}, + "outputs": [], + "source": [ + "pc_fit_resamp.H, pc_fit_resamp.phase_parameter_1" ] }, { @@ -213,6 +224,16 @@ "pc_fit.__dict__" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f05c9eb", + "metadata": {}, + "outputs": [], + "source": [ + "pc.__dict__" + ] + }, { "cell_type": "code", "execution_count": null, @@ -264,6 +285,26 @@ "phase_param" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b21f899", + "metadata": {}, + "outputs": [], + "source": [ + "pc_fit.H, pc_fit.H_err" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1043bc2e", + "metadata": {}, + "outputs": [], + "source": [ + "pc.H, pc.H_err" + ] + }, { "cell_type": "code", "execution_count": null, @@ -276,25 +317,35 @@ "ax1 = plt.subplot(gs[0, 0])\n", "ax2 = plt.subplot(gs[0, 1])\n", "\n", + "label = \"resample test\"\n", "for x, ax, ax_lab in zip([abs_mag, phase_param], [ax1, ax2], [\"H\", \"phase_parameter_1\"]):\n", " ax.hist(x, bins=\"auto\", histtype=\"step\")\n", " mean = np.mean(x)\n", " std = np.std(x)\n", - " print(mean, std)\n", - " ax.axvline(mean)\n", + " print(label, ax_lab, mean, std)\n", + " ax.axvline(mean, label=label)\n", " ax.axvspan(mean - std, mean + std, alpha=0.3)\n", " ax.set_xlabel(ax_lab)\n", "\n", - "for x, x_err, ax in zip([\"H\", \"phase_parameter_1\"], [\"H_err\", \"phase_parameter_1_err\"], [ax1, ax2]):\n", - " x = getattr(pc, x)\n", - " x_err = getattr(pc, x_err)\n", + "for _pc, col, label in zip([pc, pc_fit_resamp], [\"k\", \"r\"], [\"single fit\", \"resample\"]):\n", + " for x_name, x_err_name, ax in zip(\n", + " [\"H\", \"phase_parameter_1\"], [\"H_err\", \"phase_parameter_1_err\"], [ax1, ax2]\n", + " ):\n", + " x = getattr(_pc, x_name)\n", + " x_err = getattr(_pc, x_err_name)\n", + "\n", + " if hasattr(x, \"unit\"):\n", + " x = x.value\n", + " if hasattr(x_err, \"unit\"):\n", + " x_err = x_err.value\n", + "\n", + " print(label, x_name, x, x_err)\n", "\n", - " if hasattr(x, \"unit\"):\n", - " x = x.value\n", - " x_err = x_err.value\n", + " ax.axvline(x, color=col, label=label)\n", + " ax.axvspan(x - x_err, x + x_err, alpha=0.3, color=col)\n", "\n", - " ax.axvline(x, color=\"k\")\n", - " ax.axvspan(x - x_err, x + x_err, alpha=0.3, color=\"k\")\n", + "ax1.legend()\n", + "ax2.legend()\n", "\n", "plt.show()" ] From f88f64fcda66fb60b3454458d4b97489a2bc0ca1 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 3 Dec 2024 20:48:17 +0000 Subject: [PATCH 09/10] more PhaseCurve tests --- src/adler/science/PhaseCurve.py | 11 ++- tests/adler/science/test_PhaseCurve.py | 113 ++++++++++++++++++++++++- 2 files changed, 119 insertions(+), 5 deletions(-) diff --git a/src/adler/science/PhaseCurve.py b/src/adler/science/PhaseCurve.py index ec3fd7e..1187634 100644 --- a/src/adler/science/PhaseCurve.py +++ b/src/adler/science/PhaseCurve.py @@ -160,6 +160,10 @@ def InitModelDict(self, model_dict): # clean the input dictionary del_keys = [] + # make sure that there is not an existing model_function parameter + if "model_function" in model_dict: + del_keys.append("model_function") + # remove keys that aren't PhaseCurve parameters for key, value in model_dict.items(): if not hasattr(self, key): del_keys.append(key) @@ -361,8 +365,13 @@ def FitModel(self, phase_angle, reduced_mag, mag_err=None, fitter=None, resample # else: # setattr(model_fit, "{}_err".format(x), fit_errs[i]) for i, x in enumerate(param_names[fit_mask]): - setattr(model_fit, "{}_err".format(x), fit_errs[i]) + # setattr(model_fit, "{}_err".format(x), fit_errs[i]) # TODO: return uncertainties with units if units are passed - see MC resample code below + p = getattr(model_fit, x) + if hasattr(p, "unit") and (p.unit is not None): + setattr(model_fit, "{}_err".format(x), fit_errs[i] * p.unit) + else: + setattr(model_fit, "{}_err".format(x), fit_errs[i]) # else: ### TODO log covariance is None error here - no uncertainties diff --git a/tests/adler/science/test_PhaseCurve.py b/tests/adler/science/test_PhaseCurve.py index b3d8fd7..5c2c3c9 100644 --- a/tests/adler/science/test_PhaseCurve.py +++ b/tests/adler/science/test_PhaseCurve.py @@ -1,8 +1,11 @@ -from adler.science.PhaseCurve import PhaseCurve -from numpy.testing import assert_array_equal -import pytest import numpy as np import astropy.units as u +from numpy.testing import assert_array_equal, assert_almost_equal +import pytest + +from adler.science.PhaseCurve import PhaseCurve, ADLER_SBPY_DICT +from adler.utilities.tests_utilities import get_test_data_filepath +from adler.objectdata.AdlerPlanetoid import AdlerPlanetoid def test_PhaseCurve_init(): @@ -145,4 +148,106 @@ def test_PhaseCurve_FitModel_HG_bounds(): assert pc2.phase_parameter_1_err is not None -# TODO: test absmag +def test_PhaseCurve_FitModel_resample(): + + np.random.seed(0) # set the seed to ensure reproducibility + resample = 100 + + # load a test object + ssoid = "6098332225018" # good MBA test object + filt = "r" + test_db_path = get_test_data_filepath("testing_database.db") + planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, test_db_path, filter_list=[filt]) + sso = planetoid.SSObject_in_filter(filt) + obs = planetoid.observations_in_filter(filt) + + # set up the initial phasecurve model + H = sso.H + G12 = sso.G12 + pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name="HG12_Pen16") + + # get the observations + alpha = np.array(getattr(obs, "phaseAngle")) * u.deg + mag_err = np.array(getattr(obs, "magErr")) * u.mag + red_mag = np.array(getattr(obs, "reduced_mag")) * u.mag + + # do a single fit + pc_fit = pc.FitModel(alpha, red_mag, mag_err) + pc_fit = PhaseCurve().InitModelSbpy(pc_fit) + print(pc_fit.__dict__) + + # use the resample function within FitModel + pc_fit_resamp = pc.FitModel(alpha, red_mag, mag_err, resample=resample) + pc_fit_resamp = PhaseCurve().InitModelSbpy(pc_fit_resamp) + print(pc_fit_resamp.__dict__) + + # check that the fits are different + tolerance = 0.001 + for x in ["H", "H_err", "phase_parameter_1", "phase_parameter_1_err"]: + x1 = getattr(pc_fit_resamp, x) + x2 = getattr(pc_fit, x) + + if hasattr(x1, "unit") and (x1.unit is not None): + x1 = x1.value + if hasattr(x2, "unit") and (x2.unit is not None): + x2 = x2.value + + print(x, np.abs(x1 - x2)) + assert np.abs(x1 - x2) > tolerance + + # the resampled fit should have larger uncertainties + if "err" in x: + assert x1 > x2 + + +def test_set_models(): + + for model in ["HG", "HG1G2", "HG12", "HG12_Pen16", "LinearPhaseFunc"]: + + # create a model + pc = PhaseCurve(model_name=model) + + # check which phase parameters the model should have + phase_params = ADLER_SBPY_DICT[pc.model_name] + for p in phase_params: + assert getattr(pc.model_function, phase_params[p]) + + +def test_ReturnInit(): + + # define a PhaseCurve object + pc = PhaseCurve(H=18.0, phase_parameter_1=0.15, model_name="HG") + + # convert the PhaseCurve object to dict + pc_dict = pc.ReturnModelDict() + print(pc_dict) + + assert pc_dict["H"] == 18.0 + assert pc_dict["phase_parameter_1"] == 0.15 + assert pc_dict["model_name"] == "HG" + + # create PhaseCurve from dict + print(type(pc_dict)) + pc2 = PhaseCurve().InitModelDict(pc_dict) + + assert pc2.H == 18.0 + assert pc2.phase_parameter_1 == 0.15 + assert pc2.model_name == "HG" + + # create sbpy model from PhaseCurve + pc_sbpy = pc.model_function + print(pc_sbpy) + + assert pc_sbpy.H == 18.0 + assert pc_sbpy.G == 0.15 + + # create PhaseCurve from sbpy model + pc3 = PhaseCurve().InitModelSbpy(pc_sbpy) + + assert pc3.H == 18.0 + assert pc3.phase_parameter_1 == 0.15 + assert pc3.model_name == "HG" + + +# TODO: test absmag and units +# TODO: test FitModel with weighted data From 602eb980ef4be8fb46f4fd3bb5b381dbc1db413c Mon Sep 17 00:00:00 2001 From: jrob93 Date: Tue, 3 Dec 2024 20:49:37 +0000 Subject: [PATCH 10/10] testing nb --- notebooks/mc_phase_fit_testing.ipynb | 179 ++++----------------------- 1 file changed, 21 insertions(+), 158 deletions(-) diff --git a/notebooks/mc_phase_fit_testing.ipynb b/notebooks/mc_phase_fit_testing.ipynb index 78584e3..cb309e1 100644 --- a/notebooks/mc_phase_fit_testing.ipynb +++ b/notebooks/mc_phase_fit_testing.ipynb @@ -75,7 +75,8 @@ "metadata": {}, "outputs": [], "source": [ - "resample = 100" + "resample = 100\n", + "np.random.seed(0) # set the seed to ensure reproducibility" ] }, { @@ -91,7 +92,7 @@ "\n", "adler_cols = AdlerData(ssoid, planetoid.filter_list)\n", "\n", - "for filt in [\"g\", \"r\"]:\n", + "for filt in [\"r\", \"g\"]:\n", " # get observations and define a phase curve model to fit\n", "\n", " sso = planetoid.SSObject_in_filter(filt)\n", @@ -101,7 +102,7 @@ " G12 = sso.G12\n", " pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name=\"HG12_Pen16\")\n", "\n", - " alpha = (np.array(getattr(obs, \"phaseAngle\")) * u.deg,)\n", + " alpha = np.array(getattr(obs, \"phaseAngle\")) * u.deg\n", " mag_err = np.array(getattr(obs, \"magErr\")) * u.mag\n", " red_mag = np.array(getattr(obs, \"reduced_mag\")) * u.mag\n", "\n", @@ -125,10 +126,11 @@ " )\n", " ax1.legend()\n", "\n", + " # Resample the data within the uncertainties and fit a new model\n", " mc_models = []\n", " for i in range(resample):\n", " _red_mag = np.random.normal(loc=np.array(red_mag), scale=np.array(mag_err)) * u.mag\n", - " ax1.scatter(alpha[0].value, _red_mag.value, s=1, c=\"k\", alpha=0.1, zorder=5)\n", + " ax1.scatter(alpha.value, _red_mag.value, s=1, c=\"k\", alpha=0.1, zorder=5)\n", "\n", " _pc_fit = pc.FitModel(\n", " np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n", @@ -153,117 +155,6 @@ "plt.show()" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "5d0ef23a-ea32-4e32-95ed-a344c4f2760b", - "metadata": {}, - "outputs": [], - "source": [ - "pc_fit.__dict__" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8fa4202e-e7f4-4438-a41e-dc7582312bca", - "metadata": {}, - "outputs": [], - "source": [ - "_pc_fit.param_names" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "18e04044-8007-4b41-85d5-a9b1ad720223", - "metadata": {}, - "outputs": [], - "source": [ - "mc_models[0].__dict__.keys()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f43eb646-cf4e-4103-a9d4-7f8b6c8d2738", - "metadata": {}, - "outputs": [], - "source": [ - "pc_fit_resamp = pc.FitModel(np.array(getattr(obs, \"phaseAngle\")) * u.deg, red_mag, mag_err, resample=resample)\n", - "pc_fit_resamp = PhaseCurve().InitModelSbpy(pc_fit_resamp)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "871dc090-7c0e-41d7-bbc9-9a0ec0232d47", - "metadata": {}, - "outputs": [], - "source": [ - "pc_fit_resamp.__dict__" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3867d1de", - "metadata": {}, - "outputs": [], - "source": [ - "pc_fit_resamp.H, pc_fit_resamp.phase_parameter_1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "71cce253-96f2-44a9-923b-9d4af3c50db4", - "metadata": {}, - "outputs": [], - "source": [ - "pc_fit.__dict__" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6f05c9eb", - "metadata": {}, - "outputs": [], - "source": [ - "pc.__dict__" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "919f0e21", - "metadata": {}, - "outputs": [], - "source": [ - "mag_err" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "321547a5-28ff-416a-b370-93e4f55d4370", - "metadata": {}, - "outputs": [], - "source": [ - "np.median(mag_err)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ea39b9b0", - "metadata": {}, - "outputs": [], - "source": [ - "mc_models[0].H, mc_models[0].phase_parameter_1" - ] - }, { "cell_type": "code", "execution_count": null, @@ -271,6 +162,7 @@ "metadata": {}, "outputs": [], "source": [ + "# extract all phase parameters from the resampling test\n", "abs_mag = np.array([x.H.value for x in mc_models])\n", "phase_param = np.array([x.phase_parameter_1 for x in mc_models])" ] @@ -278,31 +170,14 @@ { "cell_type": "code", "execution_count": null, - "id": "040867f2", - "metadata": {}, - "outputs": [], - "source": [ - "phase_param" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8b21f899", - "metadata": {}, - "outputs": [], - "source": [ - "pc_fit.H, pc_fit.H_err" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1043bc2e", + "id": "f43eb646-cf4e-4103-a9d4-7f8b6c8d2738", "metadata": {}, "outputs": [], "source": [ - "pc.H, pc.H_err" + "# use the resample function within FitModel\n", + "np.random.seed(0) # set the seed to ensure reproducibility\n", + "pc_fit_resamp = pc.FitModel(np.array(getattr(obs, \"phaseAngle\")) * u.deg, red_mag, mag_err, resample=resample)\n", + "pc_fit_resamp = PhaseCurve().InitModelSbpy(pc_fit_resamp)" ] }, { @@ -312,11 +187,14 @@ "metadata": {}, "outputs": [], "source": [ + "# plot the distributions of all FitModel methods\n", + "\n", "fig = plt.figure()\n", "gs = gridspec.GridSpec(1, 2)\n", "ax1 = plt.subplot(gs[0, 0])\n", "ax2 = plt.subplot(gs[0, 1])\n", "\n", + "# plot the results of resampling the data and storing the fits\n", "label = \"resample test\"\n", "for x, ax, ax_lab in zip([abs_mag, phase_param], [ax1, ax2], [\"H\", \"phase_parameter_1\"]):\n", " ax.hist(x, bins=\"auto\", histtype=\"step\")\n", @@ -327,7 +205,9 @@ " ax.axvspan(mean - std, mean + std, alpha=0.3)\n", " ax.set_xlabel(ax_lab)\n", "\n", - "for _pc, col, label in zip([pc, pc_fit_resamp], [\"k\", \"r\"], [\"single fit\", \"resample\"]):\n", + "# first plot results for the regular FitModel result (no resampling)\n", + "# then plot the resulting FitModel with resample results (should be similar to the resample test)\n", + "for _pc, col, label in zip([pc, pc_fit_resamp], [\"k\", \"r\"], [\"FitModel single\", \"FitModel resample\"]):\n", " for x_name, x_err_name, ax in zip(\n", " [\"H\", \"phase_parameter_1\"], [\"H_err\", \"phase_parameter_1_err\"], [ax1, ax2]\n", " ):\n", @@ -353,30 +233,13 @@ { "cell_type": "code", "execution_count": null, - "id": "219b3e89", - "metadata": {}, - "outputs": [], - "source": [ - "pc.__dict__" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "aa2bb33a", + "id": "c6a0cbce", "metadata": {}, "outputs": [], "source": [ - "_pc.__dict__" + "# resampling leads to a slight shift in fitted values\n", + "# the resampled distribution has a slightly larger \"uncertainty\"" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c6a0cbce", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": {