Skip to content

Commit

Permalink
Merge branch 'develop'
Browse files Browse the repository at this point in the history
  • Loading branch information
@Gabriel-p
Gabriel-p committed Jul 31, 2024
2 parents 09204a9 + ba26ec8 commit 70c0781
Show file tree
Hide file tree
Showing 24 changed files with 1,617 additions and 966 deletions.
176 changes: 128 additions & 48 deletions asteca/cluster.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,59 +119,61 @@ def _load(self):
dim_count = 0

if self.ra is not None:
self.ra_v = self.obs_df[self.ra].values
self.dec_v = self.obs_df[self.dec].values
self.ra_v = np.array(self.obs_df[self.ra], dtype=float)
self.dec_v = np.array(self.obs_df[self.dec], dtype=float)
print(f"(RA, DEC) : ({self.ra}, {self.dec})")
dim_count += 1

# TODO: change _p to _v
if self.magnitude is not None:
self.mag_p = self.obs_df[self.magnitude].values
self.mag_p = np.array(self.obs_df[self.magnitude], dtype=float)
if self.e_mag is not None:
self.e_mag_p = self.obs_df[self.e_mag].values
self.e_mag_p = np.array(self.obs_df[self.e_mag], dtype=float)
else:
raise ValueError("Magnitude uncertainty is required")
print(f"Magnitude : {self.magnitude} [{self.e_mag}]")
dim_count += 1

if self.color is not None:
self.colors_p = [self.obs_df[self.color].values]
self.colors_p = [np.array(self.obs_df[self.color], dtype=float)]
if self.e_color is not None:
self.e_colors_p = [self.obs_df[self.e_color].values]
self.e_colors_p = [np.array(self.obs_df[self.e_color], dtype=float)]
else:
raise ValueError("Color uncertainty is required")
print(f"Color : {self.color} [{self.e_color}]")
dim_count += 1
if self.color2 is not None:
self.colors_p.append(self.obs_df[self.color2].values)
self.colors_p.append(np.array(self.obs_df[self.color2], dtype=float))
if self.e_color2 is not None:
self.e_colors_p.append(self.obs_df[self.e_color2].values)
self.e_colors_p.append(
np.array(self.obs_df[self.e_color2], dtype=float)
)
else:
raise ValueError("Color2 uncertainty is required")
print(f"Color2 : {self.color2} [{self.e_color2}]")
dim_count += 1

if self.plx is not None:
self.plx_v = self.obs_df[self.plx].values
self.plx_v = np.array(self.obs_df[self.plx], dtype=float)
if self.e_plx is not None:
self.e_plx_v = self.obs_df[self.e_plx].values
self.e_plx_v = np.array(self.obs_df[self.e_plx], dtype=float)
else:
raise ValueError("Parallax uncertainty is required")
print(f"plx : {self.plx} [{self.e_plx}]")
dim_count += 1

if self.pmra is not None:
self.pmra_v = self.obs_df[self.pmra].values
self.pmra_v = np.array(self.obs_df[self.pmra], dtype=float)
if self.e_pmra is not None:
self.e_pmra_v = self.obs_df[self.e_pmra].values
self.e_pmra_v = np.array(self.obs_df[self.e_pmra], dtype=float)
else:
raise ValueError("pmRA uncertainty is required")
print(f"pmRA : {self.pmra} [{self.e_pmra}]")
dim_count += 1
if self.pmde is not None:
self.pmde_v = self.obs_df[self.pmde].values
self.pmde_v = np.array(self.obs_df[self.pmde], dtype=float)
if self.e_pmra is not None:
self.e_pmde_v = self.obs_df[self.e_pmde].values
self.e_pmde_v = np.array(self.obs_df[self.e_pmde], dtype=float)
else:
raise ValueError("pmDE uncertainty is required")
print(f"pmDE : {self.pmra} [{self.e_pmde}]")
Expand All @@ -182,11 +184,10 @@ def _load(self):
def get_center(
self,
algo: str = "knn_5d",
radec_c: tuple | None = None,
pms_c: tuple | None = None,
data_2d: str = "radec",
radec_c: tuple[float, float] | None = None,
pms_c: tuple[float, float] | None = None,
plx_c: float | None = None,
N_cluster: int | None = None,
N_clust_min: int = 25,
) -> None:
"""Estimate center coordinates for the cluster
Expand All @@ -197,22 +198,23 @@ def get_center(
(k=N_clust_min) nearest stars to an estimate of the center in proper motions
and (ra, dec, plx), if given.
:param algo: Algorithm used to estimate center values, one of (``knn_5d``),
defaults to ``knn_5d``
``kde_2d``: Estimates the center value using a Kernel Density Estimator (KDE)
given a two dimensional array determined by the ``data_2d`` argument.
:param algo: Algorithm used to estimate center values, one of
(``knn_5d``, ``kde_2d``), defaults to ``knn_5d``
:type algo: str
:param data_2d: String indicating the data to be used to estimate
the center value, either: ``radec`` or ``pms``, defaults to ``radec``
:type data_2d: str
:param radec_c: Estimated value for the (RA, DEC) center, defaults to ``None``
:type radec_c: tuple | None
:type radec_c: tuple[float, float] | None
:param pms_c: Estimated value for the (pmRA, pmDE) center, defaults to ``None``
:type pms_c: tuple | None
:type pms_c: tuple[float, float] | None
:param plx_c: Estimated value for the plx center, defaults to ``None``
:type plx_c: float | None
:param N_cluster: Estimated number of members, defaults to ``None``
:type N_cluster: int | None
:param N_clust_min: Minimum number of cluster members, defaults to ``25``
:type N_clust_min: int
:raises ValueError: If the ``knn_5d`` algorithm is selected and any of these
attributes ``(ra, dec, pmra, pmde, plx)`` are missing from the
:raises ValueError: If required data is missing from the
:py:class:`Cluster <asteca.cluster.Cluster>` object
"""

Expand All @@ -225,8 +227,13 @@ def get_center(
+ "to be defined"
)

# To galactic coordinates
# To galactic coordinates (not optional, always use instead of equatorial)
glon, glat = cp.radec2lonlat(self.ra_v, self.dec_v)
lonlat_c = None
if radec_c is not None:
lonlat_c = cp.radec2lonlat(radec_c[0], radec_c[1])

# Remove nans
X = np.array([glon, glat, self.pmra_v, self.pmde_v, self.plx_v])
# Reject nan values and extract clean data
_, X_no_nan = cp.reject_nans(X)
Expand All @@ -238,22 +245,94 @@ def get_center(
pmRA,
pmDE,
plx,
xy_c=radec_c,
xy_c=lonlat_c,
vpd_c=pms_c,
plx_c=plx_c,
N_cluster=N_cluster,
N_clust_min=N_clust_min,
N_clust_min=self.N_clust_min,
N_clust_max=self.N_clust_max,
)
ra_c, dec_c = cp.lonlat2radec(x_c, y_c)

print("\nCenter coordinates found:")
print("radec_c : ({:.4f}, {:.4f})".format(ra_c, dec_c))
print("pms_c : ({:.3f}, {:.3f})".format(pmra_c, pmde_c))
print("plx_c : {:.3f}".format(plx_c))
print("\nCenter coordinates found:")
print("radec_c : ({:.4f}, {:.4f})".format(ra_c, dec_c))
print("pms_c : ({:.3f}, {:.3f})".format(pmra_c, pmde_c))
print("plx_c : {:.3f}".format(plx_c))

self.radec_c = [ra_c, dec_c]
self.pms_c = [pmra_c, pmde_c]
self.plx_c = plx_c

elif algo == "kde_2d":
if data_2d == "radec":
if any([_ is None for _ in (self.ra, self.dec)]):
raise ValueError("Data for (ra, dec) data is required")
c_str = "radec_c"
x, y = self.ra_v, self.dec_v
elif data_2d == "pms":
if any([_ is None for _ in (self.pmra, self.pmde)]):
raise ValueError("Data for (pmra, pmde) data is required")
c_str = "pms_c"
x, y = self.pmra_v, self.pmde_v
else:
raise ValueError(f"Algorithm '{data_2d}' argument not recognized")

# Remove rows containing any NaNs
array_2d = np.array([x, y])
x, y = array_2d[:, ~np.isnan(array_2d).any(axis=0)]

self.radec_c = [ra_c, dec_c]
self.pms_c = [pmra_c, pmde_c]
self.plx_c = plx_c
x_c, y_c = cp.get_2D_center(x, y)

print("\nCenter coordinates found:")
print("{} : ({:.4f}, {:.4f})".format(c_str, x_c, y_c))
if c_str == "radec_c":
self.radec_c = [x_c, y_c]
if c_str == "pms_c":
self.pms_c = [x_c, y_c]

else:
raise ValueError(f"Selected method '{algo}' not recognized")

# def _get_radius(self, algo: str = "field_dens") -> None:
# """Estimate the cluster radius

# - ``ripley``: Originally introduced with the ``fastMP`` membership method
# in `Perren et al.
# (2023) <https://academic.oup.com/mnras/article/526/3/4107/7276628>`__.
# Requires ``(ra, dec)`` and their center estimates.
# - ``density``: Simple algorithm that counts the number of stars within the
# cluster region (center+radius) and subtracts the expected number of field
# stars within that region. Requires the number of cluster members to be
# defined.

# :param algo: Algorithm used to estimate center values, one of
# (``ripley, density``); defaults to ``ripley``
# :type algo: str

# :raises ValueError: If required attributes are missing from the
# :py:class:`Cluster <asteca.cluster.Cluster>` object
# """

# if algo == "field_dens":
# for k in ("ra", "dec", "radec_c"):
# if hasattr(self, k) is False:
# raise ValueError(f"'{k}' must be present as a 'cluster' attribute")
# elif algo == "king":
# for k in ("ra", "dec", "radec_c"):
# if hasattr(self, k) is False:
# raise ValueError(f"'{k}' must be present as a 'cluster' attribute")

# xv, yv = self.ra_v, self.dec_v
# xy_center = self.radec_c

# if algo == "field_dens":
# radius = cp.fdens_radius(xv, yv, list(xy_center))
# elif algo == "king":
# radius = cp.king_radius(xv, yv, xy_center, self.N_cluster)
# else:
# raise ValueError(f"Unrecognized method '{algo}")

# print(f"\nRadius : {radius}")
# self.radius = radius

def get_nmembers(self, algo: str = "ripley", eq_to_gal: bool = True) -> None:
"""Estimate the number of members for the cluster
Expand All @@ -264,7 +343,8 @@ def get_nmembers(self, algo: str = "ripley", eq_to_gal: bool = True) -> None:
Requires ``(ra, dec, pmra, pmde, plx)`` and their center estimates.
- ``density``: Simple algorithm that counts the number of stars within the
cluster region (center+radius) and subtracts the expected number of field
stars within that region. Requires the radius of the cluster to be defined.
stars within that region. Requires the `(RA, DEC)` center and the radius of
the cluster to be defined.
:param algo: Algorithm used to estimate center values, one of
(``ripley, density``); defaults to ``ripley``
Expand All @@ -277,14 +357,8 @@ def get_nmembers(self, algo: str = "ripley", eq_to_gal: bool = True) -> None:
:raises ValueError: If required attributes are missing from the
:py:class:`Cluster <asteca.cluster.Cluster>` object
"""
if algo == "ripley":
for k in ("ra", "dec", "pmra", "pmde", "plx", "radec_c", "pms_c", "plx_c"):
if hasattr(self, k) is False:
raise ValueError(f"'{k}' must be present as a 'cluster' attribute")
elif algo == "density":
for k in ("ra", "dec", "radec_c", "radius"):
if hasattr(self, k) is False:
raise ValueError(f"'{k}' must be present as a 'cluster' attribute")
if algo not in ("ripley", "density"):
raise ValueError(f"'Argument algo={algo}' not recognized")

xv, yv = self.ra_v, self.dec_v
xy_center = self.radec_c
Expand All @@ -294,6 +368,9 @@ def get_nmembers(self, algo: str = "ripley", eq_to_gal: bool = True) -> None:
xy_center = cp.radec2lonlat(*xy_center)

if algo == "ripley":
for k in ("ra", "dec", "pmra", "pmde", "plx", "radec_c", "pms_c", "plx_c"):
if hasattr(self, k) is False:
raise ValueError(f"'{k}' must be present as a 'cluster' attribute")
N_cluster = nm.ripley_nmembs(
xv,
yv,
Expand All @@ -305,6 +382,9 @@ def get_nmembers(self, algo: str = "ripley", eq_to_gal: bool = True) -> None:
self.plx_c,
)
elif algo == "density":
for k in ("ra", "dec", "radec_c", "radius"):
if hasattr(self, k) is False:
raise ValueError(f"'{k}' must be present as a 'cluster' attribute")
N_cluster = nm.density_nmembs(xv, yv, xy_center, self.radius)

if N_cluster < self.N_clust_min:
Expand Down
14 changes: 7 additions & 7 deletions asteca/likelihood.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,17 +14,15 @@ class Likelihood:
:param my_cluster: :py:class:`Cluster <asteca.cluster.Cluster>` object with the
loaded data for the observed cluster
:type my_cluster: Cluster
:param lkl_name: Currently only the Poisson likelihood ratio defined in
:param lkl_name: Currently only the Poisson likelihood ratio (``plr``) defined in
`Tremmel et al. (2013)
<https://ui.adsabs.harvard.edu/abs/2013ApJ...766...19T/abstract)>`__
is accepted, defaults to ``plr``
:type lkl_name: str
:param bin_method: Bin method used to split the color-magnitude diagram into cells
(`Hess diagram <https://en.wikipedia.org/wiki/Hess_diagram>`__); one of:
``knuth, fixed, bayes_blocks, manual``. If ``manual``
is selected, a list containing an array of edge values for the magnitude,
followed by one or two arrays (depending on the number of colors defined) for
the color(s), also with edge values, defaults to ``knuth``
``knuth, fixed, bayes_blocks``, ``fixed`` uses (15, 10) bins in magnitude and
color(s) respectively; defaults to ``knuth``
:type bin_method: str
:raises ValueError: If any of the attributes is not recognized as a valid option
Expand All @@ -37,13 +35,13 @@ def __init__(
self.lkl_name = lkl_name
self.bin_method = bin_method

likelihoods = ("plr", "bins_distance")
likelihoods = ("plr", "bins_distance", "chisq")
if self.lkl_name not in likelihoods:
raise ValueError(
f"'{self.lkl_name}' not recognized. Should be one of {likelihoods}"
)

bin_methods = ("knuth", "fixed", "bayes_blocks", "manual")
bin_methods = ("knuth", "fixed", "bayes_blocks")
if self.bin_method not in bin_methods:
raise ValueError(
f"Binning '{self.bin_method}' not recognized. "
Expand Down Expand Up @@ -80,3 +78,5 @@ def get(self, synth_clust: np.array) -> float:
# return lpriv.mean_dist(self, synth_clust)
if self.lkl_name == "bins_distance":
return lpriv.bins_distance(self, synth_clust)
if self.lkl_name == "chisq":
return lpriv.chi_square(self, synth_clust)
Loading

0 comments on commit 70c0781

Please sign in to comment.