import sys
sys.path.insert(0, '../')

import subprocess
import pickle
from agasc_gaia import cross_match as xm
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from astropy.table import join
from pathlib import Path

from agasc_gaia import cross_match as xm, gaia_model as gm, config, datasets
from agasc_gaia import star_report

# to silence warnings since 2024.1
import warnings
warnings.simplefilter("ignore")

REPORTS_DIR = Path("/Users/javierg/SAO/Notebooks/agasc/gaia-magnitudes-2023/reports")

agasc_difficult = xm.get_agasc_gaia_x_match_difficult()
agasc_gaia_matches = xm.get_agasc_gaia_x_match()
# agasc_gaia_matches_all = xm.get_agasc_gaia_x_match_all()
agasc_summary = datasets.get_agasc_summary()

agasc_difficult = agasc_difficult.group_by(['group'])
arr = np.array(agasc_difficult['agasc_id'])
agasc_difficult['group_size'] = np.concatenate([
    np.repeat(np.unique(arr[i:j]).shape[0], j - i)
    for i, j in zip(agasc_difficult.groups.indices[:-1], agasc_difficult.groups.indices[1:])
])

# some difficult stars have a best match with NaN p_value
# this is NOT a problem, because agasc_difficult includes all _possible_ matches, not just the best ones
# These will get removed downstream
(
    np.count_nonzero(np.isnan(agasc_gaia_matches['mag_pred'])),  # should be zero
    np.count_nonzero(np.isnan(agasc_difficult['p_value']) & agasc_difficult['best_match'])
)

(0, 29)

# agasc_difficult = agasc_difficult[agasc_difficult['best_match']]

i = np.searchsorted(agasc_summary["agasc_id"], agasc_gaia_matches["agasc_id"])
cols = [
    "mag_aca",
    "mag_aca_err",
    "mag_aca_obs",
    "mag_aca_err_obs",
    "mag_catid",
    "mag_band",
    "random_index",
    "guide",
    "acq"
]
for col in cols:
    agasc_gaia_matches[col] = agasc_summary[col][i]

cols = [
    "pm_ra",
    "pm_dec",
    "epoch",
    "ra",
    "dec",
]
for col in cols:
    agasc_gaia_matches[f"{col}_agasc"] = agasc_summary[col][i]

new_cols = [f"{col}_gaia" for col in cols]
agasc_gaia_matches.rename_columns(cols, new_cols)

# fixing some numbers in AGASC (assuming that entries with -9999 have no proper motion)
agasc_gaia_matches['pm_dec_agasc'][agasc_gaia_matches['pm_dec_agasc'] == -9999] = 0
agasc_gaia_matches['pm_ra_agasc'][agasc_gaia_matches['pm_ra_agasc'] == -9999] = 0

with open(config.FILES["gaia_model"], "rb") as fh:
    gaia_model = pickle.load(fh)

agasc_gaia_matches["mag_aca_pred"] = gaia_model.predict(
    agasc_gaia_matches.to_pandas(),
    with_instrument_bias=False
)

agasc_gaia_matches["mag_aca_pred_err"] = gaia_model.uncertainty(
    agasc_gaia_matches.to_pandas(),
    with_instrument=False
)

Proper Motion

Proper motion is a nice way to see whether the cross-match is generally fine, because proper motion is not used directly in the cross-match. It is used to shift positions to a common epoch, but it is not considered in the matching probability.

The vertical line is all the stars with no PM in agasc and with PM in Gaia

fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)
plt.sca(axes[0])

plt.title("PM RA")
x=np.array(agasc_gaia_matches['pm_ra_agasc'], dtype=float)
y=np.array(agasc_gaia_matches['pm_ra_gaia'], dtype=float)
bins = np.linspace(-200, 200, 200)
sns.histplot(
    x=x, y=y,
    bins=[bins, bins],
)
plt.xlabel("AGASC 1p7")
plt.ylabel("AGASC 1p8")

plt.sca(axes[1])
plt.title("PM DEC")
x=np.array(agasc_gaia_matches['pm_dec_agasc'], dtype=float)
y=np.array(agasc_gaia_matches['pm_dec_gaia'], dtype=float)
bins = np.linspace(-200, 200, 200)
sns.histplot(
    x=x, y=y,
    bins=[bins, bins],
)
plt.xlabel("AGASC 1p7")
plt.ylabel("AGASC 1p8")

plt.tight_layout()

fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)
plt.sca(axes[0])
sel = np.abs(agasc_gaia_matches['pm_dec_agasc']) == 0
plt.hist(agasc_gaia_matches['pm_dec_gaia'][sel], bins=np.linspace(-800, 800, 100));
plt.sca(axes[1])
sel = np.abs(agasc_gaia_matches['pm_ra_agasc']) == 0
plt.hist(agasc_gaia_matches['pm_ra_gaia'][sel], bins=np.linspace(-800, 800, 100));

plt.yscale('log')
plt.suptitle('Proper motion Updates from Gaia')

Text(0.5, 0.98, 'Proper motion Updates from Gaia')

fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)
plt.sca(axes[0])
plt.title("PM RA")
sel = np.abs(agasc_gaia_matches['pm_ra_agasc']) > 0
x=np.array(agasc_gaia_matches['pm_ra_agasc'][sel], dtype=float)
y=np.array(agasc_gaia_matches['pm_ra_gaia'][sel], dtype=float)
bins = np.linspace(-200, 200, 200)
sns.histplot(
    x=x, y=y,
    bins=[bins, bins],
)
plt.axhline(10, color='k', linestyle='--')
plt.axhline(-10, color='k', linestyle='--')
x = np.linspace(-200, 200, 100)
plt.plot(x, x+20, color='k', linestyle='--')
plt.plot(x, x-20, color='k', linestyle='--')
plt.xlabel("AGASC 1p7")
plt.ylabel("AGASC 1p8")

plt.sca(axes[1])
plt.title("PM DEC")
sel = np.abs(agasc_gaia_matches['pm_dec_agasc']) > 0
x=np.array(agasc_gaia_matches['pm_dec_agasc'][sel], dtype=float)
y=np.array(agasc_gaia_matches['pm_dec_gaia'][sel], dtype=float)
bins = np.linspace(-200, 200, 200)
sns.histplot(
    x=x, y=y,
    bins=[bins, bins],
)
plt.axhline(10, color='k', linestyle='--')
plt.axhline(-10, color='k', linestyle='--')
x = np.linspace(-200, 200, 100)
plt.plot(x, x+20, color='k', linestyle='--')
plt.plot(x, x-20, color='k', linestyle='--')
plt.xlabel("AGASC 1p7")
plt.ylabel("AGASC 1p8")
plt.tight_layout()

sel = (
    (
        (np.abs(agasc_gaia_matches['pm_dec_agasc']) > 0)
        & (np.abs(agasc_gaia_matches['pm_dec_gaia']) < 10)
        & (np.abs(agasc_gaia_matches['pm_dec_gaia'] - agasc_gaia_matches['pm_dec_agasc']) > 20)
    )
    | (
        (np.abs(agasc_gaia_matches['pm_ra_agasc']) > 0)
        & (np.abs(agasc_gaia_matches['pm_ra_gaia']) < 10)
        & (np.abs(agasc_gaia_matches['pm_ra_gaia'] - agasc_gaia_matches['pm_ra_agasc']) > 20)
    )
)
pm_outliers = agasc_gaia_matches[sel]
pm_outliers['d_pm'] = np.sqrt(
    (pm_outliers['pm_ra_gaia'] - pm_outliers['pm_ra_agasc'])**2
    + (pm_outliers['pm_dec_gaia'] - pm_outliers['pm_dec_agasc'])**2
)
pm_outliers['d_pm'].format = "{:.1f}"
pm_outliers['p_value'].format = "{:.3f}"
pm_outliers['d2d'].format = "{:.2f}"
pm_outliers['d_mag'].format = "{:.2f}"
pm_outliers.sort('d_pm', reverse=True)
len(pm_outliers)

3265

fig, axes = plt.subplot_mosaic(
    [
        ['A', '.'],
        ['B', 'C']
    ],
    # sharex=True, sharey=True,
    figsize=(8, 6),
)
plt.sca(axes['A'])
sns.histplot(x=pm_outliers['d_pm'])
plt.xlabel('')
plt.sca(axes['B'])
sns.histplot(x=pm_outliers['d_pm'], y=pm_outliers['p_value'])
plt.xlabel(r'$\Delta PM$')
plt.sca(axes['C'])
sns.histplot(y=pm_outliers['p_value'])
plt.ylabel('')
plt.suptitle('p-value $\Delta PM$ distribution')
plt.tight_layout()

fig, axes = plt.subplot_mosaic(
    [
        ['A', '.'],
        ['B', 'C']
    ],
    # sharex=True, sharey=True,
    figsize=(8, 6),
)
plt.sca(axes['A'])
sns.histplot(x=pm_outliers['d2d'])
plt.xlabel('')
plt.yscale('log')
plt.sca(axes['B'])
sns.histplot(x=pm_outliers['d2d'], y=pm_outliers['d_mag'])
plt.xlabel('d2d')
plt.sca(axes['C'])
sns.histplot(y=pm_outliers['d_mag'])
plt.ylabel('')
plt.xscale('log')

plt.suptitle('PM outliers d2d-$\Delta mag$ distribution')
plt.tight_layout()

agasc_id = pm_outliers['agasc_id'][pm_outliers['d2d'] == np.max(pm_outliers['d2d'])][0]
report = star_report.Report(agasc_id)
report.show_in_notebook()

AGASC 930611512 epoch in AGASC: 2000.0, epoch in summary: 1991.59

AGASC entry

Table length=1

agasc_id	ra	dec	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	epoch	tycho_id	gsc2.3
int32	float64	float64	float32	float32	float32	int16	int16	float32	str12	str10
930611512	101.28444394	-37.39901464	12.34	11.95	--	31	-9999	2000.0	7100-39-1	S3DP000378

Matches

Table length=1

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
930611512	5577278955302942848	12.34	12.32	-0.02	1.90	0.0012	0.0061	1.000	-1.912	1.835	11.97	True

Match Candidates

Table length=2

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
930611512	5577278955302942848	12.34	12.32	-0.02	1.90	0.0012	0.0061	1.000	-1.912	1.835	11.97	True
930611512	5577278951007783424	12.34	19.05	6.70	13.85	0.0000	0.0000	0.000	-1.382	3.807	18.98	False

AGASC Neighborhood (within 20 arcsec)

Table length=1

agasc_id	ra	dec	d2d	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	pos_catid	pm_catid	mag_catid
			arcsec
int32	float64	float64	float64	float32	float32	float32	int16	int16	uint8	uint8	uint8
930611512	101.28444394	-37.39901464	0.05	12.34	11.95	--	31	-9999	5	5	5

# the largest pm outlier
i = 0
agasc_id = pm_outliers['agasc_id'][i]
report = star_report.Report(agasc_id)
report.show_in_notebook()

AGASC 349707128 epoch in AGASC: 2000.0, epoch in summary: 1991.54

AGASC entry

Table length=1

agasc_id	ra	dec	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	epoch	tycho_id	gsc2.3
int32	float64	float64	float32	float32	float32	int16	int16	float32	str12	str10
349707128	294.96609948	36.01583202	11.56	11.49	--	278	-403	2000.0	2668-879-1	N2GJ000523

Matches

Table length=1

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
349707128	2048222903655254912	11.56	11.44	-0.12	0.16	1.9190	0.6987	1.000	2.449	-0.270	11.31	True

Match Candidates

Table length=12

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
349707128	2048222903655254912	11.56	11.44	-0.12	0.16	1.9190	0.6987	1.000	2.449	-0.270	11.31	True
349707128	2048222903658885376	11.56	17.61	6.05	3.24	0.0000	0.0000	0.000	-3.582	-6.793	19.72	False
349707128	2048222697485766784	11.56	20.95	9.39	17.23	0.0000	0.0000	0.000	--	--	20.95	False
349707128	2048222903655622272	11.56	19.80	8.24	6.12	0.0000	0.0000	0.000	-2.441	-4.730	20.28	False
349707128	2048222903658884992	11.56	19.06	7.50	6.43	0.0000	0.0000	0.000	-5.891	-6.789	19.88	False
349707128	2048222903643368832	11.56	19.25	7.70	8.31	0.0000	0.0000	0.000	-1.080	-0.110	20.48	False
349707128	2048222903655259520	11.56	14.94	3.38	9.77	0.0000	0.0000	0.000	-3.904	-0.028	14.49	False
349707128	2048222903658885248	11.56	13.73	2.18	10.52	0.0000	0.0000	0.000	-4.348	-6.238	13.71	False
349707128	2048222903655260800	11.56	15.37	3.81	10.90	0.0000	0.0000	0.000	-0.690	-3.924	15.375	False
349707128	2048222903642138496	11.56	19.05	7.50	10.64	0.0000	0.0000	0.000	-3.686	-4.805	19.3	False
349707128	2048222903642130432	11.56	19.28	7.73	11.95	0.0000	0.0000	0.000	-3.701	-5.312	19.3	False
349707128	2048222903642149760	11.56	19.70	8.15	13.66	0.0000	0.0000	0.000	-1.597	-5.020	19.75	False

AGASC Neighborhood (within 20 arcsec)

Table length=1

agasc_id	ra	dec	d2d	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	pos_catid	pm_catid	mag_catid
			arcsec
int32	float64	float64	float64	float32	float32	float32	int16	int16	uint8	uint8	uint8
349707128	294.96609948	36.01583202	0.03	11.56	11.49	--	278	-403	5	5	5

generate_reports = True
np.random.seed(349707128)
if generate_reports:
    N = 100
    idx = np.arange(len(pm_outliers))
    sel = np.zeros(len(pm_outliers), dtype=bool)
    sel[:N] = True
    print(f"Making report with {np.count_nonzero(sel)} stars")

    description = """
    <p>This list includes the 100 PM outliers with the largest proper motion differences </p>

    <p> Note that the stars with the largest proper motion difference are not the ones with the worst p-value. </p>

    <p>
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>
    """

    reports_dir = REPORTS_DIR / "pm-outliers/large-pm-diff"
    star_report.make_report_list(
        data=pm_outliers[['agasc_id', 'p_value', 'd2d', 'd_mag', 'd_pm']][sel],
        path=reports_dir,
        title='Proper Motion Outliers',
        description=description,
        overwrite=True
    )

Making report with 100 stars

100it [04:04,  2.45s/it]

if generate_reports:
    N = 100
    idx = np.arange(len(pm_outliers))
    sel = np.zeros(len(pm_outliers), dtype=bool)
    sel[np.argsort(pm_outliers['p_value'])[:N]] = True
    print(f"Making report with {np.count_nonzero(sel)} stars")

    description = """
    <p>This list includes the 100 PM outliers with the worst p-value </p>

    <p> Note that the stars with the largest proper motion difference are not the ones with the worst p-value. </p>

    <p>
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>
    """

    reports_dir = REPORTS_DIR / "pm-outliers" / "worst-p-value"
    star_report.make_report_list(
        data=pm_outliers[['agasc_id', 'p_value', 'd2d', 'd_mag', 'd_pm']][sel],
        path=reports_dir,
        title='Proper Motion Outliers',
        description=description,
        overwrite=True
    )

Making report with 100 stars

100it [03:58,  2.38s/it]

if generate_reports:
    idx = np.arange(len(pm_outliers))
    sel = np.zeros(len(pm_outliers), dtype=bool)
    sel[np.random.choice(idx[~sel], replace=False, size=len(pm_outliers)//40)] = True
    print(f"Making report with {np.count_nonzero(sel)} stars")

    description = """
    <p>This list includes 1/40th of the proper motion outliers randomly selected. </p>

    <p> Note that the stars with the largest proper motion difference are not the ones with the worst p-value. </p>

    <p>
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>
    """

    reports_dir = REPORTS_DIR / "pm-outliers" / "random"
    star_report.make_report_list(
        data=pm_outliers[['agasc_id', 'p_value', 'd2d', 'd_mag', 'd_pm']][sel],
        path=reports_dir,
        title='Proper Motion Outliers',
        description=description,
        overwrite=True
    )

Making report with 81 stars

81it [02:59,  2.22s/it]

Catalog Magnitude Outliers

These are stars with large differences between AGASC 1.7 and AGASC 1.8 magnitude.

# the astype is needed or I get an endianness error
agasc_gaia_matches["mag_aca"] = agasc_gaia_matches["mag_aca"].astype(np.float64)

sns.histplot(
    x=agasc_gaia_matches["mag_aca"],
    y=agasc_gaia_matches["mag_aca_pred"] - agasc_gaia_matches["mag_aca"],
    bins=(np.linspace(5, 15, 101), np.linspace(-2, 2, 101)),
    cbar_kws=dict(label='Number of stars'),
)

<Axes: xlabel='mag_aca', ylabel='mag_aca_pred'>

fig, axes = plt.subplot_mosaic(
    [
        ['.', 'all', 'all', '.'],
        ['bright', 'bright', 'faint', 'faint'],
    ],
    # sharex=True, sharey=True,
    figsize=(12, 8),
)
plt.sca(axes['all'])
plt.hist(
    agasc_gaia_matches["mag_aca_pred"] - agasc_gaia_matches["mag_aca"],
    bins=np.linspace(-6, 6, 100),
    histtype="step",
)
plt.yscale('log')
plt.xlabel('Predicted - Catalog ACA Magnitude')
plt.title("all")

plt.sca(axes['bright'])
sel = (agasc_gaia_matches["mag_aca"] < 8.5)
plt.hist(
    agasc_gaia_matches["mag_aca_pred"][sel] - agasc_gaia_matches["mag_aca"][sel],
    bins=np.linspace(-6, 6, 100),
    histtype="step",
)
plt.yscale('log')
plt.xlabel('Predicted - Catalog ACA Magnitude')
plt.title('mag_aca < 8.5')

plt.sca(axes['faint'])
sel = (agasc_gaia_matches["mag_aca"] > 9.5) & (agasc_gaia_matches["mag_aca"] < 10.3)
plt.hist(
    agasc_gaia_matches["mag_aca_pred"][sel] - agasc_gaia_matches["mag_aca"][sel],
    bins=np.linspace(-6, 6, 100),
    histtype="step",
)
plt.yscale('log')
plt.title('mag_aca > 9.5')
plt.xlabel('Predicted - Catalog ACA Magnitude')

plt.suptitle('Predicted - Catalog ACA Magnitude')

plt.tight_layout()

sel = (
    (agasc_gaia_matches["mag_aca"] < 10.3)
    & (np.abs(agasc_gaia_matches["mag_aca_pred"] - agasc_gaia_matches["mag_aca"]) > 2.8)
)
mag_outliers = agasc_gaia_matches[sel]

mag_outliers['d_mag_aca'] = mag_outliers["mag_aca_pred"] - mag_outliers["mag_aca"]
mag_outliers['d_mag_aca'].format = "{:.2f}"
mag_outliers['d_mag'].format = "{:.2f}"
mag_outliers['mag_aca'].format = "{:.2f}"
mag_outliers['mag_aca_pred'].format = "{:.2f}"
mag_outliers['mag_aca_obs'].format = "{:.2f}"
mag_outliers['d2d'].format = "{:.3f}"
mag_outliers['p_value'].format = "{:.4f}"

mag_outliers['abs_d_mag_aca'] = np.abs(mag_outliers['d_mag_aca'])


mag_outliers.sort('mag_aca')
# mag_outliers.sort('abs_d_mag_aca', reverse=True)

np.random.seed(1174545944)
if generate_reports:
    print(f"Making report with {len(mag_outliers)} stars")

    description = """
    <p>This list includes the stars with abs(mag_aca_1p8 - mag_aca_1p7) > 2.8 </p>

    <p>
    Note that there are several stars with observed mag_aca. In all cases, the observed mag_aca
    is closer to the 1p8 one.
    </p>

    <p>
    The purpose of this report is to check whether outliers can be caused by misidentification.
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>

    <h2> Notable examples:</h2>
    <ul>
      <li>
        <a href="report_879886528.html"> 879886528 </a>. An example of duplicate
        record in AGASC that just happens to be in this dataset.
      </li>
      <li>
        <a href="report_796801904.html"> 796801904 </a>. An example of duplicate
        record in AGASC that just happens to be in this dataset.
      </li>
      <li>
        <a href="report_412493456.html"> 412493456 </a>. A tycho2 pair of stars
        with the same position in Tycho2, but resolved in Gaia. This is the brighter of the two.
      </li>
      <li>
        <a href="report_298718560.html"> 298718560 </a>. One star in AGASC that is
        actually two stars 2 arcsec apart in Gaia. The AGASC star is cross-matched with the brighter
        one, but strictly speaking it should be a blend of the two or a new record is needed.
      </li>
      <li>
        <a href="report_549978408.html"> 549978408 </a>. An example of duplicate
        record in AGASC that just happens to be in this dataset.
      </li>

    </ul>

    <h2> Stars </h2>
    """

    reports_dir = REPORTS_DIR / "mag_aca-outliers"
    star_report.make_report_list(
        data=mag_outliers[['agasc_id', 'p_value', 'd2d', 'd_mag', 'd_mag_aca', 'mag_aca', 'mag_aca_pred', 'mag_aca_obs']],
        path=reports_dir,
        title='Catalog Magnitude Outliers',
        description=description,
        overwrite=True
    )

Making report with 125 stars

125it [05:06,  2.45s/it]

agasc_gaia_matches[(agasc_gaia_matches["mag_aca_pred"] - agasc_gaia_matches["mag_aca"]) > 5]

Table length=1611

agasc_id	ra_agasc	dec_agasc	ra_2016	dec_2016	pm_ra_agasc	pm_dec_agasc	epoch_agasc	mag	mag_catid	mag_band	tycho_id	gaia_id	ra_gaia	dec_gaia	ra_error	dec_error	pm_ra_gaia	pm_dec_gaia	epoch_gaia	g_mag	bp_mag	rp_mag	g_mag_error	bp_mag_error	rp_mag_error	range_mag_g_fov	std_dev_mag_g_fov	mad_mag_g_fov	phot_variable_flag	phot_proc_mode	non_single_star	has_mag	ave_gaia_mag	d2d	d2d_2016	mag_pred	d_mag	p_match_gaia_pm	p_match_agasc_pm	use_agasc_pm	p_match	log_p_match	best_match	p_value	log_p_value	p_relative	log_p_relative	idx	mag_aca	mag_aca_err	mag_aca_obs	mag_aca_err_obs	random_index	guide	acq	mag_aca_pred	mag_aca_pred_err
int64	float64	float64	float64	float64	int16	int16	float32	float32	uint8	int16	bytes12	int64	float64	float64	float64	float64	float16	float16	float16	float16	float16	float16	float16	float16	float16	float32	float32	float16	bool	int16	int16	int16	float16	float16	float16	float16	float16	float32	float32	bool	float32	float16	bool	float32	float16	float32	float16	int64	float64	int16	float32	float32	int64	bool	bool	float64	float64
133304	0.89048	2.13766	0.89048	2.13766	0	0	1982.0	15.47	1	0	--	2739048842507452800	0.8904937433606026	2.13763802887449	1.033177375793457	0.9319427013397217	--	--	2016.0	19.39	17.55	16.05	0.007935	0.009346	0.004116	--	--	--	False	0	0	7	17.66	0.09326	0.09326	18.05	2.572	0.0025324202	0.0025324202	False	0.0025324202	2.596	True	0.009785407	2.01	1.0	-0.0	383	14.316000938415527	54	--	--	8923154	False	False	19.436804224803925	0.025585260300380372
269488	4.58166	1.76817	4.58166	1.76817	0	0	1982.0	15.01	1	0	--	2547955370867550592	4.581674296276818	1.7681665860575724	2.0221760272979736	1.3425837755203247	--	--	2016.0	19.81	17.56	16.08	0.00927	0.015434	0.003736	--	--	--	False	0	0	7	17.81	0.0529	0.0529	18.1	3.086	0.0015535132	0.0015535132	False	0.0015535132	2.809	True	0.007106245	2.148	1.0	-0.0	2479	13.856000900268555	52	--	--	10772157	False	False	19.861705826058106	0.025585260300380372
270128	4.14213	2.32568	4.14213	2.32568	0	0	1983.0	14.45	1	1	--	2548239182306286848	4.142067143773359	2.325660823209443	1.2723010778427124	0.8178788423538208	--	--	2016.0	19.44	17.03	15.63	0.007652	0.01584	0.003914	--	--	--	False	0	0	7	17.38	0.2365	0.2365	17.33	2.885	0.0012634854	0.0012634854	False	0.0012634854	2.898	True	0.0062418072	2.205	1.0	-0.0	2561	13.904999732971191	50	--	--	9624791	False	False	19.503698638399417	0.025585260300380372
398304	5.31811	1.97925	5.31811	1.97925	0	0	1983.0	14.42	1	1	--	2547311709888343680	5.318126887248783	1.9792420490650975	0.8330377340316772	0.5575965642929077	--	--	2016.0	19.31	17.0	15.6	0.006844	0.014946	0.003115	--	--	--	False	0	0	7	17.3	0.06714	0.06714	17.27	2.854	0.0019283851	0.0019283851	False	0.0019283851	2.715	True	0.008114262	2.09	1.0	-0.0	3620	13.875	51	--	--	8220216	False	False	19.378698638399417	0.025585260300380372
526248	2.21825	3.19358	2.21825	3.19358	0	0	1983.0	14.42	1	1	--	2740732023009424512	2.2181985152813493	3.193578059567273	1.4910318851470947	1.236511468887329	--	--	2016.0	19.6	17.14	15.766	0.008804	0.01637	0.00336	--	--	--	False	0	0	7	17.5	0.1852	0.1852	17.45	3.035	0.0012827248	0.0012827248	False	0.0012827248	2.893	True	0.0062553785	2.203	1.0	-0.0	4457	13.875	53	--	--	1939182	False	False	19.66468480275236	0.025585260300380372
786480	7.48246	3.35412	7.48246	3.35412	0	0	1983.0	14.38	1	1	--	2553848100356689408	7.482475728296715	3.3541514356613793	1.9535853862762451	1.158104658126831	--	--	2016.0	19.52	17.17	15.69	0.01209	0.01464	0.00635	--	--	--	False	0	0	7	17.45	0.1265	0.1265	17.44	3.064	0.0014324953	0.0014324953	False	0.0014324953	2.844	True	0.006720357	2.172	1.0	-0.0	7668	13.835000038146973	53	--	--	2780191	False	False	19.564830826058106	0.025585260300380372
794496	6.32801	2.60568	6.32801	2.60568	0	0	1982.0	14.8	1	0	--	2547742684087792128	6.3280271145577025	2.6057047286506854	1.5639628171920776	0.805314302444458	--	--	2016.0	19.36	17.27	15.805	0.007305	0.01148	0.003582	--	--	--	False	0	0	7	17.47	0.1082	0.1082	17.8	3.002	0.0015730215	0.0015730215	False	0.0015730215	2.803	True	0.007139939	2.146	1.0	-0.0	8811	13.646000862121582	54	--	--	9012161	False	False	19.413157085819982	0.025585260300380372
796216	5.15835	4.55606	5.15835	4.55606	0	0	1983.0	14.99	1	1	--	2741160214071044224	5.158358291993044	4.556056282729221	3.6331419944763184	1.7934240102767944	--	--	2016.0	19.86	17.33	15.875	0.0079	0.01587	0.00418	--	--	--	False	0	0	7	17.69	0.03262	0.03262	17.66	2.664	0.0024221775	0.0024221775	False	0.0024221775	2.615	True	0.009572603	2.02	1.0	-0.0	9078	14.444999694824219	46	--	--	13445884	False	False	19.914692192164942	0.025585260300380372
796544	6.53406	3.9248	6.53406	3.9248	0	0	1983.0	15.2	1	1	--	2548736299001732864	6.534034909054615	3.9248407244736665	6.735021591186523	3.3833651542663574	--	--	2016.0	20.17	17.55	16.12	0.01212	0.01816	0.005554	--	--	--	False	0	0	7	17.94	0.1721	0.1721	17.89	2.693	0.0018842678	0.0018842678	False	0.0018842678	2.725	True	0.008052563	2.094	1.0	-0.0	9087	14.654999732971191	51	--	--	18126124	False	False	20.2333810891525	0.025585260300380372
797144	6.1578	4.53175	6.1578	4.53175	0	0	1983.0	15.4	1	1	--	2554972969471402624	6.15776437772795	4.5317318434594664	3.3212082386016846	2.421332359313965	--	--	2016.0	20.5	18.22	16.62	0.01398	0.0425	0.04462	--	--	--	False	0	0	7	18.45	0.1436	0.1436	18.45	3.055	0.0013964162	0.0013964162	False	0.0013964162	2.855	True	0.006583949	2.182	1.0	-0.0	9108	14.854999542236328	50	--	--	6844910	False	False	20.52842661896111	0.025585260300380372
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
895354584	265.46998	-24.70229	265.46998	-24.70229	0	0	1987.0	12.61	1	6	--	4068200586367730944	265.4699955669825	-24.702334281418175	0.12390629947185516	0.09086712449789047	-1.082	-4.547	2016.0	17.52	15.49	14.39	0.004505	0.0233	0.0099	--	--	--	False	0	0	7	15.805	0.0868	0.1674	15.65	3.04	0.0015666754	0.0013402438	False	0.0015666754	2.805	True	0.007106245	2.148	0.9999824	7.6e-06	27041651	12.029999732971191	28	--	--	13006272	False	False	17.645037422319973	0.025585260300380372
898642872	274.42264	-29.13874	274.42264	-29.13874	0	0	1987.0	12.71	1	6	--	4049788611227556864	274.4226541309113	-29.138864614211403	0.5207995176315308	0.4581608772277832	--	--	2016.0	17.05	13.4	13.15	0.01315	0.01108	0.007477	--	--	--	False	0	0	7	14.53	0.451	0.451	14.05	1.339	0.0040144934	0.0040144934	False	0.0040144934	2.396	True	0.013515946	1.869	0.99960977	0.0001695	29085609	12.130000114440918	40	--	--	7684488	False	False	17.39643253259858	0.025585260300380372
899554776	281.17974	-25.35941	281.17974	-25.35941	0	0	1987.0	13.34	1	6	--	4073399588285893888	281.1797922846329	-25.35941963297792	0.1266280561685562	0.11990871280431747	-3.402	-7.766	2016.0	17.48	15.29	15.03	0.010475	0.01608	0.005238	--	--	--	False	0	0	7	15.94	0.3293	0.1736	15.516	2.178	0.0019428418	0.0034959507	False	0.0019428418	2.71	True	0.008155835	2.088	0.9993285	0.0002918	29533494	12.760000228881836	19	--	--	6870744	False	False	17.831651188373364	0.025585260300380372
1093934200	261.95326	-46.92814	261.95326	-46.92814	0	0	1987.0	12.84	1	6	--	5951803883085737216	261.953247423593	-46.92819910314982	0.7768323421478271	0.5084171891212463	--	--	2016.0	17.02	13.55	13.77	0.01385	0.01289	0.006016	--	--	--	False	0	0	7	14.78	0.215	0.215	14.164	1.324	0.013012245	0.013012245	False	0.013012245	1.886	True	0.027914647	1.554	0.9972403	0.0012	38347099	12.260000228881836	41	--	--	11417706	False	False	17.513507985341267	0.025585260300380372
1124871688	135.30385	-52.96399	135.30385	-52.96399	0	0	1977.0	13.76	1	0	--	5323472192280364672	135.30324978626697	-52.96378585241954	0.1181718111038208	0.1138806864619255	-2.172	2.363	2016.0	18.1	14.45	14.914	0.00293	0.1385	0.774	--	--	--	False	1	0	7	15.82	1.376	1.495	13.9	0.1365	0.0043408214	0.0027803017	False	0.0043408214	2.363	True	0.014138967	1.85	0.8890504	0.0511	38930720	12.606000900268555	54	--	--	15931639	False	False	18.67119430620698	0.025585260300380372
1143735664	256.56151	-52.50011	256.56151	-52.50011	0	0	1987.0	12.56	1	6	--	5924318978646258304	256.56150138684734	-52.50016118498638	0.17903579771518707	0.15503396093845367	-5.71	-6.293	2016.0	16.94	14.12	14.13	0.008606	0.01409	0.005505	--	--	--	False	0	0	7	15.055	0.1469	0.1852	14.55	1.986	0.0048744986	0.004423809	False	0.0048744986	2.312	True	0.015454233	1.811	0.99994886	2.223e-05	42242782	11.980000495910645	40	--	--	9712869	False	False	17.3675065523912	0.025585260300380372
1165230448	80.59997	-66.08427	80.59997	-66.08427	0	0	1975.0	13.45	1	0	--	4660544569975324672	80.60002386779149	-66.0844124114311	0.33361029624938965	0.5462325811386108	1.982	-0.3604	2016.0	17.06	14.93	15.03	0.003304	0.00395	0.0007353	--	--	--	False	1	0	7	15.67	0.4978	0.5186	14.64	1.1875	0.004151973	0.0036087306	False	0.004151973	2.38	True	0.0137928445	1.86	0.91403735	0.03903	42650388	12.296000480651855	54	--	--	6766439	False	False	17.519733481073022	0.025585260300380372
1173755360	152.55394	-67.00502	152.55394	-67.00502	0	0	1987.0	12.34	1	6	--	5244972387821935104	152.55396855330667	-67.00501883010408	1.727086067199707	2.1844708919525146	--	--	2016.0	18.28	15.32	13.695	0.02032	0.02293	0.00853	--	--	--	False	0	0	7	15.766	0.04037	0.04037	15.69	3.346	0.0012363117	0.0012363117	False	0.0012363117	2.908	True	0.0061075906	2.215	1.0	0.0	42938544	11.760000228881836	47	--	--	2700905	False	False	18.303914143331113	0.025585260300380372
1183201912	238.04446	-62.51303	238.04446	-62.51303	0	0	1988.0	12.36	1	6	--	5832448803513250048	238.04445937082093	-62.51311478793788	0.19261832535266876	0.17968599498271942	-5.9	-4.03	2016.0	17.19	14.336	14.586	0.00578	0.02518	0.0217	--	--	--	False	0	0	7	15.37	0.2527	0.3052	14.77	2.414	0.002004754	0.0016109607	False	0.002004754	2.697	True	0.008368163	2.078	0.9564528	0.01933	45628466	11.779999732971191	40	--	--	1921302	False	False	17.69314285614928	0.025585260300380372
1211240168	197.19732	-67.64373	197.19732	-67.64373	0	0	1987.0	10.53	1	6	--	5857187746464509696	197.19711332414371	-67.64376190572705	0.8254939317703247	0.7541971802711487	--	--	2016.0	17.22	11.88	11.87	0.02086	0.02484	0.01767	--	--	--	False	0	0	7	13.66	0.3054	0.3054	13.0	2.47	0.0015038083	0.0015038083	False	0.0015038083	2.822	True	0.0069091073	2.16	1.0	0.0	46477408	9.949999809265137	43	--	--	7019467	True	True	17.63900121571831	0.025585260300380372

agasc_id = 549462880
agasc_gaia_matches[agasc_gaia_matches["agasc_id"] == agasc_id][
    'mag_aca_pred', 'mag_aca', 'mag_aca_pred_err', 'mag_aca_err', 'mag_aca_obs', 'mag_aca_err_obs', 'g_mag', 'bp_mag', 'rp_mag'
].pprint(max_width=-1)
report = star_report.Report(agasc_id)
report.show_in_notebook()

   mag_aca_pred        mag_aca         mag_aca_pred_err   mag_aca_err mag_aca_obs mag_aca_err_obs g_mag bp_mag rp_mag
----------------- ------------------ -------------------- ----------- ----------- --------------- ----- ------ ------
20.13661647019506 14.846000671386719 0.025585260300380372          63          --              -- 20.06  17.08  15.72

AGASC entry

Table length=1

agasc_id	ra	dec	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	epoch	tycho_id	gsc2.3
int32	float64	float64	float32	float32	float32	int16	int16	float32	str12	str10
549462880	242.79292	63.95485	16.00	14.85	--	-9999	-9999	1983.0	--	N4M8000079

Matches

Table length=1

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
549462880	1629627694387723136	16.00	17.70	1.70	0.06	0.0088	0.0225	1.000	--	--	20.06	True

Match Candidates

Table length=1

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
549462880	1629627694387723136	16.00	17.70	1.70	0.06	0.0088	0.0225	1.000	--	--	20.06	True

AGASC Neighborhood (within 20 arcsec)

Table length=1

agasc_id	ra	dec	d2d	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	pos_catid	pm_catid	mag_catid
			arcsec
int32	float64	float64	float64	float32	float32	float32	int16	int16	uint8	uint8	uint8
549462880	242.79292	63.95485	0.00	16.00	14.85	--	-9999	-9999	6	0	1

Observed Magnitude Outliers

observed = agasc_gaia_matches[~agasc_gaia_matches["mag_aca_obs"].mask]

gaia_model = gm.get_gaia_model()
observed["mag_aca_pred"] = gaia_model.predict(observed.to_pandas())
observed["d_mag_aca"] = observed["mag_aca_obs"] - observed["mag_aca"]
observed.sort("d_mag_aca", reverse=True)

observed_outliers = observed[np.abs(observed["mag_aca_obs"] - observed["mag_aca_pred"]) > 1.]
#observed_outliers[["agasc_id", "mag_aca_obs", "mag_aca", "d_mag_aca"]].pprint(max_lines=-1)

plt.hist(
    observed["mag_aca_obs"] - observed["mag_aca_pred"],
    histtype="step",
    bins=np.linspace(-5, 5, 301),
)
plt.yscale("log")
# plt.xlim(-2, 2)

with open("mag_aca_obs_outliers.txt", "w") as fh:
    for agasc_id in observed_outliers["agasc_id"]:
        fh.write(f"{agasc_id}\n")

proc = subprocess.run(
    [
        "agasc-update-magnitudes",
        "--report",
        "--agasc-id-file=mag_aca_obs_outliers.txt",
        f"--output-dir={REPORTS_DIR / 'mag_aca_obs-outliers' / 'supplement_files'}",
        f"--reports-dir={REPORTS_DIR / 'mag_aca_obs-outliers'}",
        "--log-level=info"
    ],
    capture_output=True
)

reports_dir = REPORTS_DIR / 'mag_aca_obs-outliers'
if (reports_dir / "latest").exists():
    name = (reports_dir / "latest").resolve()
    name.rename(name.parent / "supplement_reports")
    (reports_dir / "latest").unlink()

np.random.seed(549462880)
if generate_reports:
    print(f"Making report with {len(observed_outliers)} stars")

    description = """
    <p>This list includes the stars with abs(mag_aca_1p8 - mag_aca_obs) > 1 </p>

    <p>
    The purpose of this report is to check whether outliers can be caused by misidentification.
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>

    <h2> Notable examples:</h2>
    <ul>
      <li></li>

    </ul>

    <h2> Stars </h2>
    """

    reports_dir = REPORTS_DIR / "mag_aca_obs-outliers" / "agasc_gaia_reports"
    star_report.make_report_list(
        data=observed_outliers[['agasc_id', 'p_value', 'd2d', 'd_mag', 'd_mag_aca', 'mag_aca', 'mag_aca_pred', 'mag_aca_obs']],
        path=reports_dir,
        title='Catalog Magnitude Outliers',
        description=description,
        overwrite=True
    )

Making report with 76 stars

76it [00:16,  4.69it/s]

"Difficult" Stars

There is a very important point to make

"Difficult" stars are grouped with all the stars it could be confused with. The matches in each group are recomputed to guarantee there are no two repeated AGASC or Gaia IDs. The process starts by selecting the match with the highest probability, removing the corresponding AGASC and Gaia IDs from the candidate matches, and repeating until there are no candidate matches left.

For each of these groups, we defined "latest_pos_cat" as the POS_CATID value with the highest precedence from the AGASC entries in the group. The precedence is (in decreasing order) [5, 6, 4, 3, 2, 1]. The catalog precedence IS NOT considered when recomputing the matches.

The entries with POS_CATID different than latest_pos_cat can be considered duplicates, but this is not guaranteed to be the case.

For example, AGASC 102499594 and 102499593 are two stars in Tycho2 that form a binary system (?). AGASC 102499594 is a star in GSC2.3 that lies right in between, and has a magnitude consistent with it being a blend of the other two. Based on catalog precedence, 102499594 appears to be a duplicate. Gaia IDs 3151414218873077760 and 3151414218874713600 are two resolved stars in Gaia that are matched to AGASC 102499594 and 102499593 respectively. AGASC 102499594 is matched to the nearest star, which is 13 arcsec away, and this match is discarded as background based on p-value.

One can ask the question of whether there are any possible duplicates that ended up being the best match. The answer is no:

agasc_difficult[
    (agasc_difficult['pos_catid'] != agasc_difficult['latest_pos_cat'])
    & (agasc_difficult['best_match'])
]

Table length=0

agasc_id	gaia_id	best_match	d2d	d_mag	p_match	p_value	p_relative	idx	mag	mag_band	mag_catid	pos_catid	best_match_0	latest_pos_cat	best_match_1	best_match_2	group	group_size
int64	int64	bool	float16	float16	float32	float32	float32	int64	float32	int16	uint8	uint8	bool	int64	bool	bool	int64	int64

# AGASC_GAIA_X_MATCH_ALL = star_report.AGASC_GAIA_X_MATCH_ALL
# SUMMARY = star_report.SUMMARY
# import importlib
# importlib.reload(star_report)
# star_report.AGASC_GAIA_X_MATCH_ALL = AGASC_GAIA_X_MATCH_ALL
# star_report.SUMMARY = SUMMARY

report = star_report.Report([
    181023536,
    181023537,
    181023560,
    181023552,
    181023544,
    181023568
])
report.show_in_notebook()
report.save_html(REPORTS_DIR / "difficult")

AGASC 181023536 epoch in AGASC: 2000.0, epoch in summary: 1991.74
AGASC 181023537 epoch in AGASC: 2000.0, epoch in summary: 1991.74
AGASC 181023560 epoch in AGASC: 2000.0, epoch in summary: 1991.87

AGASC entry

Table length=6

agasc_id	ra	dec	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	epoch	tycho_id	gsc2.3
int32	float64	float64	float32	float32	float32	int16	int16	float32	str12	str10
181023544	123.05355	17.648	5.59	4.89	--	74	-119	2000.0	--	N908000475
181023537	123.05319049	17.64774942	6.24	6.23	--	80	-129	2000.0	1381-1638-2	0.0
181023536	123.05319049	17.64774942	5.30	5.36	--	80	-129	2000.0	1381-1638-1	N908000476
181023552	123.05312	17.64772	5.19	4.48	--	-9999	-9999	2000.0	--	N908000477
181023560	123.05493923	17.64823567	5.85	5.87	--	86	-91	2000.0	1381-1641-1	N908000478
181023568	123.05483	17.64742	6.30	5.59	--	-9999	-9999	2000.0	--	N908000479

Matches

Table length=3

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
181023536	657244521593509376	5.30	5.56	0.26	3.74	0.2341	0.1796	0.949	--	--	5.867	True
181023560	657244521593508992	5.85	6.33	0.48	1.53	0.0003	0.0026	1.000	37.000	-150.875	5.934	True
181023537	657244586015485440	6.24	5.51	-0.73	1.38	0.0002	0.0018	0.175	118.312	-148.375	5.52	True

Match Candidates

Table length=18

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
181023552	657244586015485440	5.19	5.33	0.14	0.49	0.3338	0.2249	1.000	118.312	-148.375	5.52	False
181023536	657244521593509376	5.30	5.56	0.26	3.74	0.2341	0.1796	0.949	--	--	5.867	True
181023536	657244586015485440	5.30	5.51	0.21	0.98	0.0126	0.0275	0.051	118.312	-148.375	5.52	False
181023544	657244521593509376	5.59	5.38	-0.21	2.62	0.0100	0.0241	0.984	--	--	5.867	False
181023537	657244521593509376	6.24	5.56	-0.69	2.84	0.0010	0.0051	0.825	--	--	5.867	False
181023560	657244521593508992	5.85	6.33	0.48	1.53	0.0003	0.0026	1.000	37.000	-150.875	5.934	True
181023537	657244586015485440	6.24	5.51	-0.73	1.38	0.0002	0.0018	0.175	118.312	-148.375	5.52	True
181023552	657244521593509376	5.19	5.38	0.19	2.43	0.0002	0.0015	0.000	--	--	5.867	False
181023544	657244586015485440	5.59	5.33	-0.26	2.28	0.0002	0.0015	0.016	118.312	-148.375	5.52	False
181023568	657244521593508992	6.30	6.13	-0.17	3.90	0.0000	0.0001	0.981	37.000	-150.875	5.934	False
181023568	657244521593509376	6.30	5.38	-0.92	4.05	0.0000	0.0000	0.019	--	--	5.867	False
181023560	657244521593509376	5.85	5.56	-0.29	5.57	0.0000	0.0000	0.000	--	--	5.867	False
181023544	657244521593508992	5.59	6.13	0.54	5.27	0.0000	0.0000	0.000	37.000	-150.875	5.934	False
181023537	657244521593508992	6.24	6.33	0.08	7.12	0.0000	0.0000	0.000	37.000	-150.875	5.934	False
181023568	657244586015485440	6.30	5.33	-0.97	6.32	0.0000	0.0000	0.000	118.312	-148.375	5.52	False
181023560	657244586015485440	5.85	5.51	-0.34	7.11	0.0000	0.0000	0.000	118.312	-148.375	5.52	False
181023552	657244521593508992	5.19	6.13	0.94	7.02	0.0000	0.0000	0.000	37.000	-150.875	5.934	False
181023536	657244521593508992	5.30	6.33	1.03	7.14	0.0000	0.0000	0.000	37.000	-150.875	5.934	False

AGASC Neighborhood (within 20 arcsec)

Table length=6

agasc_id	ra	dec	d2d	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	pos_catid	pm_catid	mag_catid
			arcsec
int32	float64	float64	float64	float32	float32	float32	int16	int16	uint8	uint8	uint8
181023544	123.05355	17.648	1.17	5.59	4.89	--	74	-119	2	2	1
181023537	123.05319049	17.64774942	2.14	6.24	6.23	--	80	-129	5	5	5
181023536	123.05319049	17.64774942	2.14	5.30	5.36	--	80	-129	5	5	5
181023552	123.05312	17.64772	2.18	5.19	4.48	--	-9999	-9999	6	0	1
181023560	123.05493923	17.64823567	4.11	5.85	5.87	--	86	-91	5	5	5
181023568	123.05483	17.64742	4.22	6.30	5.59	--	-9999	-9999	6	0	1

report = star_report.Report([
    120468440, 120468441, 120468448, 120468456, 120468464
])
report.show_in_notebook()

AGASC 120468440 epoch in AGASC: 2000.0, epoch in summary: 1991.59
AGASC 120468441 epoch in AGASC: 2000.0, epoch in summary: 1991.55

AGASC entry

Table length=5

agasc_id	ra	dec	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	epoch	tycho_id	gsc2.3
int32	float64	float64	float32	float32	float32	int16	int16	float32	str12	str10
120468456	226.8868	9.2253	7.40	6.67	--	-205	7	2000.0	--	N5RX000015
120468441	226.88670436	9.22505169	7.54	7.48	--	-189	25	2000.0	919-1659-2	0.0
120468448	226.88646	9.22542	6.80	6.09	--	-9999	-9999	2000.0	--	N5RX000014
120468464	226.88658	9.22494	7.50	6.79	--	-9999	-9999	2000.0	--	N5RX000016
120468440	226.88728703	9.22604991	7.35	7.32	--	-190	29	2000.0	919-1659-1	N5RX000013

Matches

Table length=2

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
120468440	1167718049194619264	7.35	7.55	0.20	0.03	1.4049	0.5671	1.000	-194.875	22.922	7.188	True
120468441	1167718049194619136	7.54	7.74	0.20	0.02	1.3548	0.5534	1.000	-190.500	30.391	7.38	True

Match Candidates

Table length=10

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
120468440	1167718049194619264	7.35	7.55	0.20	0.03	1.4049	0.5671	1.000	-194.875	22.922	7.188	True
120468441	1167718049194619136	7.54	7.74	0.20	0.02	1.3548	0.5534	1.000	-190.500	30.391	7.38	True
120468464	1167718049194619136	7.50	7.52	0.02	0.58	0.2970	0.2097	1.000	-190.500	30.391	7.38	False
120468456	1167718049194619136	7.40	7.52	0.12	0.94	0.0250	0.0433	0.997	-190.500	30.391	7.38	False
120468448	1167718049194619136	6.80	7.52	0.72	1.50	0.0001	0.0013	0.989	-190.500	30.391	7.38	False
120468456	1167718049194619264	7.40	7.33	-0.07	3.13	0.0001	0.0008	0.003	-194.875	22.922	7.188	False
120468441	1167718049194619264	7.54	7.55	0.01	4.14	0.0000	0.0001	0.000	-194.875	22.922	7.188	False
120468448	1167718049194619264	6.80	7.33	0.53	3.64	0.0000	0.0000	0.011	-194.875	22.922	7.188	False
120468440	1167718049194619136	7.35	7.74	0.39	4.12	0.0000	0.0000	0.000	-190.500	30.391	7.38	False
120468464	1167718049194619264	7.50	7.33	-0.17	4.64	0.0000	0.0000	0.000	-194.875	22.922	7.188	False

AGASC Neighborhood (within 20 arcsec)

Table length=5

agasc_id	ra	dec	d2d	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	pos_catid	pm_catid	mag_catid
			arcsec
int32	float64	float64	float64	float32	float32	float32	int16	int16	uint8	uint8	uint8
120468456	226.8868	9.2253	0.22	7.40	6.67	--	-205	7	2	2	1
120468441	226.88670436	9.22505169	1.09	7.54	7.48	--	-189	25	5	5	5
120468448	226.88646	9.22542	1.74	6.80	6.09	--	-9999	-9999	6	0	1
120468464	226.88658	9.22494	1.89	7.50	6.79	--	-9999	-9999	6	0	1
120468440	226.88728703	9.22604991	3.13	7.35	7.32	--	-190	29	5	5	5

report = star_report.Report([102499592, 102499593, 102499594])
report.show_in_notebook()

AGASC 102499593 epoch in AGASC: 2000.0, epoch in summary: 1991.90
AGASC 102499594 epoch in AGASC: 2000.0, epoch in summary: 1991.90

AGASC entry

Table length=3

agasc_id	ra	dec	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	epoch	tycho_id	gsc2.3
int32	float64	float64	float32	float32	float32	int16	int16	float32	str12	str10
102499592	116.68226885	11.09789718	10.57	10.72	--	-9999	-9999	1991.25	--	N88Z000037
102499593	116.68230917	11.09787733	11.08	11.22	--	-9999	-9999	2000.0	782-161-2	0.0
102499594	116.68230917	11.09787733	11.09	11.41	--	-9999	-9999	2000.0	782-161-3	0.0

Matches

Table length=3

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
102499593	3151414218873077760	11.08	11.06	-0.02	0.09	3.0041	0.9532	1.000	-0.856	0.586	11.09	True
102499594	3151414218874713600	11.09	10.87	-0.22	0.04	1.2962	0.5388	0.999	-0.781	0.512	10.94	True
102499592	3151414214578577280	10.57	17.81	7.25	13.02	0.0000	0.0000	0.000	3.031	3.064	17.69	True

Match Candidates

Table length=10

agasc_id	gaia_id	mag_1p7	mag_pred	d_mag	d2d	p_match	p_value	p_relative	pm_ra	pm_dec	g_mag	baseline
int64	int64	float32	float16	float16	float16	float32	float32	float32	float16	float16	float16	bool
102499593	3151414218873077760	11.08	11.06	-0.02	0.09	3.0041	0.9532	1.000	-0.856	0.586	11.09	True
102499594	3151414218874713600	11.09	10.87	-0.22	0.04	1.2962	0.5388	0.999	-0.781	0.512	10.94	True
102499592	3151414218874713600	10.57	10.87	0.30	0.95	0.0084	0.0221	0.743	-0.781	0.512	10.94	False
102499592	3151414218873077760	10.57	11.06	0.49	0.95	0.0029	0.0107	0.257	-0.856	0.586	11.09	False
102499594	3151414218873077760	11.09	11.06	-0.03	1.92	0.0012	0.0058	0.001	-0.856	0.586	11.09	False
102499593	3151414218874713600	11.08	10.87	-0.21	1.97	0.0004	0.0030	0.000	-0.781	0.512	10.94	False
102499594	3151414214578577280	11.09	17.81	6.73	12.28	0.0000	0.0000	0.000	3.031	3.064	17.69	False
102499592	3151414214578577280	10.57	17.81	7.25	13.02	0.0000	0.0000	0.000	3.031	3.064	17.69	True
102499593	3151414214578577280	11.08	17.81	6.74	13.85	0.0000	0.0000	0.000	3.031	3.064	17.69	False
102499593	3151414283303787648	11.08	18.08	6.99	14.05	0.0000	0.0000	0.000	-1.754	-7.258	17.94	False

AGASC Neighborhood (within 20 arcsec)

Table length=3

agasc_id	ra	dec	d2d	mag	mag_aca	mag_aca_obs	pm_ra	pm_dec	pos_catid	pm_catid	mag_catid
			arcsec
int32	float64	float64	float64	float32	float32	float32	int16	int16	uint8	uint8	uint8
102499592	116.68226885	11.09789718	0.02	10.57	10.72	--	-9999	-9999	5	0	5
102499593	116.68230917	11.09787733	0.18	11.08	11.22	--	-9999	-9999	5	0	5
102499594	116.68230917	11.09787733	0.18	11.09	11.41	--	-9999	-9999	5	0	5

np.random.seed(120468440)

diff_1 = join(agasc_difficult, agasc_gaia_matches, keys=['agasc_id', 'gaia_id'], join_type='left', table_names=['difficult', 'base'])
diff_1.rename_columns(
    [col for col in diff_1.colnames if col[-5:] == '_base'],
    [col[:-5] for col in diff_1.colnames if col[-5:] == '_base'],
)
diff_1['best_match'][diff_1['best_match'].mask] = False

diff_1['d_mag'].format = "{:.2f}"
diff_1['mag_aca'].format = "{:.2f}"
diff_1['mag_aca_pred'].format = "{:.2f}"
diff_1['mag_aca_obs'].format = "{:.2f}"
diff_1['d2d'].format = "{:.3f}"
diff_1['p_value'].format = "{:.4f}"

diff_1.sort(['group', 'mag_aca'])
diff_1 = diff_1.group_by('group')

hand_picked_ids = np.asarray(diff_1['group'][np.in1d(diff_1['agasc_id'], [11931168, 980817608, 5253584, 77996162, 23877032, 102499593, 102499592])])
group_idx = np.unique(np.concatenate([
    np.random.choice(np.unique(diff_1['group']), size=50),
    np.asarray(np.unique(diff_1['group'][diff_1['group_size'] >= 3])[:50]),
    hand_picked_ids
]))
diff = diff_1.groups[group_idx]

generate_reports = True
if generate_reports:
    print(f"Making report with {len(diff.groups)} groups of stars")

    description = """
    <p>This list includes "difficult" stars. </p>

    <p>
    Difficult stars are the ones that had matching collisions. In other words, the same Gaia star
    was matched to more than one AGASC star. Difficult stars are divided in equivalence classes,
    based on which stars it collided. Often these stars are duplicates in AGASC, but not always.
    </p>

    <p>
    The purpose of this report is to check whether outliers can be caused by misidentification.
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>

    <h2> Notable examples:</h2>
    <ul>
      <li>
        <a href="report_181023536_181023537_181023560_181023552_181023544_181023568.html"> 181023536 et al </a>.
        The most difficult one. Six AGASC stars and three Gaia stars. Some AGASC stars might be duplicates.
        Not all Gaia stars have proper motion.
      </li>
      <li>
        <a href="report_11931168.html"> 11931168 </a>. Two stars in AGASC matched to two in Gaia, but
        they are difficult. The table values are masked. NEED TO CHECK.
      </li>
      <li>
        <a href="report_980817608.html"> 980817608 </a>. Three stars in AGASC, but only one in Gaia.
        Two of the AGASC stars are from GSC2.3 and are separated by 6 arcsec. The other star is a
        Tycho2 star right in between. The Gaia star matches the Tycho2 star. Sounds like a spurious
        binary system identification in GSC2.3, which should correspond to the single Tycho star.
      </li>
      <li>
        <a href="report_5253584.html"> 5253584 </a>. Three stars in AGASC matched to one in Gaia.
        Two of the three are Tycho2 stars at the exact same location. The third is an older catalog.
      </li>
      <li>
        <a href="report_77996162.html"> 77996162 </a>. Four stars in AGASC matched to two in Gaia.
      </li>
      <li>
        <a href="report_23877032.html"> 23877032 </a>. Three stars in AGASC matched to one in Gaia.
        Two are duplicates.
      </li>
      <li>
        <a href="report_102499593.html"> 102499593 </a>. Two stars in AGASC matched to three in Gaia.
        
      </li>

      [102499592, 102499593, 102499594]. One star is a duplicate. Fortunately, this is the one that
      is matched to a third Gaia star, and is removed based on p-value. It could have happened that
      the duplicate gets a better p-value than one of the other two, so one of the other two is not
      updated, and the duplicate is updated instead. Duplicates must be removed from the match candidates.
    </ul>

    <h2> Stars </h2>
    """

    reports_dir = REPORTS_DIR / "difficult"
    star_report.make_report_list_by_group(
        data=diff.groups[:][['agasc_id', 'best_match', 'group_size', 'group']],
        path=reports_dir,
        title='Difficult Stars',
        description=description,
        overwrite=True
    )

Making report with 102 groups of stars

102it [00:24,  4.09it/s]

Candidates with low p-value

bins = np.linspace(0, 1, 101)
sns.histplot(
    agasc_gaia_matches['p_value'],
    stat="density",
    bins=bins,
    label="all",
)
sns.histplot(
    agasc_gaia_matches['p_value'][agasc_gaia_matches['guide'] | agasc_gaia_matches['acq']],
    stat="density",
    bins=bins,
    label="candidates",
)
plt.legend()
plt.xlabel('p-value')
plt.title("p-value distribution (candidates)")
plt.yscale('log')

n_cand = np.count_nonzero((agasc_gaia_matches['guide'] | agasc_gaia_matches['acq']))
n_cand_marginal = np.count_nonzero((agasc_gaia_matches['p_value'] < 0.02) & (agasc_gaia_matches['guide'] | agasc_gaia_matches['acq']))

print(f"{n_cand_marginal} marginal candidates out of {n_cand}")

443208

candidates = agasc_gaia_matches[
    (agasc_gaia_matches['p_value'] < 0.022) & (agasc_gaia_matches['guide'] | agasc_gaia_matches['acq'])
]
candidates.sort('p_value')
candidates = candidates[:100]

np.random.seed(120468464)
generate_reports = True
if generate_reports:
    print(f"Making report with {len(candidates)} stars")

    description = """
    <p>This list includes stars that are guide or acq candidates with low p-value. </p>

    <p>
    The purpose of this report is to check whether outliers can be caused by misidentification.
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>

    <h2> Notable examples:</h2>
    <ul>
      <li>
        <a href="report_303961568.html"> 303961568 </a>. Gaia star is missing from plot. NEED TO CHECK.
      </li>
    </ul>

    <h2> Stars </h2>
    """

    reports_dir = "/Users/javierg/SAO/Notebooks/agasc/gaia-magnitudes-2023/reports/candidates-p-value-0.01-0.022"
    star_report.make_report_list(
        data=candidates[['agasc_id', 'p_value', 'd2d', 'd_mag', 'mag_aca', 'mag_aca_pred', 'mag_aca_obs']],
        path=reports_dir,
        title='Candidates with low p-value',
        description=description,
        overwrite=True
    )

Making report with 100 stars

100it [03:48,  2.28s/it]

agasc_gaia_matches_all = xm.get_agasc_gaia_x_match_difficult_fixed()
agasc_gaia_matches_all = agasc_gaia_matches_all[agasc_gaia_matches_all['best_match']]

i = np.searchsorted(agasc_summary["agasc_id"], agasc_gaia_matches_all["agasc_id"])
cols = [
    "mag_aca",
    "mag_aca_err",
    "mag_aca_obs",
    "mag_aca_err_obs",
    "mag_catid",
    "mag_band",
    "random_index",
    "guide",
    "acq"
]
for col in cols:
    agasc_gaia_matches_all[col] = agasc_summary[col][i]

cols = [
    "pm_ra",
    "pm_dec",
    "epoch",
    "ra",
    "dec",
]
for col in cols:
    agasc_gaia_matches_all[f"{col}_agasc"] = agasc_summary[col][i]

bins = np.linspace(0, 1, 101)
sns.histplot(
    agasc_gaia_matches_all['p_value'],
    stat="density",
    bins=bins,
    label="all",
)
sns.histplot(
    agasc_gaia_matches_all['p_value'][agasc_gaia_matches_all['guide'] | agasc_gaia_matches_all['acq']],
    stat="density",
    bins=bins,
    label="candidates",
)
plt.legend()
plt.xlabel('p-value')
plt.title("p-value distribution (candidates)")
plt.yscale('log')

n_cand = np.count_nonzero((agasc_gaia_matches_all['guide'] | agasc_gaia_matches_all['acq']))
n_cand_marginal = np.count_nonzero((agasc_gaia_matches_all['p_value'] < 0.02) & (agasc_gaia_matches_all['guide'] | agasc_gaia_matches_all['acq']))

print(f"{n_cand_marginal} marginal candidates out of {n_cand}")

1633 marginal candidates out of 444377

np.count_nonzero(
    (agasc_gaia_matches_all['guide'] | agasc_gaia_matches_all['acq'])
    & (agasc_gaia_matches_all["p_value"] < 0.005)
    & (agasc_gaia_matches_all["p_value"] > 0.003)
)

157

bad_candidates = agasc_gaia_matches_all[
    (agasc_gaia_matches_all['guide'] | agasc_gaia_matches_all['acq'])
    & (agasc_gaia_matches_all["p_value"] < 0.005)
    & (agasc_gaia_matches_all["p_value"] > 0.003)
]

np.random.seed(120468464)
generate_reports = True
if generate_reports:
    print(f"Making report with {len(bad_candidates)} stars")

    description = """
    <p>This list includes stars that are guide or acq candidates with p-value between 0.003 and 0.005. </p>

    <p>
    The purpose of this report is to check whether outliers can be caused by misidentification.
    When looking at the report for a given outlier, consider that if the true match is a star with no
    proper motion in Gaia, then there should be an AGASC star matched to this Gaia star. Is there one?
    </p>

    <p>If you double-click on a report's figure, it will zoom out and show you all AGASC stars around.</p>

    <h2> Notable examples:</h2>
    <ul>
      <li>
        <a href="report_303961568.html"> 303961568 </a>. Gaia star is missing from plot. NEED TO CHECK.
      </li>
    </ul>

    <h2> Stars </h2>
    """

    reports_dir = "/Users/javierg/SAO/Notebooks/agasc/gaia-magnitudes-2023/reports/candidates-p-value-0.003-0.005"
    star_report.make_report_list(
        data=bad_candidates[['agasc_id', 'p_value', 'd2d', 'd_mag', 'mag_aca', 'mag_aca_obs']],
        path=reports_dir,
        title='Candidates with low p-value',
        description=description,
        overwrite=True
    )

Making report with 157 stars

157it [06:47,  2.60s/it]