Skip to content

Catching bootstrap_fix up with mm_rework #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mwhosek merged 8 commits into from
Jan 9, 2026
Merged

Catching bootstrap_fix up with mm_rework #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
82585b4
attempt to fix plot_stars_nfilt
Nov 18, 2025
69af17b
update startable_subset() with bootstrap errors
skterry Nov 22, 2025
2b0bccb
Remove unecessary plx functions and add motion model documentation
Dec 4, 2025
a5c988a
Merge branch 'mm_rework' of https://github.com/MovingUniverseLab/flys…
Dec 4, 2025
9aa4da7
missing tab fix
Dec 9, 2025
6de890b
test debugging
hustonm Dec 16, 2025
0663092
minor debugs and new helper functions
Jan 9, 2026
3e07730
Merge branch 'mm_rework' of https://github.com/MovingUniverseLab/flys…
Jan 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions flystar/align.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1288,6 +1288,8 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
ye=ye_trans_arr[:,boot_idx],
me=me_trans_arr[:,boot_idx],
t=np.tile(t_boot, (len(ref_table),1)))
if 'motion_model_used' in ref_table.columns:
star_table['motion_model_input'] = ref_table['motion_model_used']

# Now, do proper motion calculation, making sure to fix t0 to the
# orig value (so we can get a reasonable error on x0, y0)
Expand Down
7 changes: 4 additions & 3 deletions flystar/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -438,16 +438,17 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False):
combined astrometry + uncombined photometry table.
"""
# Multiples: ['x', 'y', 'm', 'name_in_list', 'xe', 'ye', 'me', 't',
# 'x_orig', 'y_orig', 'm_orig', 'xe_orig', 'ye_orig', 'me_orig', 'used_in_trans']
# 'x_orig', 'y_orig', 'm_orig', 'xe_orig', 'ye_orig', 'me_orig', 'used_in_trans',
# 'xe_boot','ye_boot','me_boot']
# Single: ['name', 'm0', 'm0_err', 'use_in_trans', 'ref_orig', 'n_detect',
# 'x0', 'vx', 'y0', 'vy', 'x0_err', 'vx_err', 'y0_err', 'vy_err', 't0']
# Don't include n_vfit

new_tab = copy.deepcopy(tab)
#new_tab.remove_column('n_fit')
new_tab.remove_column('n_detect')
for col in ['x','y','m','xe','ye','me','t','x_orig','y_orig','m_orig',
'xe_orig','ye_orig','me_orig','used_in_trans']:
for col in ['x','y','m','name_in_list','xe','ye','me','t','x_orig','y_orig','m_orig',
'xe_orig','ye_orig','me_orig','used_in_trans','xe_boot','ye_boot','me_boot']:
new_tab[col] = tab[col][:,idx]

new_tab.combine_lists('m', weights_col='me', sigma=3, ismag=True)
Expand Down
77 changes: 52 additions & 25 deletions flystar/motion_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,27 +26,47 @@ class MotionModel(ABC):
optional_param_names = []

def __init__(self, *args, **kwargs):
# TODO: do we need this?
'''for param in self.fitter_param_names:
param_var = getattr(self, param)
if not isinstance(param_var, (list, np.ndarray)):
setattr(self, param, np.array([param_var]))'''
"""
Make a motion model object. This object defines the fitter and fixed parameters,
and if needed stores metadata such as RA and Dec for Parallax,
for the given motion model and contains functions to fit these values to data
and apply the values to compute expected positions at given times. Each instance
corresponds to a given motion model, not an individual star, and thus the fit
values are only input/returned in functions and not stored in the object.
"""
return

def get_pos_at_time(self, params, t):
"""
Position calculator for a single star using a given motion model and input
model parameters and times.
"""
#return x, y
pass

def get_batch_pos_at_time(self, t):
"""
Position calculator for a set of stars using a given motion model and input
model parameters and times.
"""
#return x, y, x_err, y_err
pass

def run_fit(self, t, x, y, xe, ye, t0, weighting='var',
use_scipy=True, absolute_sigma=True):
"""
Run a single fit of the data to the motion model and return the best parameters.
This function is used by the overall fit_motion_model function once for a basic fit
or several times for a bootstrap fit.
"""
# Run a single fit (used both for overall fit + bootstrap iterations)
pass

def get_weights(self, xe, ye, weighting='var'):
"""
Get the weights for each data point for fitting. Options are 'var' (default)
and 'std'.
"""
if weighting=='std':
return 1./xe, 1./ye
elif weighting=='var':
Expand All @@ -56,6 +76,9 @@ def get_weights(self, xe, ye, weighting='var'):
return 1./xe**2, 1./ye**2

def scale_errors(self, errs, weighting='var'):
"""
Rescale the fit result errors as needed, according to the weighting scheme used.
"""
if weighting=='std':
return np.array(errs)**2
elif weighting=='var':
Expand All @@ -70,6 +93,7 @@ def fit_motion_model(self, t, x, y, xe, ye, t0, bootstrap=0, weighting='var',
Fit the input positions on the sky and errors
to determine new parameters for this motion model (MM).
Best-fit parameters will be returned along with uncertainties.
Optionally, bootstrap error estimation can be performed.
"""
params, param_errs = self.run_fit(t, x, y, xe, ye, t0=t0, weighting=weighting,
use_scipy=use_scipy, absolute_sigma=absolute_sigma)
Expand Down Expand Up @@ -99,8 +123,7 @@ def fit_motion_model(self, t, x, y, xe, ye, t0, bootstrap=0, weighting='var',

def get_chi2(self, fit_params, fixed_params, t, x, y, xe, ye, reduced=False):
"""
Get the chi^2 value for the current MM and
the input data.
Get the chi^2 value for the input motion model parameters and data.
"""
x_pred, y_pred = self.get_pos_at_time(fit_params, fixed_params, t)
chi2x = np.sum((x-x_pred)**2 / xe**2)
Expand Down Expand Up @@ -345,8 +368,7 @@ class Parallax(MotionModel):
"""
Motion model for linear proper motion + parallax

Requires RA, Dec, and PA parameters (degrees) for parallax calculation.
RA, Dec in J2000
Requires RA & Dec (J2000) for parallax calculation.
Optional PA is counterclockwise offset of the image y-axis from North.
Optional obs parameter describes observer location, default is 'earth'.
"""
Expand Down Expand Up @@ -450,11 +472,13 @@ def fit_func(use_t, x0,vx, y0,vy, pi):
param_errors = [x0_err, vx_err, y0_err, vy_err, pi_err]
return params, param_errors

"""
Check that everything is set up properly for motion models to run and their
required metadata.
"""

def validate_motion_model_dict(motion_model_dict, startable, default_motion_model):
"""
Check that everything is set up properly for motion models to run and their
required metadata.
"""

# Collect names of all motion models that might get used.
all_motion_model_names = ['Fixed']
if default_motion_model is not None:
Expand All @@ -478,11 +502,12 @@ def validate_motion_model_dict(motion_model_dict, startable, default_motion_mode

return motion_model_dict

"""
Get all the motion model parameters for a given motion_model_name.
Optionally, include fixed and error parameters (included by default).
"""

def get_one_motion_model_param_names(motion_model_name, with_errors=True, with_fixed=True):
"""
Get all the motion model parameters for a given motion_model_name.
Optionally, include fixed and error parameters (included by default).
"""
mod = eval(motion_model_name)
list_of_parameters = []
list_of_parameters += getattr(mod, 'fitter_param_names')
Expand All @@ -492,11 +517,12 @@ def get_one_motion_model_param_names(motion_model_name, with_errors=True, with_f
list_of_parameters += [par+'_err' for par in getattr(mod, 'fitter_param_names')]
return list_of_parameters

"""
Get all the motion model parameters for all models given in motion_model_list.
Optionally, include fixed and error parameters (included by default).
"""

def get_list_motion_model_param_names(motion_model_list, with_errors=True, with_fixed=True):
"""
Get all the motion model parameters for all models given in motion_model_list.
Optionally, include fixed and error parameters (included by default).
"""
list_of_parameters = []
all_motion_models = [eval(mm) for mm in np.unique(motion_model_list).tolist()]
for aa in range(len(all_motion_models)):
Expand All @@ -512,11 +538,12 @@ def get_list_motion_model_param_names(motion_model_list, with_errors=True, with_

return np.unique(list_of_parameters).tolist()

"""
Get all the motion model parameters for all models defined in this module.
Optionally, include fixed and error parameters (included by default).
"""

def get_all_motion_model_param_names(with_errors=True, with_fixed=True):
"""
Get all the motion model parameters for all models defined in this module.
Optionally, include fixed and error parameters (included by default).
"""
list_of_parameters = []
all_motion_models = MotionModel.__subclasses__()
for aa in range(len(all_motion_models)):
Expand Down
182 changes: 0 additions & 182 deletions flystar/parallax.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,39 +64,6 @@ def parallax_in_direction(RA, Dec, mjd, obsLocation='earth', PA=0):
return pvec


def dparallax_dt_in_direction(RA, Dec, mjd, obsLocation='earth'):
"""
R.A. in degrees. (J2000)
Dec. in degrees. (J2000)
MJD

Equations following MulensModel.
Time derivative --> units are yr^-1

"""
# print('parallax_in_direction: len(t) = ', len(mjd))
# Munge inputs into astropy format.
times = Time(mjd + 2400000.5, format='jd', scale='tdb')
coord = SkyCoord(RA, Dec, unit=(units.deg, units.deg))

direction = coord.cartesian.xyz.value
north = np.array([0., 0., 1.])
_east_projected = np.cross(north, direction) / np.linalg.norm(np.cross(north, direction))
_north_projected = np.cross(direction, _east_projected) / np.linalg.norm(np.cross(direction, _east_projected))

obs_posvel = get_observer_barycentric(obsLocation, times, velocity=True)[1]
sun_posvel = get_body_barycentric_posvel('Sun', times)[1]
sun_obs_vel = sun_posvel - obs_posvel
vel = sun_obs_vel.xyz.T.to(units.au / units.year)

e = np.dot(vel, _east_projected)
n = np.dot(vel, _north_projected)

dpvec_dt = np.array([e.value, n.value]).T

return dpvec_dt


def get_observer_barycentric(body, times, min_ephem_step=1, velocity=False):
"""
Get the barycentric position of a satellite or other Solar System body
Expand Down Expand Up @@ -180,152 +147,3 @@ def get_observer_barycentric(body, times, min_ephem_step=1, velocity=False):
return obs_pos


def sun_position(mjd, radians=False):
"""

NAME:
SUNPOS

PURPOSE:
To compute the RA and Dec of the Sun at a given date.

INPUTS:
mjd - The modified Julian date of the day (and time), scalar or vector

OUTPUTS:
ra:
| The right ascension of the sun at that date in DEGREES
| double precision, same number of elements as jd
dec:
The declination of the sun at that date in DEGREES
elong:
Ecliptic longitude of the sun at that date in DEGREES.
obliquity:
the obliquity of the ecliptic, in DEGREES

OPTIONAL INPUT KEYWORD:
RADIAN [def=False] - If this keyword is set to True, then all output variables
are given in Radians rather than Degrees

NOTES:
Patrick Wallace (Rutherford Appleton Laboratory, UK) has tested the
accuracy of a C adaptation of the sunpos.pro code and found the
following results. From 1900-2100 SUNPOS gave 7.3 arcsec maximum
error, 2.6 arcsec RMS. Over the shorter interval 1950-2050 the figures
were 6.4 arcsec max, 2.2 arcsec RMS.

The returned RA and Dec are in the given date's equinox.

Procedure was extensively revised in May 1996, and the new calling
sequence is incompatible with the old one.
METHOD:
Uses a truncated version of Newcomb's Sun. Adapted from the IDL
routine SUN_POS by CD Pike, which was adapted from a FORTRAN routine
by B. Emerson (RGO).
EXAMPLE:
(1) Find the apparent RA and Dec of the Sun on May 1, 1982

| IDL> jdcnv, 1982, 5, 1,0 ,jd ;Find Julian date jd = 2445090.5
| IDL> sunpos, jd, ra, dec
| IDL> print,adstring(ra,dec,2)
| 02 31 32.61 +14 54 34.9

The Astronomical Almanac gives 02 31 32.58 +14 54 34.9 so the error
in SUNPOS for this case is < 0.5".

(2) Find the apparent RA and Dec of the Sun for every day in 1997

| IDL> jdcnv, 1997,1,1,0, jd ;Julian date on Jan 1, 1997
| IDL> sunpos, jd+ dindgen(365), ra, dec ;RA and Dec for each day

MODIFICATION HISTORY:

* Written by Michael R. Greason, STX, 28 October 1988.
* Accept vector arguments, W. Landsman - April,1989
* Eliminated negative right ascensions - MRG, Hughes STX, 6 May 1992.
* Rewritten using the 1993 Almanac. Keywords added. MRG, HSTX, 10 February 1994.
* Major rewrite, improved accuracy, always return values in degrees - W. Landsman May, 1996
* Added /RADIAN keyword; W. Landsman; August, 1997
* Converted to IDL V5.0; W. Landsman; September 1997
* Converted to python; J. R. Lu; August 2016
"""
# form time in Julian centuries from 1900.0
t_obj = Time(mjd, format='mjd')
t = (t_obj.jd - 2415020.0) / 36525.0

# form sun's mean longitude
l = (279.696678 + ((36000.768925 * t) % 360.0)) * 3600.0

# allow for ellipticity of the orbit (equation of centre)
# using the Earth's mean anomaly ME
me = 358.475844 + ((35999.049750 * t) % 360.0)
ellcor = (6910.1 - 17.2 * t) * np.sin(np.radians(me)) + 72.3 * np.sin(
np.radians(2.0 * me))
l = l + ellcor

# allow for the Venus perturbations using the mean anomaly of Venus MV
mv = 212.603219 + ((58517.803875 * t) % 360.0)
vencorr = 4.8 * np.cos(np.radians(299.1017 + mv - me)) + \
5.5 * np.cos(np.radians(148.3133 + 2.0 * mv - 2.0 * me)) + \
2.5 * np.cos(np.radians(315.9433 + 2.0 * mv - 3.0 * me)) + \
1.6 * np.cos(np.radians(345.2533 + 3.0 * mv - 4.0 * me)) + \
1.0 * np.cos(np.radians(318.1500 + 3.0 * mv - 5.0 * me))
l += vencorr

# Allow for the Mars perturbations using the mean anomaly of Mars MM
mm = 319.529425 + ((19139.858500 * t) % 360.0)
marscorr = 2.0 * np.cos(np.radians(343.8883 - 2.0 * mm + 2.0 * me)) + \
1.8 * np.cos(np.radians(200.4017 - 2.0 * mm + me))
l += marscorr

# Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ
mj = 225.328328 + ((3034.6920239 * t) % 360.0)
jupcorr = 7.2 * np.cos(np.radians(179.5317 - mj + me)) + \
2.6 * np.cos(np.radians(263.2167 - mj)) + \
2.7 * np.cos(np.radians(87.1450 - 2.0 * mj + 2.0 * me)) + \
1.6 * np.cos(np.radians(109.4933 - 2.0 * mj + me))
l += jupcorr

# Allow for the Moons perturbations using the mean elongation of
# the Moon from the Sun D
d = 350.7376814 + ((445267.11422 * t) % 360.0)
mooncorr = 6.5 * np.sin(np.radians(d))
l += mooncorr

# Allow for long period terms
longterm = + 6.4 * np.sin(np.radians(231.19 + 20.20 * t))
l += longterm
l = (l + 2592000.0) % 1296000.0
longmed = l / 3600.0

# Allow for Aberration
l -= 20.5

# Allow for Nutation using the longitude of the Moons mean node OMEGA
omega = 259.183275 - ((1934.142008 * t) % 360.0)
l -= 17.2 * np.sin(np.radians(omega))

# Form the True Obliquity
oblt = 23.452294 - 0.0130125 * t + (
9.2 * np.cos(np.radians(omega))) / 3600.0

# Form Right Ascension and Declination
l = l / 3600.0
l_rad = np.radians(l)
oblt_rad = np.radians(oblt)
ra = np.arctan2(np.sin(l_rad) * np.cos(oblt_rad), np.cos(l_rad))

if (len(ra) > 1):
neg = np.where(ra < 0.0)[0]
ra[neg] = ra[neg] + 2.0 * math.pi

dec = np.arcsin(np.sin(l_rad) * np.sin(oblt_rad))

if radians:
oblt = oblt_rad
longmed = np.radians(longmed)
else:
ra = np.degrees(ra)
dec = np.degrees(dec)

return ra, dec, longmed, oblt
Loading