Skip to content

Commit

Permalink
Junk the old, slow, wrong broeg_weights
Browse files Browse the repository at this point in the history
  • Loading branch information
@mwcraig
mwcraig committed Oct 15, 2023
1 parent 7da2beb commit c0d4e0c
Showing 1 changed file with 36 additions and 148 deletions.
184 changes: 36 additions & 148 deletions stellarphot/differential_photometry/broeg_weights.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,161 +2,35 @@
import numpy as np


def broeg_weights(photometry, filter=None, max_iterations=100, convergence_threshold=0.001):
def broeg_weights2(mags,
errors,
max_iterations=10,
convergence_threshold=0.0001,
error_factor=1.0,
sigma_nI=0.0,):
"""
Calculate the Broeg weights for a given set of instrumental magnitudes and errors.
Parameters
----------
photometry : `astropy.table.Table`
The photometry table.
mags : array-like
The instrumental magnitudes.
filter : str, optional
The filter to use. If not specified, and only one filter is present, that
filter is used. If more than one filter is present an error is raised.
errors : array-like
The instrumental magnitude errors.
max_iterations : int, optional
The maximum number of iterations to perform.
convergence_threshold : float, optional
The convergence threshold for the iterative algorithm.
Returns
-------
array-like
The Broeg weights.
Notes
-----
This implementation of the algorithm described in C. Broeg's
algorithm ('A new algorithm for differential photometry:
computing an optimum artificial comparison
star', 2005, http://adsabs.harvard.edu/abs/2005AN....326..134B) is based
heavily on the LEMON package at https://github.com/vterron/lemon/tree/master
Specifically, see line 474 of diffphot.py in that package.
"""
# Rather than following the LEMON approach of calculating intial weights from
# instrumental magnitudes, we use the instrumental magnitude errors directly,
# as Broeg suggests in the paper.

# What filters are present in the data?
filters = set(photometry['filter'])

# If a filter is specified, check that it is present in the data
if filter is not None:
if filter not in filters:
raise ValueError('Filter {} not present in photometry table'.format(filter))

else:
# If no filter is specified, and only one filter is present, use that filter
if len(filters) == 1:
filter = filters.pop()
else:
raise ValueError('More than one filter present in photometry table. You must specify a filter.')

use_phot = photometry[photometry['filter'] == filter]

# Get number of stars in the photometry table. Do this after selecting the desired filter
# in case there are observations in filters other than the one the user wants to use.
star_ids = sorted(set(use_phot['star_id']))
n_stars = len(star_ids)

if n_stars < 3:
raise ValueError('Photometry table must contain at least three stars')

# Estimate the initial error by aggregating the error for each star over all
# of the measurements of the error.
use_phot_grouped = use_phot.group_by('star_id')

avg_phot_grouped = use_phot_grouped.groups.aggregate(np.nanmean)

# import pdb; pdb.set_trace()

# Calculate the initial Broeg weights from the instrumental errors
weights = 1.0 / avg_phot_grouped['mag_error']**2

# Weights need to be normalized
weights /= np.sum(weights)

weights_list = [weights]
light_curves = []

# This will be slow, but easy to implement for now
for i in range(max_iterations):
weights = np.zeros_like(weights)
for idx, a_star in enumerate(star_ids):
# Calculate the weighted mean magnitude of the comparison stars, excluding
# the current star
weight_table = Table([star_ids, weights_list[-1]], names=['star_id', 'weight'])
comp_mag = calc_comp_mag(use_phot, weight_table, exclude_star_id=a_star)

this_star = use_phot['star_id'] == a_star
# Calculate the difference between the two
delta_mag = use_phot[this_star]['mag_inst'] - comp_mag

# Calculate the new weight for the current star
if np.isnan(delta_mag).all() or np.nanstd(delta_mag) == 0:
print(f"Setting weight to 0 for star {a_star}")
weights[idx] = 0.0
else:
weights[idx] = 1.0 / np.nanstd(delta_mag)**2
weights_list.append(weights / np.sum(weights))
return star_ids, weights_list


def calc_comp_mag(photometry, weights_table, exclude_star_id=None):
"""
Calculate the weighted mean magnitude of the comparison stars.
The threshold of difference in weights at which to stop iterating.
Parameters
----------
error_factor : float, optional
The factor by which to multiply the errors to get the initial weights.
photometry : `astropy.table.Table`
The instrumental magnitudes of the comparison stars.
weights : array-like
The Broeg weights.
exclude_star_id : str, optional
The ID of a star to exclude from the calculation.
Returns
-------
float
The weighted mean magnitude of the comparison stars.
"""
# If a star is to be excluded, set its weight to zero
if exclude_star_id is not None:
weights_table = weights_table.copy()
# Set the weight of the excluded star to zero...
excluded_star_index = weights_table['star_id'] == exclude_star_id
weights_table['weight'][excluded_star_index] = 0.0
# ...and renormalize the weights
weights_table['weight'] /= np.sum(weights_table['weight'])

joined = join(photometry, weights_table, keys='star_id', join_type='left')
joined['weighted_mag'] = joined['mag_inst'] * joined['weight']
joined = joined['BJD', 'weighted_mag']
agged = joined.group_by('BJD').groups.aggregate(np.nansum)
# Calculate the weighted mean magnitude
return agged['weighted_mag']


def broeg_weights2(mags, errors, max_iterations=10, convergence_threshold=0.0001):
"""
Calculate the Broeg weights for a given set of instrumental magnitudes and errors.
Parameters
----------
mags : array-like
The instrumental magnitudes.
errors : array-like
The instrumental magnitude errors.
sigma_nI : float, optional
The additive term to add to the errors to get the initial weights.
Returns
-------
Expand All @@ -166,13 +40,12 @@ def broeg_weights2(mags, errors, max_iterations=10, convergence_threshold=0.0001
Notes
-----
This implementation of the algorithm described in C. Broeg's
This is an implementation of the algorithm described in C. Broeg at al's
algorithm ('A new algorithm for differential photometry:
computing an optimum artificial comparison
star', 2005, http://adsabs.harvard.edu/abs/2005AN....326..134B) is based
heavily on the LEMON package at
star', 2005, http://adsabs.harvard.edu/abs/2005AN....326..134B).
"""
weight_vec = 1.0 / (errors.mean(axis=1)**2)
weight_vec = 1.0 / (error_factor * errors.mean(axis=1)**2 + sigma_nI**2)
n_stars = mags.shape[0]

weight_vecs = [weight_vec / weight_vec.sum()]
Expand All @@ -184,22 +57,37 @@ def broeg_weights2(mags, errors, max_iterations=10, convergence_threshold=0.0001
weight_vecs = weight_vecs[:-1]
break

# Make weights into an array
weights = np.array(
[
weight_vec for _ in range(n_stars)
]
)

# Diagonal entries are zero because a star cannot be compared
# to itself.
weights[np.diag_indices(n_stars)] = 0.0

# Normalize the weights -- note the normalization factor
# is different for each star.
weights_norm_factor = 1 / weights.sum(axis=1)

# Create a diagonal matrix of the normalization factors
weights_norm_matrix = np.zeros([n_stars, n_stars])
weights_norm_matrix[np.diag_indices(n_stars)] = weights_norm_factor

# With the previous definitions, the comparison magnitudes, i.e.
# the magnitudes of the artificial comparison star, are given by
# the following matrix multiplication:
comp_mags = weights_norm_matrix @ weights @ mags

# Find the differential magnitude for each star...
delta_mags = mags - comp_mags
# ...and save it.
light_curves.append(delta_mags)
# print(f"Iteration {i}: {weight_vec}")
# print(f"\tNormalized weights: {weights_norm_matrix @ weights}")
# print(f"Delta mags: {delta_mags}")

# Calculate the new weights, which are the inverse of the
# variance of the differential magnitudes.
weight_vec = 1 / delta_mags.std(axis=1)**2
weight_vecs.append(weight_vec / weight_vec.sum())

Expand Down

0 comments on commit c0d4e0c

Please sign in to comment.