Skip to content

Commit

Permalink
unit tests and demo nb
Browse files Browse the repository at this point in the history
  • Loading branch information
@jrob93
jrob93 committed May 13, 2024
1 parent 1eeed56 commit 83d917e
Show file tree
Hide file tree
Showing 2 changed files with 233 additions and 178 deletions.
213 changes: 35 additions & 178 deletions notebooks/outlier_detection.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,16 @@
"res = obs_r.reduced_mag - pc.ReducedMag(obs_r.phaseAngle * u.degree).value"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "81e04bad",
"metadata": {},
"outputs": [],
"source": [
"res"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -77,13 +87,13 @@
"\n",
"x = getattr(obs_r, x_plot)\n",
"y = getattr(obs_r, y_plot)\n",
"xerr = obs_r.magErr\n",
"yerr = obs_r.magErr\n",
"\n",
"fig = plt.figure()\n",
"gs = gridspec.GridSpec(1, 1)\n",
"ax1 = plt.subplot(gs[0, 0])\n",
"\n",
"ax1.errorbar(x, y, xerr, fmt=\"o\")\n",
"ax1.errorbar(x, y, yerr, fmt=\"o\")\n",
"\n",
"# plot the phase curve model\n",
"alpha = np.linspace(0, np.amax(obs_r.phaseAngle)) * u.deg\n",
Expand Down Expand Up @@ -116,13 +126,13 @@
"x_plot = \"phaseAngle\"\n",
"\n",
"x = getattr(obs_r, x_plot)\n",
"xerr = obs_r.magErr\n",
"yerr = obs_r.magErr\n",
"\n",
"fig = plt.figure()\n",
"gs = gridspec.GridSpec(1, 1)\n",
"ax1 = plt.subplot(gs[0, 0])\n",
"\n",
"ax1.errorbar(x, res, xerr, fmt=\"o\")\n",
"ax1.errorbar(x, res, yerr, fmt=\"o\")\n",
"\n",
"ax1.axhline(0, c=\"k\")\n",
"\n",
Expand All @@ -143,186 +153,55 @@
{
"cell_type": "code",
"execution_count": null,
"id": "16c712ff",
"metadata": {},
"outputs": [],
"source": [
"utils.sigma_clip(res, {\"cenfunc\": \"median\"})"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8983c9dd",
"metadata": {},
"outputs": [],
"source": [
"utils.sigma_clip(res)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15bdb74e",
"metadata": {},
"outputs": [],
"source": [
"utils.outlier_std(res, res)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4369f850",
"metadata": {},
"outputs": [],
"source": [
"df_obs_r = pd.DataFrame(obs_r.__dict__)\n",
"df_obs_r = df_obs_r.sort_values(\"midpointMjdTai\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "890a9eff",
"metadata": {},
"outputs": [],
"source": [
"df_obs_r"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c24dd432",
"id": "13e25b85",
"metadata": {},
"outputs": [],
"source": [
"# step through all observations and check if photometry is outlying\n",
"for i in range(1, len(df_obs_r) + 1):\n",
" # print(df_obs_r.iloc[:i].index)\n",
" _df = df_obs_r.iloc[:i]\n",
"\n",
" _res = _df[\"reduced_mag\"] - pc.ReducedMag(np.array(_df[\"phaseAngle\"]) * u.degree).value\n",
" # print(_res)\n",
"\n",
" # print(utils.outlier_std(_res,res))\n",
" print(utils.outlier_std(_res, _res))"
"# return a list of flags for outlying objects"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2d21af47",
"id": "16c712ff",
"metadata": {},
"outputs": [],
"source": [
"# step through all observations and check if new photometry is outlying\n",
"\n",
"fig = plt.figure()\n",
"gs = gridspec.GridSpec(1, 1)\n",
"ax1 = plt.subplot(gs[0, 0])\n",
"\n",
"for i in range(2, len(df_obs_r) + 1):\n",
" # print(df_obs_r.iloc[:i].index)\n",
" _df = df_obs_r.iloc[:i]\n",
"\n",
" _res = _df[\"reduced_mag\"] - pc.ReducedMag(np.array(_df[\"phaseAngle\"]) * u.degree).value\n",
" _res = np.array(_res)\n",
"\n",
" # print(_res)\n",
" # print(np.delete(_res,-1),[_res[-1]])\n",
"\n",
" outlier_flag = utils.outlier_std([_res[-1]], np.delete(_res, -1))\n",
" print(outlier_flag, np.abs(_res[-1]), 3 * np.std(np.delete(_res, -1)))\n",
"\n",
" s1 = ax1.scatter(\n",
" _df.iloc[-1][\"phaseAngle\"],\n",
" _res[-1],\n",
" c=_df.iloc[-1][\"midpointMjdTai\"],\n",
" vmin=np.amin(df_obs_r[\"midpointMjdTai\"]),\n",
" vmax=np.amax(df_obs_r[\"midpointMjdTai\"]),\n",
" )\n",
"\n",
" if outlier_flag[0]:\n",
" ax1.scatter(_df.iloc[-1][\"phaseAngle\"], _res[-1], c=\"r\", marker=\"x\")\n",
"\n",
"ax1.axhline(0, c=\"k\")\n",
"plt.colorbar(s1)\n",
"ax1.invert_yaxis()\n",
"\n",
"plt.show()"
"# astropy sigma_clip normally uses the median as the central function\n",
"utils.sigma_clip(res, {\"cenfunc\": \"median\"})"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4231b25",
"id": "8983c9dd",
"metadata": {},
"outputs": [],
"source": [
"# step through all observations and check if new photometry is outlying\n",
"\n",
"fig = plt.figure()\n",
"gs = gridspec.GridSpec(1, 1)\n",
"ax1 = plt.subplot(gs[0, 0])\n",
"\n",
"std_list = []\n",
"for i in range(2, len(df_obs_r) + 1):\n",
" # print(df_obs_r.iloc[:i].index)\n",
" _df = df_obs_r.iloc[:i]\n",
"\n",
" _res = _df[\"reduced_mag\"] - pc.ReducedMag(np.array(_df[\"phaseAngle\"]) * u.degree).value\n",
" _res = np.array(_res)\n",
"\n",
" # print(_res)\n",
" # print(np.delete(_res,-1),[_res[-1]])\n",
"\n",
" outlier_flag = utils.outlier_std([_res[-1]], np.delete(_res, -1))\n",
" print(outlier_flag, np.abs(_res[-1]), 3 * np.std(np.delete(_res, -1)))\n",
"\n",
" s1 = ax1.scatter(\n",
" _df.iloc[-1][\"midpointMjdTai\"],\n",
" _res[-1],\n",
" c=_df.iloc[-1][\"phaseAngle\"],\n",
" vmin=np.amin(df_obs_r[\"phaseAngle\"]),\n",
" vmax=np.amax(df_obs_r[\"phaseAngle\"]),\n",
" )\n",
"\n",
" if outlier_flag[0]:\n",
" ax1.scatter(_df.iloc[-1][\"midpointMjdTai\"], _res[-1], c=\"r\", marker=\"x\")\n",
"\n",
" std_list.append(np.std(np.delete(_res, -1)))\n",
"\n",
"std_list = np.array(std_list)\n",
"\n",
"ax1.plot(np.array(df_obs_r[\"midpointMjdTai\"])[1:], std_list * 3)\n",
"\n",
"ax1.axhline(0, c=\"k\")\n",
"plt.colorbar(s1)\n",
"ax1.invert_yaxis()\n",
"\n",
"plt.show()"
"# assuming that the model is the ground truth, we use zero as the centroid for the residuals\n",
"utils.sigma_clip(res, {\"cenfunc\": utils.zero_func})"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8910801b",
"id": "15bdb74e",
"metadata": {},
"outputs": [],
"source": [
"std_list"
"# use the standard deviation of the residuals to identify outliers\n",
"utils.outlier_std(res, res, std_cut=3.0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "46e31e45",
"id": "50276ce5",
"metadata": {},
"outputs": [],
"source": [
"df_obs_r[\"midpointMjdTai\"]"
"# use a simple threshold value for residuals to find outliers\n",
"utils.outlier_diff(res, diff_cut=1.0)"
]
},
{
Expand All @@ -332,7 +211,9 @@
"metadata": {},
"outputs": [],
"source": [
"# consider the residual compared to the uncertainty"
"# consider the residual compared to the uncertainty of the measurement\n",
"std_err = 5\n",
"utils.outlier_sigma_diff(res, yerr, std_err)"
]
},
{
Expand All @@ -346,13 +227,13 @@
"x_plot = \"phaseAngle\"\n",
"\n",
"x = getattr(obs_r, x_plot)\n",
"xerr = obs_r.magErr\n",
"yerr = obs_r.magErr\n",
"\n",
"fig = plt.figure()\n",
"gs = gridspec.GridSpec(1, 1)\n",
"ax1 = plt.subplot(gs[0, 0])\n",
"\n",
"ax1.errorbar(x, res, xerr, fmt=\"o\")\n",
"ax1.errorbar(x, res, yerr, fmt=\"o\")\n",
"\n",
"ax1.axhline(0, c=\"k\")\n",
"\n",
Expand All @@ -362,11 +243,7 @@
" ax1.axhline(res_std * i, ls=\":\", c=\"C{}\".format(i), label=\"{} sigma\".format(i))\n",
" ax1.axhline(-res_std * i, ls=\":\", c=\"C{}\".format(i))\n",
"\n",
"std_err = 5\n",
"mask = np.abs(res) > (std_err * xerr)\n",
"ax1.scatter(x[mask], res[mask], c=\"r\", marker=\"x\", zorder=3)\n",
"\n",
"mask2 = utils.outlier_sigma_diff(res, xerr, std_err)\n",
"mask2 = utils.outlier_sigma_diff(res, yerr, std_err)\n",
"ax1.scatter(x[mask], res[mask], edgecolor=\"r\", facecolor=\"none\", marker=\"o\", zorder=3)\n",
"\n",
"ax1.invert_yaxis()\n",
Expand All @@ -384,27 +261,7 @@
"metadata": {},
"outputs": [],
"source": [
"# in general the residuals can be much larger than the photometric uncertainty"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad0e0c40",
"metadata": {},
"outputs": [],
"source": [
"res"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3ce23237",
"metadata": {},
"outputs": [],
"source": [
"xerr"
"# NB that for phase curve models, residuals can be much larger than the photometric uncertainty!"
]
},
{
Expand Down
Loading

0 comments on commit 83d917e

Please sign in to comment.