From 90461a6a84d5157ffc2620e56688da2dc53ca95e Mon Sep 17 00:00:00 2001 From: Lingfeng Wei Date: Fri, 20 Feb 2026 00:19:50 -0800 Subject: [PATCH 01/12] Fixed error on data type and mag_lim dimensions --- flystar/align.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/flystar/align.py b/flystar/align.py index db37954..4de4c37 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -255,7 +255,7 @@ def fix_iterable_conditions(self): if self.mag_lim is None: self.mag_lim = np.repeat([[None, None]], len(self.star_lists), axis=0) - elif (len(self.mag_lim) == 2): + elif (len(self.mag_lim) == 2) and (np.ndim(self.mag_lim) == 1): self.mag_lim = np.repeat([self.mag_lim], len(self.star_lists), axis=0) assert len(self.mag_lim) == len(self.star_lists) @@ -2894,8 +2894,8 @@ def trans_initial_guess(ref_list, star_list, trans_args, motion_model_dict, mode if mode == 'name': # First trim the two lists down to only those that don't contain # the "ignore_contains" string. - idx_r = np.flatnonzero(np.char.find(ref_list['name'], ignore_contains) == -1) - idx_s = np.flatnonzero(np.char.find(star_list['name'], ignore_contains) == -1) + idx_r = np.flatnonzero(np.char.find(ref_list['name'].astype(str), ignore_contains) == -1) + idx_s = np.flatnonzero(np.char.find(star_list['name'].astype(str), ignore_contains) == -1) # Match the star names name_matches, ndx_r, ndx_s = np.intersect1d(ref_list['name'][idx_r], From 3427370df838d7025a8396b057c24ede22311167 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Sat, 1 Aug 2026 22:54:38 -0700 Subject: [PATCH 02/12] Update initial guess verbose printing --- flystar/align.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/flystar/align.py b/flystar/align.py index 4de4c37..36c3878 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -2954,7 +2954,11 @@ def trans_initial_guess(ref_list, star_list, trans_args, motion_model_dict, mode trans.mag_offset = 0 if verbose > 1: - print('init guess: ', trans.px.parameters, trans.py.parameters) + # print('init guess: ', trans.px.parameters, trans.py.parameters) + print('Initial guess:') + print(f'{trans.px.parameters=}') + print(f'{trans.py.parameters=}') + print(f'{trans.mag_offset=}') warnings.filterwarnings('default', category=AstropyUserWarning) From 1d33eaf1dcf97bb92668eac0544272720264a5cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Sat, 1 Aug 2026 23:58:39 -0700 Subject: [PATCH 03/12] Added script for comparing branches --- flystar/tests/compare_branches.py | 53 +++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 flystar/tests/compare_branches.py diff --git a/flystar/tests/compare_branches.py b/flystar/tests/compare_branches.py new file mode 100644 index 0000000..4908449 --- /dev/null +++ b/flystar/tests/compare_branches.py @@ -0,0 +1,53 @@ +import pickle +import flystar +import matplotlib.pyplot as plt +from flystar import align, transforms, motion_model +from flystar.plots import plot_stars + +branch = 'mm_rework' # 'mm_compare' or 'mm_rework' + +test_data_path = f'{flystar.__path__[0]}/tests/test_data' + +with open(f'{test_data_path}/my_gaia.pkl', 'rb') as f: + my_gaia = pickle.load(f) +with open(f'{test_data_path}/list_of_starlists.pkl', 'rb') as f: + list_of_starlists = pickle.load(f) +ra_deg, dec_deg = 18.0, -30.0 +my_gaia.remove_column('motion_model_used') +if branch == 'mm_compare': + msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=3, + dr_tol=[0.2, 0.1, 0.08], dm_tol=[5,5,5], + outlier_tol=[None, None, 3], mag_lim=[6, 20], + trans_class=transforms.PolyTransform, + trans_args=[{'order': 1}, {'order': 1}, {'order': 1}], + motion_models=['Fixed'], + fixed_params_dict = {'ra':ra_deg, 'dec':dec_deg, 'pa':0.0, 'obsLocation':'earth'}, + use_ref_new=True, + update_ref_orig=False, + mag_trans=True, + trans_weights='both,std', + init_guess_mode='name', verbose=3) +elif branch == 'mm_rework': + my_gaia['motion_model_input'] = ['Fixed'] * len(my_gaia) + msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=3, + dr_tol=[0.2, 0.1, 0.08], dm_tol=[5,5,5], + outlier_tol=[None, None, 3], mag_lim=[6, 20], + trans_class=transforms.PolyTransform, + trans_args=[{'order': 1}, {'order': 1}, {'order': 1}], + default_motion_model='Fixed', + # motion_model_dict = {'Parallax': motion_model.Parallax(RA=ra_deg, Dec=dec_deg, PA=0.0, obsLocation='earth')}, + use_ref_new=True, + update_ref_orig=False, + mag_trans=True, + trans_weights='both,std', + init_guess_mode='name', verbose=3) + +msc.fit() + +with open(f'{test_data_path}/ref_table_old.pkl', 'wb') as f: + pickle.dump(msc.ref_table, f) + +for i in range(msc.ref_table['x'].shape[1]): + plt.scatter(msc.ref_table['x'][:, i], msc.ref_table['y'][:, i]) +plt.show() +plot_stars(msc.ref_table, msc.ref_table['name'][:3]) \ No newline at end of file From 05156696a0fcc3665a48f33c4f1034450df865be Mon Sep 17 00:00:00 2001 From: Lingfeng Wei Date: Sun, 2 Aug 2026 01:19:23 -0700 Subject: [PATCH 04/12] Parallax testing --- flystar/align.py | 8 +++----- flystar/tests/compare_branches.py | 20 ++++++++++---------- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/flystar/align.py b/flystar/align.py index 36c3878..3614ded 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -506,11 +506,11 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar dy=(star_t['y'] - star_r['y']) * 1e3, dm=(star_t['m'] - star_r['m']), xo=star_s['x'], yo=star_s['y'], mo=star_s['m'])) - + idx_lis, idx_ref, dr, dm = match.match(star_list_T['x'], star_list_T['y'], star_list_T['m'], ref_list['x'], ref_list['y'], ref_list['m'], dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose) - + if self.verbose > 1: print( ' Match 2: After trans, found ', len(idx_lis), ' matches out of ', len(star_list_T), '. If match count is low, check dr_tol, dm_tol.' ) @@ -858,7 +858,7 @@ def update_ref_table_aggregates(self, keep_orig=None, n_boot=0): fit_star_idxs = [idx for idx in range(len(self.ref_table)) if idx not in keep_orig] else: fit_star_idxs = None - #pdb.set_trace() + # Figure out whether motion fits are necessary all_fixed = np.all(self.ref_table['motion_model_input']=='Fixed') if all_fixed: @@ -1240,7 +1240,6 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot m=starlist_boot['m'], mref=ref_boot['m'], weights=weight, mag_trans=self.mag_trans) #print(jj) - #pdb.set_trace() # Apply transformation to *all* orig positions in this epoch. Need to make a new # FLYSTAR starlist object with the original positions for this. We don't @@ -1378,7 +1377,6 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot col[idx_good] = data_dict[ff] self.ref_table.add_column(col) - #pdb.set_trace() print('===============================') print('Done with bootstrap') diff --git a/flystar/tests/compare_branches.py b/flystar/tests/compare_branches.py index 4908449..f764309 100644 --- a/flystar/tests/compare_branches.py +++ b/flystar/tests/compare_branches.py @@ -14,13 +14,14 @@ list_of_starlists = pickle.load(f) ra_deg, dec_deg = 18.0, -30.0 my_gaia.remove_column('motion_model_used') +# my_gaia['motion_model_input'] = 'Fixed' if branch == 'mm_compare': msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=3, dr_tol=[0.2, 0.1, 0.08], dm_tol=[5,5,5], outlier_tol=[None, None, 3], mag_lim=[6, 20], trans_class=transforms.PolyTransform, trans_args=[{'order': 1}, {'order': 1}, {'order': 1}], - motion_models=['Fixed'], + motion_models=['Fixed', 'Parallax'], fixed_params_dict = {'ra':ra_deg, 'dec':dec_deg, 'pa':0.0, 'obsLocation':'earth'}, use_ref_new=True, update_ref_orig=False, @@ -28,14 +29,13 @@ trans_weights='both,std', init_guess_mode='name', verbose=3) elif branch == 'mm_rework': - my_gaia['motion_model_input'] = ['Fixed'] * len(my_gaia) msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=3, dr_tol=[0.2, 0.1, 0.08], dm_tol=[5,5,5], outlier_tol=[None, None, 3], mag_lim=[6, 20], trans_class=transforms.PolyTransform, trans_args=[{'order': 1}, {'order': 1}, {'order': 1}], - default_motion_model='Fixed', - # motion_model_dict = {'Parallax': motion_model.Parallax(RA=ra_deg, Dec=dec_deg, PA=0.0, obsLocation='earth')}, + default_motion_model='Parallax', + motion_model_dict = {'Parallax': motion_model.Parallax(RA=ra_deg, Dec=dec_deg, PA=0.0, obsLocation='earth')}, use_ref_new=True, update_ref_orig=False, mag_trans=True, @@ -44,10 +44,10 @@ msc.fit() -with open(f'{test_data_path}/ref_table_old.pkl', 'wb') as f: - pickle.dump(msc.ref_table, f) +# with open(f'{test_data_path}/ref_table_old.pkl', 'wb') as f: +# pickle.dump(msc.ref_table, f) -for i in range(msc.ref_table['x'].shape[1]): - plt.scatter(msc.ref_table['x'][:, i], msc.ref_table['y'][:, i]) -plt.show() -plot_stars(msc.ref_table, msc.ref_table['name'][:3]) \ No newline at end of file +# for i in range(msc.ref_table['x'].shape[1]): +# plt.scatter(msc.ref_table['x'][:, i], msc.ref_table['y'][:, i]) +# plt.show() +# plot_stars(msc.ref_table, msc.ref_table['name'][:3]) \ No newline at end of file From 8f65230ca57881e9e457494ecaf5205144fac988 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Mon, 3 Aug 2026 19:23:00 -0700 Subject: [PATCH 05/12] Change to relative imports --- flystar/align.py | 14 +++++++------- flystar/analysis.py | 6 +----- flystar/{tests => }/compare_branches.py | 4 ++-- flystar/examples.py | 6 +----- flystar/match.py | 2 +- flystar/motion_model.py | 2 +- flystar/plots.py | 2 +- flystar/startables.py | 2 +- flystar/transforms.py | 2 +- 9 files changed, 16 insertions(+), 24 deletions(-) rename flystar/{tests => }/compare_branches.py (96%) diff --git a/flystar/align.py b/flystar/align.py index 3614ded..c1295f9 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -1,10 +1,10 @@ import numpy as np -from flystar import match -from flystar import transforms -from flystar import plots -from flystar.starlists import StarList -from flystar.startables import StarTable -from flystar import motion_model +from . import match +from . import transforms +from . import plots +from .starlists import StarList +from .startables import StarTable +from . import motion_model from astropy.table import Table, Column, vstack import datetime import copy @@ -419,7 +419,7 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar star_list_T.transform_xym(trans) # trimmed, transformed else: star_list_T.transform_xy(trans) - + pdb.set_trace() # Match stars between the transformed, trimmed lists. idx1, idx2, dr, dm = match.match(star_list_T['x'], star_list_T['y'], star_list_T['m'], ref_list['x'], ref_list['y'], ref_list['m'], diff --git a/flystar/analysis.py b/flystar/analysis.py index f300723..966b913 100644 --- a/flystar/analysis.py +++ b/flystar/analysis.py @@ -1,10 +1,6 @@ import numpy as np import pylab as plt -from flystar import starlists -from flystar import startables -from flystar import align -from flystar import match -from flystar import transforms +from . import starlists, startables, align, match, transforms from astropy import table from astropy.table import Table, Column from astropy.coordinates import SkyCoord diff --git a/flystar/tests/compare_branches.py b/flystar/compare_branches.py similarity index 96% rename from flystar/tests/compare_branches.py rename to flystar/compare_branches.py index f764309..49f9731 100644 --- a/flystar/tests/compare_branches.py +++ b/flystar/compare_branches.py @@ -1,8 +1,8 @@ import pickle import flystar import matplotlib.pyplot as plt -from flystar import align, transforms, motion_model -from flystar.plots import plot_stars +import align, transforms, motion_model +from plots import plot_stars branch = 'mm_rework' # 'mm_compare' or 'mm_rework' diff --git a/flystar/examples.py b/flystar/examples.py index 8059562..0e0a042 100644 --- a/flystar/examples.py +++ b/flystar/examples.py @@ -1,8 +1,4 @@ -from flystar import transforms -from flystar import match -from flystar import align -from flystar import starlists -from flystar import plots +from . import transforms, match, align, starlists, plots import numpy as np import copy import pdb diff --git a/flystar/match.py b/flystar/match.py index bba108a..fe749f0 100644 --- a/flystar/match.py +++ b/flystar/match.py @@ -1,5 +1,5 @@ import numpy as np -from flystar import starlists, transforms, startables, align +from . import starlists, transforms, startables, align from collections import Counter from scipy.spatial import cKDTree as KDT from astropy.table import Column, Table diff --git a/flystar/motion_model.py b/flystar/motion_model.py index 0b86d07..2b9e570 100644 --- a/flystar/motion_model.py +++ b/flystar/motion_model.py @@ -1,7 +1,7 @@ import numpy as np from abc import ABC import pdb -from flystar import parallax +from . import parallax from astropy.time import Time from scipy.optimize import curve_fit import warnings diff --git a/flystar/plots.py b/flystar/plots.py index 2d65b2c..7d4b700 100755 --- a/flystar/plots.py +++ b/flystar/plots.py @@ -1,4 +1,4 @@ -from flystar import analysis, motion_model, startables +from . import analysis, motion_model, startables import pylab as py import pylab as plt import numpy as np diff --git a/flystar/startables.py b/flystar/startables.py index d75fca9..af7e6e5 100644 --- a/flystar/startables.py +++ b/flystar/startables.py @@ -9,7 +9,7 @@ import pdb import time import copy -from flystar import motion_model +from . import motion_model import pandas as pd class StarTable(Table): diff --git a/flystar/transforms.py b/flystar/transforms.py index 6cc865a..4f4410c 100755 --- a/flystar/transforms.py +++ b/flystar/transforms.py @@ -6,7 +6,7 @@ import collections import re import pdb -from flystar import motion_model +from . import motion_model class Transform2D(object): ''' From 42e80ef7ded8b30bd6e01c3e9b266e6dabcb065e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Thu, 6 Aug 2026 01:13:06 -0700 Subject: [PATCH 06/12] Match the mm_rework_lingfeng behavior --- flystar/align.py | 6 +++--- flystar/motion_model.py | 42 +++++++++++++++++++++++++++-------------- flystar/parallax.py | 2 +- flystar/startables.py | 27 ++++++++++++++++++-------- 4 files changed, 51 insertions(+), 26 deletions(-) diff --git a/flystar/align.py b/flystar/align.py index c1295f9..3914378 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -28,7 +28,7 @@ def __init__(self, list_of_starlists, ref_index=0, iters=2, default_motion_model='Fixed', motion_model_dict = {}, use_scipy=True, - absolute_sigma=False, + absolute_sigma=True, save_path=None, verbose=True): """ @@ -419,7 +419,7 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar star_list_T.transform_xym(trans) # trimmed, transformed else: star_list_T.transform_xy(trans) - pdb.set_trace() + # Match stars between the transformed, trimmed lists. idx1, idx2, dr, dm = match.match(star_list_T['x'], star_list_T['y'], star_list_T['m'], ref_list['x'], ref_list['y'], ref_list['m'], @@ -1414,7 +1414,7 @@ def __init__(self, ref_list, list_of_starlists, iters=2, default_motion_model='Fixed', motion_model_dict={}, use_scipy=True, - absolute_sigma=False, + absolute_sigma=True, save_path=None, verbose=True): diff --git a/flystar/motion_model.py b/flystar/motion_model.py index 2b9e570..3564098 100644 --- a/flystar/motion_model.py +++ b/flystar/motion_model.py @@ -61,7 +61,16 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', """ # Run a single fit (used both for overall fit + bootstrap iterations) pass - + + def calc_sigma(self, xe, ye, weighting='var'): + if weighting=='std': + return np.sqrt(np.abs(xe)), np.sqrt(np.abs(ye)) + elif weighting=='var': + return np.abs(xe), np.abs(ye) + else: + warnings.warn("Invalid weighting, using default weighting scheme var.", UserWarning) + return np.abs(xe), np.abs(ye) + def get_weights(self, xe, ye, weighting='var'): """ Get the weights for each data point for fitting. Options are 'var' (default) @@ -174,11 +183,16 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, x0,y0,x0e,y0e = x[0],y[0],xe[0],ye[0] else: - x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting) + sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting) + x_wt, y_wt = 1. / sigma_x**2, 1. / sigma_y**2 + x_wt /= np.sum(x_wt) + y_wt /= np.sum(y_wt) x0 = np.average(x, weights=x_wt) - x0e = np.sqrt(np.average((x-x0)**2,weights=x_wt)) + # x0e = np.sqrt(np.average((x-x0)**2,weights=x_wt)) + x0e = np.sum(x_wt**2 * xe**2)**0.5 # Error propagation y0 = np.average(y, weights=y_wt) - y0e = np.sqrt(np.average((y-y0)**2,weights=y_wt)) + # y0e = np.sqrt(np.average((y-y0)**2,weights=y_wt)) + y0e = np.sum(y_wt**2 * ye**2)**0.5 # Error propagation params = [x0, y0] param_errors = [x0e, y0e] @@ -226,7 +240,6 @@ def get_batch_pos_at_time(self, t, x0=[],vx=[], y0=[],vy=[], t0=[], def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, use_scipy=True, absolute_sigma=True): dt = t-t0 - x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting) if params_guess is None: params_guess = [x.mean(),0.0,y.mean(),0.0] @@ -256,8 +269,9 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, if use_scipy: def linear(t, c0, c1): return c0 + c1*t - x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=1/np.sqrt(x_wt), absolute_sigma=absolute_sigma) - y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=1/np.sqrt(y_wt), absolute_sigma=absolute_sigma) + sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting) + x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=sigma_x, absolute_sigma=absolute_sigma) + y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=sigma_y, absolute_sigma=absolute_sigma) x0, vx = x_opt y0, vy = y_opt x0e, vxe = np.sqrt(x_cov.diagonal()) @@ -339,15 +353,14 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, if not use_scipy: Warning("Acceleration model has no non-scipy fitter option. Running with scipy.") dt = t-t0 - x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting) if params_guess is None: params_guess = [x.mean(),0.0,0.0,y.mean(),0.0,0.0] def accel(t, c0,c1,c2): return c0 + c1*t + 0.5*c2*t**2 - - x_opt, x_cov = curve_fit(accel, dt, x, p0=np.array(params_guess[:3]), sigma=1/x_wt**0.5, absolute_sigma=True) - y_opt, y_cov = curve_fit(accel, dt, y, p0=np.array(params_guess[3:]), sigma=1/y_wt**0.5, absolute_sigma=True) + sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting) + x_opt, x_cov = curve_fit(accel, dt, x, p0=np.array(params_guess[:3]), sigma=sigma_x, absolute_sigma=True) + y_opt, y_cov = curve_fit(accel, dt, y, p0=np.array(params_guess[3:]), sigma=sigma_y, absolute_sigma=True) x0 = x_opt[0] y0 = y_opt[0] vx0 = x_opt[1] @@ -372,7 +385,7 @@ class Parallax(MotionModel): Optional PA is counterclockwise offset of the image y-axis from North. Optional obs parameter describes observer location, default is 'earth'. """ - n_pts_req = 4 + n_pts_req = 3 n_params=3 fitter_param_names = ['x0', 'vx', 'y0', 'vy', 'pi'] fixed_param_names = ['t0'] @@ -451,7 +464,6 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, Warning("Parallax model has no non-scipy fitter option. Running with scipy.") t_mjd = Time(t, format='decimalyear', scale='utc').mjd pvec = self.get_parallax_vector(t_mjd) - x_wt, y_wt = self.get_weights(xe,ye, weighting=weighting) def fit_func(use_t, x0,vx, y0,vy, pi): x_res = x0 + vx*(use_t-t0) + pi*pvec[0] y_res = y0 + vy*(use_t-t0) + pi*pvec[1] @@ -463,8 +475,10 @@ def fit_func(use_t, x0,vx, y0,vy, pi): idx_first, idx_last = np.argmin(t), np.argmax(t) params_guess = [x.mean(),(x[idx_last]-x[idx_first])/(t[idx_last]-t[idx_first]), y.mean(),(y[idx_last]-y[idx_first])/(t[idx_last]-t[idx_first]), 0.1] + sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting) + sigma = np.hstack([sigma_x, sigma_y]) res = curve_fit(fit_func, t, np.hstack([x,y]), - p0=params_guess, sigma = 1.0/np.hstack([x_wt,y_wt])) + p0=params_guess, sigma=sigma, absolute_sigma=absolute_sigma) x0,vx,y0,vy,pi = res[0] x0_err,vx_err,y0_err,vy_err,pi_err = self.scale_errors(np.sqrt(np.diag(res[1])), weighting=weighting) diff --git a/flystar/parallax.py b/flystar/parallax.py index 4792ec6..586bda8 100755 --- a/flystar/parallax.py +++ b/flystar/parallax.py @@ -36,7 +36,7 @@ def parallax_in_direction(RA, Dec, mjd, obsLocation='earth', PA=0): #print('parallax_in_direction: len(t) = ', len(mjd)) # Munge inputs into astropy format. - times = Time(mjd + 2400000.5, format='jd', scale='tdb') + times = Time(mjd, format='mjd', scale='tdb') coord = SkyCoord(RA, Dec, unit=(units.deg, units.deg)) direction = coord.cartesian.xyz.value diff --git a/flystar/startables.py b/flystar/startables.py index af7e6e5..185aea3 100644 --- a/flystar/startables.py +++ b/flystar/startables.py @@ -463,6 +463,11 @@ def combine_lists(self, col_name_in, weights_col=None, mask_val=None, if all(isinstance(item, int) for item in mask_lists): val_2d.mask[:, mask_lists] = True + use_lists = np.array([i for i in np.arange(self[col_name_in].data.shape[1]) if i not in mask_lists]) + else: + # Use all indices + use_lists = np.arange(self[col_name_in].data.shape[1]) + # Throw a warning if mask_lists is not a list if not isinstance(mask_lists, list): raise RuntimeError('mask_lists needs to be a list.') @@ -500,15 +505,18 @@ def combine_lists(self, col_name_in, weights_col=None, mask_val=None, # the N_lists direction (axis=1). if wgt_2d is not None: avg = np.ma.average(val_2d_clip, weights=wgt_2d, axis=1) - std = np.sqrt(np.ma.average((val_2d_clip.T - avg).T**2, weights=wgt_2d, axis=1)) + # std = np.sqrt(np.ma.average((val_2d_clip.T - avg).T**2, weights=wgt_2d, axis=1)) + std = np.ma.sqrt(1. / np.ma.sum(wgt_2d, axis=1)) # Error propagation else: avg = np.ma.mean(val_2d_clip, axis=1) - std = np.ma.std(val_2d_clip, axis=1) + # std = np.ma.std(val_2d_clip, axis=1) + std = np.ma.std(val_2d_clip, axis=1) / np.sqrt(len(use_lists)) # Error propagation # To Do: bring the previous uncertainties of stars that are detected # in only one input frame. - if (weights_col and weights_col in self.colnames) and (val_2d.shape[1] > 1): - mask_for_singles = ((~np.isnan(val_2d_clip)).sum(axis=1)==1) - std[mask_for_singles]=np.nanmean(err_2d[mask_for_singles], axis=1) + # This can be removed now as error propagation won't result in std=0 anymore. + # if (weights_col and weights_col in self.colnames) and (val_2d.shape[1] > 1): + # mask_for_singles = ((~np.isnan(val_2d_clip)).sum(axis=1)==1) + # std[mask_for_singles]=np.nanmean(err_2d[mask_for_singles], axis=1) # Save off our new AVG and STD into new columns with shape (N_stars). col_name_avg = col_name_in + '0' @@ -607,9 +615,12 @@ def fit_velocities(self, weighting='var', use_scipy=True, absolute_sigma=True, b # Define output arrays for the best-fit parameters. for col in new_col_list: # Clean/remove up old arrays. - if col in self.colnames: self.remove_column(col) - # Add column #TODO: is this good for filling??? - self.add_column(Column(data = np.full(N_stars, np.nan, dtype=float), name = col)) + # if col in self.colnames: self.remove_column(col) + # # Add column #TODO: is this good for filling??? + # self.add_column(Column(data = np.full(N_stars, np.nan, dtype=float), name=col)) + # Keep existing values + if col not in self.colnames: + self.add_column(Column(data = np.full(N_stars, np.nan, dtype=float), name=col)) # Add a column to keep track of the number of points used in a fit. self['n_fit'] = 0 From eadf9824e48344d7358169c97f152ae44e59a6c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Thu, 6 Aug 2026 02:07:52 -0700 Subject: [PATCH 07/12] Fixed remaining problems of mm_rework and matched mm_rework_lingfeng's result --- flystar/align.py | 32 ++++++++++++++++++-------------- flystar/compare_branches.py | 32 ++++++++++++++++++++++---------- flystar/startables.py | 9 +++++---- 3 files changed, 45 insertions(+), 28 deletions(-) diff --git a/flystar/align.py b/flystar/align.py index 3914378..e1bac80 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -420,10 +420,24 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar else: star_list_T.transform_xy(trans) - # Match stars between the transformed, trimmed lists. - idx1, idx2, dr, dm = match.match(star_list_T['x'], star_list_T['y'], star_list_T['m'], - ref_list['x'], ref_list['y'], ref_list['m'], - dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose) + if 'use_in_trans' in ref_list.colnames: + # Only use stars specified by "use_in_trans" column. + use_in_trans = ref_list['use_in_trans'] + # Match stars between the transformed, trimmed lists. + idx1, idx2, dr, dm = match.match( + star_list_T['x'], star_list_T['y'], star_list_T['m'], + ref_list['x'][use_in_trans], ref_list['y'][use_in_trans], ref_list['m'][use_in_trans], + dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose + ) + # Restore idx2 to the full reference list indices + idx2 = np.where(use_in_trans)[0][idx2] + else: + idx1, idx2, dr, dm = match.match( + star_list_T['x'], star_list_T['y'], star_list_T['m'], + ref_list['x'], ref_list['y'], ref_list['m'], + dr_tol=dr_tol, dm_tol=dm_tol, verbose=self.verbose + ) + if self.verbose > 1: print( ' Match 1: Found ', len(idx1), ' matches out of ', len(star_list_T), '. If match count is low, check dr_tol, dm_tol.' ) @@ -438,16 +452,6 @@ def match_and_transform(self, ref_mag_lim, dr_tol, dm_tol, outlier_tol, trans_ar idx1 = idx1[keepers] idx2 = idx2[keepers] - # Only use stars specified by "use_in_trans" column. - if 'use_in_trans' in ref_list.colnames: - keepers = np.where(ref_list[idx2]['use_in_trans'] == True)[0] - - if self.verbose > 1: - print( ' Rejected ', len(idx1) - len(keepers), ' with use_in_trans=False.' ) - - idx1 = idx1[keepers] - idx2 = idx2[keepers] - # Determine weights in the fit. weight = self.get_weights_for_lists(ref_list[idx2], star_list_T[idx1]) diff --git a/flystar/compare_branches.py b/flystar/compare_branches.py index 49f9731..d5f1e51 100644 --- a/flystar/compare_branches.py +++ b/flystar/compare_branches.py @@ -1,12 +1,24 @@ +import os +import sys +from pathlib import Path + +# 1. Calculate the absolute path to the parent directory +parent_dir = str(Path(__file__).resolve().parent) + +# 2. Temporarily inject it into Python's search path +if parent_dir not in sys.path: + sys.path.insert(0, parent_dir) + + +import os import pickle -import flystar import matplotlib.pyplot as plt import align, transforms, motion_model from plots import plot_stars branch = 'mm_rework' # 'mm_compare' or 'mm_rework' -test_data_path = f'{flystar.__path__[0]}/tests/test_data' +test_data_path = f'{os.path.expanduser("~")}/Software/flystar/flystar/tests/test_data' with open(f'{test_data_path}/my_gaia.pkl', 'rb') as f: my_gaia = pickle.load(f) @@ -16,11 +28,11 @@ my_gaia.remove_column('motion_model_used') # my_gaia['motion_model_input'] = 'Fixed' if branch == 'mm_compare': - msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=3, - dr_tol=[0.2, 0.1, 0.08], dm_tol=[5,5,5], - outlier_tol=[None, None, 3], mag_lim=[6, 20], + msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=1, + dr_tol=[0.2], dm_tol=[5], + outlier_tol=[None], mag_lim=[6, 20], trans_class=transforms.PolyTransform, - trans_args=[{'order': 1}, {'order': 1}, {'order': 1}], + trans_args=[{'order': 1}], motion_models=['Fixed', 'Parallax'], fixed_params_dict = {'ra':ra_deg, 'dec':dec_deg, 'pa':0.0, 'obsLocation':'earth'}, use_ref_new=True, @@ -29,11 +41,11 @@ trans_weights='both,std', init_guess_mode='name', verbose=3) elif branch == 'mm_rework': - msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=3, - dr_tol=[0.2, 0.1, 0.08], dm_tol=[5,5,5], - outlier_tol=[None, None, 3], mag_lim=[6, 20], + msc = align.MosaicToRef(my_gaia, list_of_starlists, iters=1, + dr_tol=[0.2], dm_tol=[5], + outlier_tol=[None], mag_lim=[6, 20], trans_class=transforms.PolyTransform, - trans_args=[{'order': 1}, {'order': 1}, {'order': 1}], + trans_args=[{'order': 1}], default_motion_model='Parallax', motion_model_dict = {'Parallax': motion_model.Parallax(RA=ra_deg, Dec=dec_deg, PA=0.0, obsLocation='earth')}, use_ref_new=True, diff --git a/flystar/startables.py b/flystar/startables.py index 185aea3..f4a04c1 100644 --- a/flystar/startables.py +++ b/flystar/startables.py @@ -464,13 +464,14 @@ def combine_lists(self, col_name_in, weights_col=None, mask_val=None, val_2d.mask[:, mask_lists] = True use_lists = np.array([i for i in np.arange(self[col_name_in].data.shape[1]) if i not in mask_lists]) - else: - # Use all indices - use_lists = np.arange(self[col_name_in].data.shape[1]) # Throw a warning if mask_lists is not a list if not isinstance(mask_lists, list): - raise RuntimeError('mask_lists needs to be a list.') + raise RuntimeError(f'mask_lists needs to be a list., not {type(mask_lists)}') + + else: + # Use all indices + use_lists = np.arange(self[col_name_in].data.shape[1]) # Decide if we are going to have weights (before we # do the expensive sigma clipping routine). Note that From 960b017b2eca9dd4affbde6e71f30faa7fa2cff2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Thu, 6 Aug 2026 15:30:27 -0700 Subject: [PATCH 08/12] Moved compare branches to tests folder --- flystar/{ => tests}/compare_branches.py | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) rename flystar/{ => tests}/compare_branches.py (86%) diff --git a/flystar/compare_branches.py b/flystar/tests/compare_branches.py similarity index 86% rename from flystar/compare_branches.py rename to flystar/tests/compare_branches.py index d5f1e51..ed55b39 100644 --- a/flystar/compare_branches.py +++ b/flystar/tests/compare_branches.py @@ -1,20 +1,8 @@ -import os -import sys -from pathlib import Path - -# 1. Calculate the absolute path to the parent directory -parent_dir = str(Path(__file__).resolve().parent) - -# 2. Temporarily inject it into Python's search path -if parent_dir not in sys.path: - sys.path.insert(0, parent_dir) - - import os import pickle import matplotlib.pyplot as plt -import align, transforms, motion_model -from plots import plot_stars +from flystar import align, transforms, motion_model +from flystar.plots import plot_stars branch = 'mm_rework' # 'mm_compare' or 'mm_rework' From 754cb7811b22a0750354a1eba822119edfc1f551 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Thu, 6 Aug 2026 16:02:29 -0700 Subject: [PATCH 09/12] Cleaned up imports to avoid cyclic imports --- flystar/align.py | 270 +++++++++++++++++++++++++++++++--------- flystar/analysis.py | 84 +------------ flystar/match.py | 208 ------------------------------- flystar/motion_model.py | 24 +--- flystar/plots.py | 111 ++++++++++++++--- flystar/transforms.py | 2 +- 6 files changed, 318 insertions(+), 381 deletions(-) diff --git a/flystar/align.py b/flystar/align.py index e1bac80..ff08ba6 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -1,10 +1,7 @@ import numpy as np -from . import match -from . import transforms -from . import plots -from .starlists import StarList -from .startables import StarTable -from . import motion_model +from flystar import match, transforms, plots, motion_model +from flystar.starlists import StarList +from flystar.startables import StarTable from astropy.table import Table, Column, vstack import datetime import copy @@ -2405,8 +2402,7 @@ def write_transform(transform, starlist, reference, N_trans, deltaMag=0, restric Xcoeff = transform.px.parameters Ycoeff = transform.py.parameters else: - print(( '{0} not yet supported!'.format(transType))) - return + raise TypeError(( '{0} not yet supported!'.format(trans_name))) # Write output _out = open(outFile, 'w') @@ -2634,7 +2630,7 @@ def position_transform_from_object(x, y, xe, ye, transform): order = transform.order else: txt = 'Transform not yet supported by position_transform_from_object' - raise StandardError(txt) + raise TypeError(txt) # How the transformation is applied depends on the type of transform. # This can be determined by the length of Xcoeff, Ycoeff @@ -2733,7 +2729,7 @@ def velocity_transform_from_object(x0, y0, x0e, y0e, vx, vy, vxe, vye, transform order = transform.order else: txt = 'Transform not yet supported by velocity_transform_from_object' - raise StandardError(txt) + raise TypeError(txt) # How the transformation is applied depends on the type of transform. # This can be determined by the length of Xcoeff, Ycoeff @@ -3024,53 +3020,6 @@ def copy_and_rename_for_ref(star_list): return ref_list -def outlier_rejection_indices(star_list, ref_list, outlier_tol, verbose=True): - """ - Determine the outliers based on the residual positions between two different - starlists and some threshold (in sigma). Return the indices of the stars - to keep (that shouldn't be rejected as outliers). - - Note that we assume that the star_list and ref_list are already transformed and - matched. - - Parameters - ---------- - star_list : StarList - starlist with 'x', 'y' - - ref_list : StarList - starlist with 'x0', 'y0' - - outlier_tol : float - Number of sigma inside which we keep stars and outside of which we - reject stars as outliers. - - Optional Parameters - -------------------- - verbose : boolean - - Returns - ---------- - keepers : nd.array - The indicies of the stars to keep. - """ - # Optionally propogate the reference positions forward in time. - xref, yref = get_pos_in_time(star_list['t'][0], ref_list) - - # Residuals - x_resid_on_old_trans = star_list['x'] - xref - y_resid_on_old_trans = star_list['y'] - yref - resid_on_old_trans = np.hypot(x_resid_on_old_trans, y_resid_on_old_trans) - - threshold = outlier_tol * resid_on_old_trans.std() - keepers = np.where(resid_on_old_trans < threshold)[0] - - if verbose > 0: - msg = ' Outlier Rejection: Keeping {0:d} of {1:d}' - print(msg.format(len(keepers), len(resid_on_old_trans))) - - return keepers - def setup_trans_info(trans_input, trans_args, N_lists, iters): """ Setup transformation info into a usable format. @@ -3199,3 +3148,210 @@ def logger(logfile, message, verbose = 9): print(message) logfile.write(message + '\n') return + + +def generic_match(sl1, sl2, init_mode='triangle', + model=transforms.PolyTransform, order_dr=(1, 1.0), + dr_final=1.0, + xy_match=(None, None, None, None, None, None, None, None), + m_match=(None, None, None, None), sigma_match=None, + n_bright=100, verbose=True, **kwargs): + """ + Finds the transformation between two starlists using the first one + as reference frame. Different matching methods can be used. If no + transformation is found, it returns an error message. + + + Parameters + sl1 : StarList + starlist used for reference frame + sl2 : StarList + starlist transformed + init_mode : str + Initial matching method. + If 'triangle', uses the blind triangle method. + If 'match_name', uses match by name + If 'load', uses the transformation from a loaded file + model : str + Transformation model to be used with the 'triangle' initial mode + poly_order : int + Order of the transformation model + order_dr : int, float [n, 2] + Combinations of polinomial order (first column) and search radius + (second column) to refine the transformation. Rows are executed in + orders + dr_final: float + Search radius used for the final matching + n_bright : int + Number of bright stars used in the initial blind triangles matching + xy_match : array + Area of the images to remove in the matching [reference catalog min x, + reference catalog max x, reference catalog min y, reference catalog max y, + transformed catalog min x, transformed catalog max x, + transformed catalog min y, transformed catalog max y]. Use None for values not used. + m_match : array + Magnitude limits of matching stars used to find transformations + [reference catalog min mag, reference catalog max mag, transformed + catalog min mag, transformed catalog max mag]. Use None for values not + used + sigma_match : array + Number of Deltap movement sigmas [0] used for sigma-cutting matched + stars for a number of times [1]. Use None for no sigma-cut. The last + polynomial order and search radius in 'order_dr' are used + transf_file : str + File name and path of the transformation file used with the 'load' + init_mode + verbose : bool, optional + Prints on screen information on the matching + + Returns + ------- + transf : Transform2D + Transformation of the second starlist respect to the first + st : StarTable + Startable of the two matched catalogs + + """ + from flystar import starlists, startables + # Check the input StarLists and transform them into astropy Tables + if not isinstance(sl1, starlists.StarList): + raise TypeError("The first catalog has to be a StarList") + if not isinstance(sl2, starlists.StarList): + raise TypeError("The second catalog has to be a StarList") + + # Find the initial transformation + if init_mode == 'triangle': # Blind triangles method + + # Prepare the reduced starlists for matching + sl1_cut = copy.deepcopy(sl1) + sl2_cut = copy.deepcopy(sl2) + sl1_cut.restrict_by_value(x_min=xy_match[0], x_max=xy_match[1], + y_min=xy_match[2], y_max=xy_match[3]) + sl2_cut.restrict_by_value(x_min=xy_match[4], x_max=xy_match[5], + y_min=xy_match[6], y_max=xy_match[7]) + sl1_cut.restrict_by_value(m_min=m_match[0], m_max=m_match[1]) + sl2_cut.restrict_by_value(m_min=m_match[2], m_max=m_match[3]) + + # Find the transformation + # TODO: test 'initial_align' with StarList input + transf = initial_align(sl1_cut, sl2_cut, briteN=n_bright, + transformModel=model, order=order_dr[0]) #order_dr[i_loop][0] ? + + elif init_mode == 'match_name': # Name match + sl1_idx_init, sl2_idx_init, _ = starlists.restrict_by_name(sl1, sl2) + transf = model(sl2['x'][sl2_idx_init], sl2['y'][sl2_idx_init], + sl1['x'][sl1_idx_init], sl1['y'][sl1_idx_init], + order=int(order_dr[0][0])) + + elif init_mode == 'load': # Load a transformation file + transf = transforms.Transform2D.from_file(kwargs['transf_file']) + + else: # None of the above + raise TypeError("Unrecognized initial matching method") + + # Restrict the matching catalogs + sl1_match = copy.deepcopy(sl1) + sl2_match = copy.deepcopy(sl2) + sl1_match.restrict_by_value(m_min=m_match[0], m_max=m_match[1]) + sl2_match.restrict_by_value(m_min=m_match[2], m_max=m_match[3]) + + # Refine the transformation + if sigma_match: + order_dr_len = len(order_dr) + + for i_loop in range(sigma_match[1]): + order_dr = np.vstack((np.array(order_dr), np.array(order_dr[-1]))) + + for i_loop in range(len(order_dr)): + + # Transform and match the catalog to the reference frame +# sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf, +# dr_tol=order_dr[i_loop][1], +# verbose=verbose) + + sl2_idx, sl1_idx = transform_and_match(sl2_match, sl1_match, transf, + dr_tol=order_dr[1], + verbose=verbose) + + # Transform the catalog to the reference frame + sl2_transf_match = transform_from_object(sl2_match, transf) + + # Sigma-rejection + if sigma_match and (i_loop >= order_dr_len): + resid = np.sqrt((sl1_match['x'][sl1_idx] - + sl2_transf_match['x'][sl2_idx])**2 + + (sl1_match['y'][sl1_idx] - + sl2_transf_match['y'][sl2_idx])**2) + sl1_idx = sl1_idx[resid <= (sigma_match[0] * np.std(resid))] + sl2_idx = sl2_idx[resid <= (sigma_match[0] * np.std(resid))] + + # Test section to observe the matching catalogs before refining the transformation + """ + from matplotlib import pyplot + + _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10,10)) + axarr.scatter(sl1_match['x'][sl1_idx], sl1_match['y'][sl1_idx]) + xlim = axarr.get_xlim() + ylim = axarr.get_ylim() + + _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10, 10)) + axarr.scatter(sl2_transf_match['x'][sl2_idx], sl2_transf_match['y'][sl2_idx]) + axarr.set_xlim(xlim) + axarr.set_ylim(ylim) + """ + + # Find a better transformation + transf, _ = find_transform(sl2_match[sl2_idx], + sl2_transf_match[sl2_idx], + sl1_match[sl1_idx], transModel=model, + order=order_dr[0], verbose=verbose) +# order=int(order_dr[i_loop][0]), verbose=verbose) + + # This section was used for testing transformations with normalized + # coordinates. Only several catalogs had reduced residuals when using + # high order polynomials (>3), some of them became unstable + """sl1_match_norm = sl1_match[sl1_idx] + sl2_match_norm = sl2_match[sl2_idx] + sl2_transf_match_norm = sl2_transf_match[sl2_idx] + mm = max(max(sl1_match_norm['x']), max(sl1_match_norm['y']), + max(sl2_transf_match_norm['x']), max(sl2_transf_match_norm['y'])) + sl1_match_norm['x'] = sl1_match_norm['x'] / mm + sl1_match_norm['y'] = sl1_match_norm['y'] / mm + sl2_match_norm['x'] = sl2_match_norm['x'] / mm + sl2_match_norm['y'] = sl2_match_norm['y'] / mm + sl2_transf_match_norm['x'] = sl2_transf_match_norm['x'] / mm + sl2_transf_match_norm['y'] = sl2_transf_match_norm['y'] / mm + transf, _ = align.find_transform(sl2_match_norm, sl2_transf_match_norm, + sl1_match_norm, transModel=model, + order=poly_order, verbose=verbose) + c_exp = np.zeros(len(transf.px._parameters)) + + for i_c in range(len(transf.px._parameters)): + c_exp[i_c] = int(transf.px._param_names[i_c][1:].split('_')[0]) +\ + int(transf.px._param_names[i_c][1:].split('_')[1]) + + c_corr = mm ** (1 - c_exp) + transf.px._parameters = transf.px._parameters * c_corr + transf.py._parameters = transf.py._parameters * c_corr""" + + # Do the final transformation and matching using + sl2_idx, sl1_idx = transform_and_match(sl2, sl1, transf, dr_tol=dr_final, + verbose=verbose) + # StarTable output + sl2_transf = transform_from_object(sl2, transf) + unames = np.array(range(len(sl1_idx))) + st = startables.StarTable(name=unames, + x=np.column_stack((np.array(sl1['x'][sl1_idx]), np.array(sl2_transf['x'][sl2_idx]))), + y=np.column_stack((np.array(sl1['y'][sl1_idx]), np.array(sl2_transf['y'][sl2_idx]))), + m=np.column_stack((np.array(sl1['m'][sl1_idx]), np.array(sl2_transf['m'][sl2_idx]))), + ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx])))) +# ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))), +# list_times=[sl1.meta['list_time'], sl2.meta['list_time']], +# list_names=[sl1.meta['list_name'], sl2.meta['list_name']]) + + for col in sl1.colnames: + if col in sl2.colnames: + if col not in ['name', 'x', 'y', 'm']: + st.add_column(Column(np.column_stack((np.array(sl1[col][sl1_idx]),np.array(sl2_transf[col][sl2_idx]))), name=col)) + + return transf, st diff --git a/flystar/analysis.py b/flystar/analysis.py index 966b913..dc8b61d 100644 --- a/flystar/analysis.py +++ b/flystar/analysis.py @@ -1,15 +1,12 @@ import numpy as np import pylab as plt -from . import starlists, startables, align, match, transforms +from flystar import starlists, match from astropy import table from astropy.table import Table, Column from astropy.coordinates import SkyCoord from astropy import units as u -from astropy.wcs import WCS from astroquery.gaia import Gaia -from astroquery.mast import Observations, Catalogs import pdb, copy -import math from scipy.stats import f ################################################## @@ -471,85 +468,6 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False): # Old codes. ################################################## -def calc_chi2(ref_mat, starlist_mat, transform, errs='both'): - """ - calculate the chi2 and reduced chi2 of the position - between two matched starlists. - Input: - ref_mat: astropy table - Reference starlist only containing matched stars that were used in the - transformation. Standard column headers are assumed. - - starlist_mat: astropy table - Transformed starlist only containing the matched stars used in - the transformation. Standard column headers are assumed. - - transform: transformation object - Transformation object of final transform. Used in chi-square - determination - - errs: string; 'both', 'reference', or 'starlist' - If both, add starlist errors in quadrature with reference errors. - - If reference, only consider reference errors. This should be used if the starlist - does not have valid errors - - If starlist, only consider starlist errors. This should be used if the reference - does not have valid errors - - Output: - chi_sq: float - chi2 = sum (diff_x**2 / xerr**2 + diff_y**2 /yerr**2) - chi_sq_red: float - reduced chi2 = chi2/ degree of freedom - deg_freedom: int - degree of freedom - - """ - diff_x = ref_mat['x'] - starlist_mat['x'] - diff_y = ref_mat['y'] - starlist_mat['y'] - - # Set errors as per user input - if errs == 'both': - xerr = np.hypot(ref_mat['xe'], starlist_mat['xe']) - yerr = np.hypot(ref_mat['ye'], starlist_mat['ye']) - elif errs == 'reference': - xerr = ref_mat['xe'] - yerr = ref_mat['ye'] - elif errs == 'starlist': - xerr = starlist_mat['xe'] - yerr = starlist_mat['ye'] - - - # For both X and Y, calculate chi-square. Combine arrays to get combined - # chi-square - chi_sq_x = diff_x**2. / xerr**2. - chi_sq_y = diff_y**2. / yerr**2. - - chi_sq = np.append(chi_sq_x, chi_sq_y) - - # Calculate degrees of freedom in transformation - num_mod_params = calc_nparam(transform) - deg_freedom = len(chi_sq) - num_mod_params - - # Calculate reduced chi-square - chi_sq = np.sum(chi_sq) - chi_sq_red = chi_sq / deg_freedom - - return chi_sq, chi_sq_red, deg_freedom - - -def calc_nparam(transformation): - """ - calculate the degree of freedom for a transformation - """ - # Read transformation: Extract X, Y coefficients from transform - if transformation.__class__.__name__ == 'four_paramNW': - nparam = 4 - elif transformation.__class__.__name__ == 'PolyTransform': - order = transformation.order - nparam = (order+1) * (order+2) - return nparam def calc_F(red_chi2_1, red_chi2_2, v1, v2): """ diff --git a/flystar/match.py b/flystar/match.py index fe749f0..55b707e 100644 --- a/flystar/match.py +++ b/flystar/match.py @@ -1,5 +1,4 @@ import numpy as np -from . import starlists, transforms, startables, align from collections import Counter from scipy.spatial import cKDTree as KDT from astropy.table import Column, Table @@ -462,210 +461,3 @@ def add_votes(votes, match1, match2): votes.flat[unique_idx] += deltas return - - -def generic_match(sl1, sl2, init_mode='triangle', - model=transforms.PolyTransform, order_dr=(1, 1.0), - dr_final=1.0, - xy_match=(None, None, None, None, None, None, None, None), - m_match=(None, None, None, None), sigma_match=None, - n_bright=100, verbose=True, **kwargs): - """ - Finds the transformation between two starlists using the first one - as reference frame. Different matching methods can be used. If no - transformation is found, it returns an error message. - - - Parameters - sl1 : StarList - starlist used for reference frame - sl2 : StarList - starlist transformed - init_mode : str - Initial matching method. - If 'triangle', uses the blind triangle method. - If 'match_name', uses match by name - If 'load', uses the transformation from a loaded file - model : str - Transformation model to be used with the 'triangle' initial mode - poly_order : int - Order of the transformation model - order_dr : int, float [n, 2] - Combinations of polinomial order (first column) and search radius - (second column) to refine the transformation. Rows are executed in - orders - dr_final: float - Search radius used for the final matching - n_bright : int - Number of bright stars used in the initial blind triangles matching - xy_match : array - Area of the images to remove in the matching [reference catalog min x, - reference catalog max x, reference catalog min y, reference catalog max y, - transformed catalog min x, transformed catalog max x, - transformed catalog min y, transformed catalog max y]. Use None for values not used. - m_match : array - Magnitude limits of matching stars used to find transformations - [reference catalog min mag, reference catalog max mag, transformed - catalog min mag, transformed catalog max mag]. Use None for values not - used - sigma_match : array - Number of Deltap movement sigmas [0] used for sigma-cutting matched - stars for a number of times [1]. Use None for no sigma-cut. The last - polynomial order and search radius in 'order_dr' are used - transf_file : str - File name and path of the transformation file used with the 'load' - init_mode - verbose : bool, optional - Prints on screen information on the matching - - Returns - ------- - transf : Transform2D - Transformation of the second starlist respect to the first - st : StarTable - Startable of the two matched catalogs - - """ - - # Check the input StarLists and transform them into astropy Tables - if not isinstance(sl1, starlists.StarList): - raise TypeError("The first catalog has to be a StarList") - if not isinstance(sl2, starlists.StarList): - raise TypeError("The second catalog has to be a StarList") - - # Find the initial transformation - if init_mode == 'triangle': # Blind triangles method - - # Prepare the reduced starlists for matching - sl1_cut = copy.deepcopy(sl1) - sl2_cut = copy.deepcopy(sl2) - sl1_cut.restrict_by_value(x_min=xy_match[0], x_max=xy_match[1], - y_min=xy_match[2], y_max=xy_match[3]) - sl2_cut.restrict_by_value(x_min=xy_match[4], x_max=xy_match[5], - y_min=xy_match[6], y_max=xy_match[7]) - sl1_cut.restrict_by_value(m_min=m_match[0], m_max=m_match[1]) - sl2_cut.restrict_by_value(m_min=m_match[2], m_max=m_match[3]) - - # Find the transformation - # TODO: test 'initial_align' with StarList input - transf = align.initial_align(sl1_cut, sl2_cut, briteN=n_bright, - transformModel=model, order=order_dr[0]) #order_dr[i_loop][0] ? - - elif init_mode == 'match_name': # Name match - sl1_idx_init, sl2_idx_init, _ = starlists.restrict_by_name(sl1, sl2) - transf = model(sl2['x'][sl2_idx_init], sl2['y'][sl2_idx_init], - sl1['x'][sl1_idx_init], sl1['y'][sl1_idx_init], - order=int(order_dr[0][0])) - - elif init_mode == 'load': # Load a transformation file - transf = transforms.Transform2D.from_file(kwargs['transf_file']) - - else: # None of the above - raise TypeError("Unrecognized initial matching method") - - # Restrict the matching catalogs - sl1_match = copy.deepcopy(sl1) - sl2_match = copy.deepcopy(sl2) - sl1_match.restrict_by_value(m_min=m_match[0], m_max=m_match[1]) - sl2_match.restrict_by_value(m_min=m_match[2], m_max=m_match[3]) - - # Refine the transformation - if sigma_match: - order_dr_len = len(order_dr) - - for i_loop in range(sigma_match[1]): - order_dr = np.vstack((np.array(order_dr), np.array(order_dr[-1]))) - - for i_loop in range(len(order_dr)): - - # Transform and match the catalog to the reference frame -# sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf, -# dr_tol=order_dr[i_loop][1], -# verbose=verbose) - - sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf, - dr_tol=order_dr[1], - verbose=verbose) - - # Transform the catalog to the reference frame - sl2_transf_match = align.transform_from_object(sl2_match, transf) - - # Sigma-rejection - if sigma_match and (i_loop >= order_dr_len): - resid = np.sqrt((sl1_match['x'][sl1_idx] - - sl2_transf_match['x'][sl2_idx])**2 + - (sl1_match['y'][sl1_idx] - - sl2_transf_match['y'][sl2_idx])**2) - sl1_idx = sl1_idx[resid <= (sigma_match[0] * np.std(resid))] - sl2_idx = sl2_idx[resid <= (sigma_match[0] * np.std(resid))] - - # Test section to observe the matching catalogs before refining the transformation - """ - from matplotlib import pyplot - - _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10,10)) - axarr.scatter(sl1_match['x'][sl1_idx], sl1_match['y'][sl1_idx]) - xlim = axarr.get_xlim() - ylim = axarr.get_ylim() - - _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10, 10)) - axarr.scatter(sl2_transf_match['x'][sl2_idx], sl2_transf_match['y'][sl2_idx]) - axarr.set_xlim(xlim) - axarr.set_ylim(ylim) - """ - - # Find a better transformation - transf, _ = align.find_transform(sl2_match[sl2_idx], - sl2_transf_match[sl2_idx], - sl1_match[sl1_idx], transModel=model, - order=order_dr[0], verbose=verbose) -# order=int(order_dr[i_loop][0]), verbose=verbose) - - # This section was used for testing transformations with normalized - # coordinates. Only several catalogs had reduced residuals when using - # high order polynomials (>3), some of them became unstable - """sl1_match_norm = sl1_match[sl1_idx] - sl2_match_norm = sl2_match[sl2_idx] - sl2_transf_match_norm = sl2_transf_match[sl2_idx] - mm = max(max(sl1_match_norm['x']), max(sl1_match_norm['y']), - max(sl2_transf_match_norm['x']), max(sl2_transf_match_norm['y'])) - sl1_match_norm['x'] = sl1_match_norm['x'] / mm - sl1_match_norm['y'] = sl1_match_norm['y'] / mm - sl2_match_norm['x'] = sl2_match_norm['x'] / mm - sl2_match_norm['y'] = sl2_match_norm['y'] / mm - sl2_transf_match_norm['x'] = sl2_transf_match_norm['x'] / mm - sl2_transf_match_norm['y'] = sl2_transf_match_norm['y'] / mm - transf, _ = align.find_transform(sl2_match_norm, sl2_transf_match_norm, - sl1_match_norm, transModel=model, - order=poly_order, verbose=verbose) - c_exp = np.zeros(len(transf.px._parameters)) - - for i_c in range(len(transf.px._parameters)): - c_exp[i_c] = int(transf.px._param_names[i_c][1:].split('_')[0]) +\ - int(transf.px._param_names[i_c][1:].split('_')[1]) - - c_corr = mm ** (1 - c_exp) - transf.px._parameters = transf.px._parameters * c_corr - transf.py._parameters = transf.py._parameters * c_corr""" - - # Do the final transformation and matching using - sl2_idx, sl1_idx = align.transform_and_match(sl2, sl1, transf, dr_tol=dr_final, - verbose=verbose) - # StarTable output - sl2_transf = align.transform_from_object(sl2, transf) - unames = np.array(range(len(sl1_idx))) - st = startables.StarTable(name=unames, - x=np.column_stack((np.array(sl1['x'][sl1_idx]), np.array(sl2_transf['x'][sl2_idx]))), - y=np.column_stack((np.array(sl1['y'][sl1_idx]), np.array(sl2_transf['y'][sl2_idx]))), - m=np.column_stack((np.array(sl1['m'][sl1_idx]), np.array(sl2_transf['m'][sl2_idx]))), - ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx])))) -# ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))), -# list_times=[sl1.meta['list_time'], sl2.meta['list_time']], -# list_names=[sl1.meta['list_name'], sl2.meta['list_name']]) - - for col in sl1.colnames: - if col in sl2.colnames: - if col not in ['name', 'x', 'y', 'm']: - st.add_column(Column(np.column_stack((np.array(sl1[col][sl1_idx]),np.array(sl2_transf[col][sl2_idx]))), name=col)) - - return transf, st diff --git a/flystar/motion_model.py b/flystar/motion_model.py index 3564098..b789cb3 100644 --- a/flystar/motion_model.py +++ b/flystar/motion_model.py @@ -1,7 +1,7 @@ import numpy as np from abc import ABC import pdb -from . import parallax +from flystar import parallax from astropy.time import Time from scipy.optimize import curve_fit import warnings @@ -83,18 +83,7 @@ def get_weights(self, xe, ye, weighting='var'): else: warnings.warn("Invalid weighting, using default weighting scheme var.", UserWarning) 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': - return errs - else: - warnings.warn("Invalid weighting, using default weighting scheme var.", UserWarning) - return errs + def fit_motion_model(self, t, x, y, xe, ye, t0, bootstrap=0, weighting='var', use_scipy=True, absolute_sigma=True): @@ -240,6 +229,8 @@ def get_batch_pos_at_time(self, t, x0=[],vx=[], y0=[],vy=[], t0=[], def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, use_scipy=True, absolute_sigma=True): dt = t-t0 + sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting) + x_wt, y_wt = 1. / sigma_x**2, 1. / sigma_y**2 if params_guess is None: params_guess = [x.mean(),0.0,y.mean(),0.0] @@ -269,14 +260,13 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, if use_scipy: def linear(t, c0, c1): return c0 + c1*t - sigma_x, sigma_y = self.calc_sigma(xe, ye, weighting=weighting) x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=sigma_x, absolute_sigma=absolute_sigma) y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=sigma_y, absolute_sigma=absolute_sigma) x0, vx = x_opt y0, vy = y_opt x0e, vxe = np.sqrt(x_cov.diagonal()) y0e, vye = np.sqrt(y_cov.diagonal()) - x0e, vxe, y0e, vye = self.scale_errors([x0e, vxe, y0e, vye], weighting=weighting) + else: # Use https://en.wikipedia.org/wiki/Weighted_least_squares#Solution scheme x = np.array(x) @@ -300,7 +290,6 @@ def linear(t, c0, c1): y0, vy = popt_y[1], popt_y[0] x0e, vxe = perr_x[1], perr_x[0] y0e, vye = perr_y[1], perr_y[0] - x0e, vxe, y0e, vye = self.scale_errors([x0e, vxe, y0e, vye], weighting=weighting) params = [x0, vx, y0, vy] param_errors = [x0e, vxe, y0e, vye] @@ -370,7 +359,6 @@ def accel(t, c0,c1,c2): x0e, vx0e, axe = np.sqrt(x_cov.diagonal()) y0e, vy0e, aye = np.sqrt(y_cov.diagonal()) - x0e, vx0e, axe, y0e, vy0e, aye = self.scale_errors([x0e, vx0e, axe, y0e, vy0e, aye], weighting=weighting) params = [x0, vx0, ax, y0, vy0, ay] param_errors = [x0e, vx0e, axe, y0e, vy0e, aye] @@ -480,7 +468,7 @@ def fit_func(use_t, x0,vx, y0,vy, pi): res = curve_fit(fit_func, t, np.hstack([x,y]), p0=params_guess, sigma=sigma, absolute_sigma=absolute_sigma) x0,vx,y0,vy,pi = res[0] - x0_err,vx_err,y0_err,vy_err,pi_err = self.scale_errors(np.sqrt(np.diag(res[1])), weighting=weighting) + x0_err,vx_err,y0_err,vy_err,pi_err = np.sqrt(res[1].diagonal()) params = [x0, vx, y0, vy, pi] param_errors = [x0_err, vx_err, y0_err, vy_err, pi_err] diff --git a/flystar/plots.py b/flystar/plots.py index 7d4b700..77642d7 100755 --- a/flystar/plots.py +++ b/flystar/plots.py @@ -1,4 +1,4 @@ -from . import analysis, motion_model, startables +from flystar import motion_model import pylab as py import pylab as plt import numpy as np @@ -17,6 +17,89 @@ from astropy.coordinates import SkyCoord from astropy import units as u + +# Moved here from analysis +def calc_chi2(ref_mat, starlist_mat, transform, errs='both'): + """ + calculate the chi2 and reduced chi2 of the position + between two matched starlists. + Input: + ref_mat: astropy table + Reference starlist only containing matched stars that were used in the + transformation. Standard column headers are assumed. + + starlist_mat: astropy table + Transformed starlist only containing the matched stars used in + the transformation. Standard column headers are assumed. + + transform: transformation object + Transformation object of final transform. Used in chi-square + determination + + errs: string; 'both', 'reference', or 'starlist' + If both, add starlist errors in quadrature with reference errors. + + If reference, only consider reference errors. This should be used if the starlist + does not have valid errors + + If starlist, only consider starlist errors. This should be used if the reference + does not have valid errors + + Output: + chi_sq: float + chi2 = sum (diff_x**2 / xerr**2 + diff_y**2 /yerr**2) + chi_sq_red: float + reduced chi2 = chi2/ degree of freedom + deg_freedom: int + degree of freedom + + """ + diff_x = ref_mat['x'] - starlist_mat['x'] + diff_y = ref_mat['y'] - starlist_mat['y'] + + # Set errors as per user input + if errs == 'both': + xerr = np.hypot(ref_mat['xe'], starlist_mat['xe']) + yerr = np.hypot(ref_mat['ye'], starlist_mat['ye']) + elif errs == 'reference': + xerr = ref_mat['xe'] + yerr = ref_mat['ye'] + elif errs == 'starlist': + xerr = starlist_mat['xe'] + yerr = starlist_mat['ye'] + + + # For both X and Y, calculate chi-square. Combine arrays to get combined + # chi-square + chi_sq_x = diff_x**2. / xerr**2. + chi_sq_y = diff_y**2. / yerr**2. + + chi_sq = np.append(chi_sq_x, chi_sq_y) + + # Calculate degrees of freedom in transformation + num_mod_params = calc_nparam(transform) + deg_freedom = len(chi_sq) - num_mod_params + + # Calculate reduced chi-square + chi_sq = np.sum(chi_sq) + chi_sq_red = chi_sq / deg_freedom + + return chi_sq, chi_sq_red, deg_freedom + + +def calc_nparam(transformation): + """ + calculate the degree of freedom for a transformation + """ + # Read transformation: Extract X, Y coefficients from transform + if transformation.__class__.__name__ == 'four_paramNW': + nparam = 4 + elif transformation.__class__.__name__ == 'PolyTransform': + order = transformation.order + nparam = (order+1) * (order+2) + return nparam + + #################################################### # Code for making diagnostic plots for astrometry # alignment @@ -226,15 +309,15 @@ def pos_diff_err_hist(ref_mat, starlist_mat, transform, nbins=25, bin_width=None chi_sq_red = np.sum(chi_sq) / deg_freedom """ # Chi-square analysis for all stars, including outliers - chi_sq, chi_sq_red, deg_freedom = analysis.calc_chi2(ref_mat, starlist_mat, + chi_sq, chi_sq_red, deg_freedom = calc_chi2(ref_mat, starlist_mat, transform, errs=errs) # Chi-square analysis for only non-outlier stars - chi_sq_good, chi_sq_red_good, deg_freedom_good = analysis.calc_chi2(ref_mat[good], + chi_sq_good, chi_sq_red_good, deg_freedom_good = calc_chi2(ref_mat[good], starlist_mat[good], transform, errs=errs) - num_mod_params = analysis.calc_nparam(transform) + num_mod_params = calc_nparam(transform) #-------------------------------------------# # Plotting @@ -262,7 +345,7 @@ def pos_diff_err_hist(ref_mat, starlist_mat, transform, nbins=25, bin_width=None py.plot(x, norm.pdf(x,mean,sigma), 'g-', linewidth=2) # Annotate reduced chi-sqared values in plot: with outliers - xstr = '$\chi^2_r$ = {0}'.format(np.round(chi_sq_red, decimals=3)) + xstr = r'$\chi^2_r$ = {0}'.format(np.round(chi_sq_red, decimals=3)) py.annotate(xstr, xy=(0.3, 0.77), xycoords='figure fraction', color='black') txt = r'$\nu$ = 2*{0} - {1} = {2}'.format(len(diff_x), num_mod_params, deg_freedom) @@ -273,7 +356,7 @@ def pos_diff_err_hist(ref_mat, starlist_mat, transform, nbins=25, bin_width=None py.annotate(xstr3, xy=(0.25, 0.80), xycoords='figure fraction', color='black') # Annotate reduced chi-sqared values in plot: without outliers - xstr = '$\chi^2_r$ = {0}'.format(np.round(chi_sq_red_good, decimals=3)) + xstr = r'$\chi^2_r$ = {0}'.format(np.round(chi_sq_red_good, decimals=3)) py.annotate(xstr, xy=(0.7, 0.8), xycoords='figure fraction', color='black') txt = r'$\nu$ = 2*{0} - {1} = {2}'.format(len(good[0]), num_mod_params, deg_freedom_good) @@ -2226,8 +2309,8 @@ def plot_chi2_dist(tab, Ndetect, motion_model_dict={}, xlim=40, n_bins=50, boot_ plt.hist(x[idx], bins=chi2_bins, histtype='step', label='X', density=True) plt.hist(y[idx], bins=chi2_bins, histtype='step', label='Y', density=True) plt.plot(chi2_xaxis, chi2.pdf(chi2_xaxis, Ndof), 'r-', alpha=0.6, - label='$\chi^2$ ' + str(round(Ndof,2)) + ' dof') - plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(round(Ndof,2))) + label=r'$\chi^2$ ' + str(round(Ndof,2)) + ' dof') + plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(round(Ndof,2))) plt.xlim(0, xlim) plt.legend() @@ -2390,8 +2473,8 @@ def plot_chi2_dist_per_filter(tab, Ndetect, motion_model_dict={}, xlim=40, n_bin plt.hist(x[idx], bins=chi2_bins, histtype='stepfilled', label='RA', density=True, color='skyblue', alpha=0.8, edgecolor='k') plt.hist(y[idx], bins=chi2_bins, histtype='stepfilled', label='DEC', density=True, color='orange', alpha=0.8, edgecolor='k') plt.plot(chi2_xaxis, chi2.pdf(chi2_xaxis, Ndof), 'r-', alpha=0.6, - label='$\chi^2$ ' + str(Ndof) + ' dof') - #plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) + label=r'$\chi^2$ ' + str(Ndof) + ' dof') + #plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) plt.title(str(filter)+' (N = '+str(len(chi2_x_list))+')', fontsize=22) plt.xlim(0, xlim) plt.ylabel(r'PDF', fontsize=28) @@ -2677,8 +2760,8 @@ def plot_chi2_dist_mag(tab, Ndetect, xlim=40, n_bins=30, boot_err=False): plt.clf() plt.hist(chi2_m[idx], bins=np.arange(xlim*10), histtype='step', density=True) plt.plot(chi2_maxis, chi2.pdf(chi2_maxis, Ndof), 'r-', alpha=0.6, - label='$\chi^2$ ' + str(Ndof) + ' dof') - plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) + label=r'$\chi^2$ ' + str(Ndof) + ' dof') + plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) plt.xlim(0, xlim) plt.legend() @@ -2726,8 +2809,8 @@ def plot_chi2_dist_mag_per_filter(tab, Ndetect, mlim=40, n_bins=30, xlim=40, fil plt.clf() plt.hist(chi2_m[idx], bins=np.arange(xlim*10), label='mag', histtype='stepfilled', density=True, color='green', alpha=0.7, edgecolor='k') plt.plot(chi2_maxis, chi2.pdf(chi2_maxis, Ndof), 'r-', alpha=0.6, - label='$\chi^2$ ' + str(Ndof) + ' dof') - #plt.title('$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) + label=r'$\chi^2$ ' + str(Ndof) + ' dof') + #plt.title(r'$N_{epoch} = $' + str(Ndetect) + ', $N_{dof} = $' + str(Ndof)) plt.xlim(0, xlim) plt.xlabel(r'$\chi^{2}$', fontsize=28) plt.ylabel(r'PDF', fontsize=28) diff --git a/flystar/transforms.py b/flystar/transforms.py index 4f4410c..6cc865a 100755 --- a/flystar/transforms.py +++ b/flystar/transforms.py @@ -6,7 +6,7 @@ import collections import re import pdb -from . import motion_model +from flystar import motion_model class Transform2D(object): ''' From 621cb3f99fb2af708f8eca67e391187caa546e27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Thu, 6 Aug 2026 16:43:23 -0700 Subject: [PATCH 10/12] Update test data path; Fixed gaia query error --- flystar/analysis.py | 2 +- flystar/tests/test_align.py | 119 ++++++++-------------- flystar/tests/{ => test_data}/A.lis | 0 flystar/tests/{ => test_data}/B.lis | 0 flystar/tests/{ => test_data}/C.lis | 0 flystar/tests/{ => test_data}/D.lis | 0 flystar/tests/{ => test_data}/E.lis | 0 flystar/tests/{ => test_data}/F.lis | 0 flystar/tests/{ => test_data}/ref.lis | 0 flystar/tests/{ => test_data}/ref_vel.lis | 0 10 files changed, 44 insertions(+), 77 deletions(-) rename flystar/tests/{ => test_data}/A.lis (100%) rename flystar/tests/{ => test_data}/B.lis (100%) rename flystar/tests/{ => test_data}/C.lis (100%) rename flystar/tests/{ => test_data}/D.lis (100%) rename flystar/tests/{ => test_data}/E.lis (100%) rename flystar/tests/{ => test_data}/F.lis (100%) rename flystar/tests/{ => test_data}/ref.lis (100%) rename flystar/tests/{ => test_data}/ref_vel.lis (100%) diff --git a/flystar/analysis.py b/flystar/analysis.py index dc8b61d..1a2ea82 100644 --- a/flystar/analysis.py +++ b/flystar/analysis.py @@ -42,7 +42,7 @@ def query_gaia(ra, dec, search_radius=30.0, table_name='gaiadr3'): search_radius *= u.arcsec Gaia.ROW_LIMIT = 50000 - gaia_job = Gaia.cone_search_async(target_coords, search_radius, table_name = table_name + '.gaia_source') + gaia_job = Gaia.cone_search_async(target_coords, radius=search_radius, table_name=table_name + '.gaia_source') gaia = gaia_job.get_results() #Change new 'SOURCE_ID' column header back to lowercase 'source_id' so all subsequent functions still work: diff --git a/flystar/tests/test_align.py b/flystar/tests/test_align.py index 2d6b0dc..d64ff9f 100644 --- a/flystar/tests/test_align.py +++ b/flystar/tests/test_align.py @@ -1,21 +1,16 @@ -from flystar import align -from flystar import starlists -from flystar import startables -from flystar import transforms -from flystar import analysis -from flystar import motion_model -from astropy.table import Table import numpy as np import pylab as plt -import pdb -import datetime -import pytest +import flystar +from astropy.table import Table +from flystar import align, starlists, transforms, analysis, motion_model + +test_data_path = f'{flystar.__path__[0]}/tests/test_data' def test_MosaicSelfRef(): """ Cross-match and align 4 starlists using the OO version of mosaic lists. """ - list_files = ['A.lis', 'B.lis', 'C.lis', 'D.lis'] + list_files = [f'{test_data_path}/{f}' for f in ['A.lis', 'B.lis', 'C.lis', 'D.lis']] lists = [starlists.StarList.from_lis_file(lf) for lf in list_files] ########## @@ -91,7 +86,7 @@ def test_MosaicSelfRef_vel_tconst(): The 4 lists are all taken at the same time (so 0 velocities should result). """ - list_files = ['A.lis', 'B.lis', 'C.lis', 'D.lis'] + list_files = [f'{test_data_path}/{f}' for f in ['A.lis', 'B.lis', 'C.lis', 'D.lis']] lists = [starlists.StarList.from_lis_file(lf) for lf in list_files] ########## @@ -149,7 +144,7 @@ def test_MosaicSelfRef_vel(): Cross-match and align 4 starlists using the OO version of mosaic lists. """ - list_files = ['A.lis', 'B.lis', 'C.lis', 'D.lis'] + list_files = [f'{test_data_path}/{f}' for f in ['A.lis', 'B.lis', 'C.lis', 'D.lis']] lists = [starlists.StarList.from_lis_file(lf) for lf in list_files] # Modify the times so that we get velocities out. @@ -215,15 +210,8 @@ def test_MosaicSelfRef_vel(): def test_MosaicToRef(): make_fake_starlists_poly1(seed=42) - ref_file = 'random_ref.fits' - list_files = ['random_0.fits', - 'random_1.fits', - 'random_2.fits', - 'random_3.fits', - 'random_4.fits', - 'random_5.fits', - 'random_6.fits', - 'random_7.fits'] + ref_file = f'{test_data_path}/random_ref.fits' + list_files = [f'{test_data_path}/random_{i}.fits' for i in range(8)] ref_list = Table.read(ref_file) @@ -272,15 +260,8 @@ def test_MosaicToRef(): def test_MosaicToRef_p0_vel(): make_fake_starlists_poly0_vel(seed=42) - ref_file = 'random_vel_ref.fits' - list_files = ['random_vel_p0_0.fits', - 'random_vel_p0_1.fits', - 'random_vel_p0_2.fits', - 'random_vel_p0_3.fits'] - #'random_vel_4.fits', - #'random_vel_5.fits', - #'random_vel_6.fits', - #'random_vel_7.fits'] + ref_file = f'{test_data_path}/random_vel_ref.fits' + list_files = [f'{test_data_path}/random_vel_p0_{i}.fits' for i in range(4)] ref_list = Table.read(ref_file) @@ -338,15 +319,8 @@ def test_MosaicToRef_p0_vel(): def test_MosaicToRef_vel(): make_fake_starlists_poly1_vel(seed=42) - ref_file = 'random_vel_ref.fits' - list_files = ['random_vel_0.fits', - 'random_vel_1.fits', - 'random_vel_2.fits', - 'random_vel_3.fits'] - #'random_vel_4.fits', - #'random_vel_5.fits', - #'random_vel_6.fits', - #'random_vel_7.fits'] + ref_file = f'{test_data_path}/random_vel_ref.fits' + list_files = [f'{test_data_path}/random_vel_{i}.fits' for i in range(4)] ref_list = Table.read(ref_file) @@ -404,15 +378,8 @@ def test_MosaicToRef_vel(): def test_MosaicToRef_acc(): make_fake_starlists_poly1_acc(seed=42) - ref_file = 'random_acc_ref.fits' - list_files = ['random_acc_0.fits', - 'random_acc_1.fits', - 'random_acc_2.fits', - 'random_acc_3.fits', - 'random_acc_4.fits', - 'random_acc_5.fits', - 'random_acc_6.fits', - 'random_acc_7.fits'] + ref_file = f'{test_data_path}/random_acc_ref.fits' + list_files = [f'{test_data_path}/random_acc_{i}.fits' for i in range(8)] ref_list = Table.read(ref_file) @@ -500,7 +467,7 @@ def make_fake_starlists_shifts(): # Save original positions as reference (1st) list. fmt = '{0:10s} {1:5.2f} 2015.0 {2:9.4f} {3:9.4f} 0 0 0 0\n' - _out = open('random_0.lis', 'w') + _out = open(f'{test_data_path}/random_0.lis', 'w') for ii in range(N_stars): _out.write(fmt.format(name[ii], m[ii], x[ii], y[ii])) _out.close() @@ -525,7 +492,7 @@ def make_fake_starlists_shifts(): mnew = m + np.random.randn(N_stars) * 0.05 - _out = open('random_shift_{0:d}.lis'.format(ss+1), 'w') + _out = open(f'{test_data_path}/random_shift_{ss+1}.lis', 'w') for ii in range(N_stars): _out.write(fmt.format(name[ii], mnew[ii], xnew[ii], ynew[ii])) _out.close() @@ -563,7 +530,7 @@ def make_fake_starlists_poly1(seed=-1): # Save original positions as reference (1st) list # in a StarList format (with velocities). - lis.write('random_ref.fits', overwrite=True) + lis.write(f'{test_data_path}/random_ref.fits', overwrite=True) ########## # Shifts @@ -614,7 +581,7 @@ def make_fake_starlists_poly1(seed=-1): new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t], names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't')) - new_lis.write('random_{0:d}.fits'.format(ss), overwrite=True) + new_lis.write(f'{test_data_path}/random_{ss}.fits', overwrite=True) return (xy_trans,mag_trans) @@ -656,7 +623,7 @@ def make_fake_starlists_poly0_vel(seed=-1): # Save original positions as reference (1st) list # in a StarList format (with velocities). - lis.write('random_vel_ref.fits', overwrite=True) + lis.write(f'{test_data_path}/random_vel_ref.fits', overwrite=True) ########## # Propogate to new times and distort. @@ -707,7 +674,7 @@ def make_fake_starlists_poly0_vel(seed=-1): new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t], names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't')) - new_lis.write('random_vel_p0_{0:d}.fits'.format(ss), overwrite=True) + new_lis.write(f'{test_data_path}/random_vel_p0_{ss}.fits', overwrite=True) return (xy_trans, mag_trans) @@ -750,7 +717,7 @@ def make_fake_starlists_poly1_vel(seed=-1): # Save original positions as reference (1st) list # in a StarList format (with velocities). - lis.write('random_vel_ref.fits', overwrite=True) + lis.write(f'{test_data_path}/random_vel_ref.fits', overwrite=True) ########## # Propogate to new times and distort. @@ -801,7 +768,7 @@ def make_fake_starlists_poly1_vel(seed=-1): new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t], names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't')) - new_lis.write('random_vel_{0:d}.fits'.format(ss), overwrite=True) + new_lis.write(f'{test_data_path}/random_vel_{ss}.fits', overwrite=True) return (xy_trans, mag_trans) @@ -856,7 +823,7 @@ def make_fake_starlists_poly1_acc(seed=-1): # Save original positions as reference (1st) list # in a StarList format (with velocities). - lis.write('random_acc_ref.fits', overwrite=True) + lis.write(f'{test_data_path}/random_acc_ref.fits', overwrite=True) ########## # Propogate to new times and distort. @@ -907,7 +874,7 @@ def make_fake_starlists_poly1_acc(seed=-1): new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t], names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't')) - new_lis.write('random_acc_{0:d}.fits'.format(ss), overwrite=True) + new_lis.write(f'{test_data_path}/random_acc_{ss}.fits', overwrite=True) return (xy_trans, mag_trans) @@ -959,7 +926,7 @@ def make_fake_starlists_poly1_par(seed=-1): # Save original positions as reference (1st) list # in a StarList format (with velocities). - lis.write('random_par_ref.fits', overwrite=True) + lis.write(f'{test_data_path}/random_par_ref.fits', overwrite=True) ########## # Propogate to new times and distort. @@ -1019,7 +986,7 @@ def make_fake_starlists_poly1_par(seed=-1): new_lis = starlists.StarList([lis['name'], md, mde, xd, xde, yd, yde, t], names=('name', 'm', 'me', 'x', 'xe', 'y', 'ye', 't')) - new_lis.write('random_par_{0:d}.fits'.format(ss), overwrite=True) + new_lis.write(f'{test_data_path}/random_par_{ss}.fits', overwrite=True) return (xy_trans, mag_trans) @@ -1036,15 +1003,15 @@ def test_MosaicToRef_hst_me(): dec = '-34:27:05.01' # Load up a Gaia catalog (queried around the RA/Dec above) - my_gaia = Table.read('mb10364_data/my_gaia.fits') + my_gaia = Table.read(f'{test_data_path}/my_gaia.fits') my_gaia['me'] = 0.01 # Gather the list of starlists. For first pass, don't modify the starlists. # Loop through the observations and read them in, in prep for alignment with Gaia epochs = [2011.83, 2012.73, 2013.81] - starlist_names = ['mb10364_data/2011_10_31_F606W_MATCHUP_XYMEEE_final.calib', - 'mb10364_data/2012_09_25_F606W_MATCHUP_XYMEEE_final.calib', - 'mb10364_data/2013_10_24_F606W_MATCHUP_XYMEEE_final.calib'] + starlist_names = [f'{test_data_path}/mb10364_data/2011_10_31_F606W_MATCHUP_XYMEEE_final.calib', + f'{test_data_path}/mb10364_data/2012_09_25_F606W_MATCHUP_XYMEEE_final.calib', + f'{test_data_path}/mb10364_data/2013_10_24_F606W_MATCHUP_XYMEEE_final.calib'] list_of_starlists = [] @@ -1090,9 +1057,9 @@ def test_bootstrap(): etc.) """ # Read in starlists for MosaicToRef - ref = Table.read('ref_vel.lis', format='ascii') - list1 = Table.read('E.lis', format='ascii') - list2 = Table.read('F.lis', format='ascii') + ref = Table.read(f'{test_data_path}/ref_vel.lis', format='ascii') + list1 = Table.read(f'{test_data_path}/E.lis', format='ascii') + list2 = Table.read(f'{test_data_path}/F.lis', format='ascii') list1 = starlists.StarList.from_table(list1) list2 = starlists.StarList.from_table(list2) @@ -1202,10 +1169,10 @@ def test_calc_vel_in_bootstrap(): import copy # Define match parameters - ref = Table.read('ref_vel.lis', format='ascii') + ref = Table.read(f'{test_data_path}/ref_vel.lis', format='ascii') - list1 = Table.read('E.lis', format='ascii') - list2 = Table.read('F.lis', format='ascii') + list1 = Table.read(f'{test_data_path}/E.lis', format='ascii') + list2 = Table.read(f'{test_data_path}/F.lis', format='ascii') list1 = starlists.StarList.from_table(list1) list2 = starlists.StarList.from_table(list2) @@ -1271,9 +1238,9 @@ def test_transform_xym(): otherwise """ #---Align 1: self.mag_Trans = False---# - ref = Table.read('ref_vel.lis', format='ascii') - list1 = Table.read('E.lis', format='ascii') - list2 = Table.read('F.lis', format='ascii') + ref = Table.read(f'{test_data_path}/ref_vel.lis', format='ascii') + list1 = Table.read(f'{test_data_path}/E.lis', format='ascii') + list2 = Table.read(f'{test_data_path}/F.lis', format='ascii') list1 = starlists.StarList.from_table(list1) list2 = starlists.StarList.from_table(list2) @@ -1369,7 +1336,7 @@ def test_MosaicToRef_mag_bug(): """ make_fake_starlists_poly1_vel() - ref_list = starlists.StarList.read('random_vel_0.fits') + ref_list = starlists.StarList.read(f'{test_data_path}/random_vel_0.fits') lists = [ref_list] msc = align.MosaicToRef(ref_list, lists, @@ -1432,7 +1399,7 @@ def test_masked_cols(): list_of_starlists = [] for ee in range(len(epochs)): - lis_file = 'mag' + epochs[ee] + '_ob150029_kp_rms_named.lis' + lis_file = f'{test_data_path}/mag{epochs[ee]}_ob150029_kp_rms_named.lis' lis = starlists.StarList.from_lis_file(lis_file) list_of_starlists.append(lis) diff --git a/flystar/tests/A.lis b/flystar/tests/test_data/A.lis similarity index 100% rename from flystar/tests/A.lis rename to flystar/tests/test_data/A.lis diff --git a/flystar/tests/B.lis b/flystar/tests/test_data/B.lis similarity index 100% rename from flystar/tests/B.lis rename to flystar/tests/test_data/B.lis diff --git a/flystar/tests/C.lis b/flystar/tests/test_data/C.lis similarity index 100% rename from flystar/tests/C.lis rename to flystar/tests/test_data/C.lis diff --git a/flystar/tests/D.lis b/flystar/tests/test_data/D.lis similarity index 100% rename from flystar/tests/D.lis rename to flystar/tests/test_data/D.lis diff --git a/flystar/tests/E.lis b/flystar/tests/test_data/E.lis similarity index 100% rename from flystar/tests/E.lis rename to flystar/tests/test_data/E.lis diff --git a/flystar/tests/F.lis b/flystar/tests/test_data/F.lis similarity index 100% rename from flystar/tests/F.lis rename to flystar/tests/test_data/F.lis diff --git a/flystar/tests/ref.lis b/flystar/tests/test_data/ref.lis similarity index 100% rename from flystar/tests/ref.lis rename to flystar/tests/test_data/ref.lis diff --git a/flystar/tests/ref_vel.lis b/flystar/tests/test_data/ref_vel.lis similarity index 100% rename from flystar/tests/ref_vel.lis rename to flystar/tests/test_data/ref_vel.lis From ab7b36269ed57381ae3301ef4c874f86faaded49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CLingfeng?= Date: Thu, 6 Aug 2026 17:40:56 -0700 Subject: [PATCH 11/12] Fixing Linear model --- flystar/align.py | 8 +-- flystar/motion_model.py | 108 ++++++++++++++++++------------------ flystar/tests/test_align.py | 7 ++- 3 files changed, 63 insertions(+), 60 deletions(-) diff --git a/flystar/align.py b/flystar/align.py index ff08ba6..eabe887 100755 --- a/flystar/align.py +++ b/flystar/align.py @@ -204,10 +204,10 @@ def = None. If not None, then this should contain an array or list of transform self.verbose = verbose # For backwards compatibility. - if self.verbose is True: - self.verbose = 9 - if self.verbose is False: - self.verbose = 0 + # if self.verbose is True: + # self.verbose = 9 + # if self.verbose is False: + # self.verbose = 0 self.N_lists = len(self.star_lists) diff --git a/flystar/motion_model.py b/flystar/motion_model.py index b789cb3..c4964d1 100644 --- a/flystar/motion_model.py +++ b/flystar/motion_model.py @@ -234,62 +234,62 @@ def run_fit(self, t, x, y, xe, ye, t0, weighting='var', params_guess=None, if params_guess is None: params_guess = [x.mean(),0.0,y.mean(),0.0] - # Handle 2-data point case - if len(np.unique(dt))==2: - if len(x)>2: # Catch case where bootstrap sends only 2 unique epochs - _,idx=np.unique(dt, return_index=True) - dt = dt[idx] - x = x[idx] - y = y[idx] - xe = xe[idx] - ye = ye[idx] - dx = np.diff(x)[0] - dy = np.diff(y)[0] - dt_diff = np.diff(dt)[0] - vx = dx / dt_diff - vy = dy / dt_diff - # TODO: still not sure about the error handling here - x0 = x[0] - dt[0]*vx # np.average(x, weights=x_wt) # - y0 = y[0] - dt[0]*vy # np.average(y, weights=y_wt) # - x0e = np.abs(dx) / 2**0.5 # np.sqrt(np.sum(xe**2)/2) # - y0e = np.abs(dy) / 2**0.5 # np.sqrt(np.sum(ye**2)/2) # - vxe = 0.0 #np.abs(vx) * np.sqrt(np.sum(xe**2/x**2)) - vye = 0.0 #np.abs(vy) * np.sqrt(np.sum(ye**2/y**2)) + # # Handle 2-data point case + # if len(np.unique(dt))==2: + # if len(x)>2: # Catch case where bootstrap sends only 2 unique epochs + # _,idx=np.unique(dt, return_index=True) + # dt = dt[idx] + # x = x[idx] + # y = y[idx] + # xe = xe[idx] + # ye = ye[idx] + # dx = np.diff(x)[0] + # dy = np.diff(y)[0] + # dt_diff = np.diff(dt)[0] + # vx = dx / dt_diff + # vy = dy / dt_diff + # # TODO: still not sure about the error handling here + # x0 = x[0] - dt[0]*vx # np.average(x, weights=x_wt) # + # y0 = y[0] - dt[0]*vy # np.average(y, weights=y_wt) # + # x0e = np.abs(dx) / 2**0.5 # np.sqrt(np.sum(xe**2)/2) # + # y0e = np.abs(dy) / 2**0.5 # np.sqrt(np.sum(ye**2)/2) # + # vxe = 0.0 #np.abs(vx) * np.sqrt(np.sum(xe**2/x**2)) + # vye = 0.0 #np.abs(vy) * np.sqrt(np.sum(ye**2/y**2)) - else: - if use_scipy: - def linear(t, c0, c1): - return c0 + c1*t - x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=sigma_x, absolute_sigma=absolute_sigma) - y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=sigma_y, absolute_sigma=absolute_sigma) - x0, vx = x_opt - y0, vy = y_opt - x0e, vxe = np.sqrt(x_cov.diagonal()) - y0e, vye = np.sqrt(y_cov.diagonal()) + # else: + if use_scipy: + def linear(t, c0, c1): + return c0 + c1*t + x_opt, x_cov = curve_fit(linear, dt, x, p0=np.array(params_guess[:2]), sigma=sigma_x, absolute_sigma=absolute_sigma) + y_opt, y_cov = curve_fit(linear, dt, y, p0=np.array(params_guess[2:]), sigma=sigma_y, absolute_sigma=absolute_sigma) + x0, vx = x_opt + y0, vy = y_opt + x0e, vxe = np.sqrt(x_cov.diagonal()) + y0e, vye = np.sqrt(y_cov.diagonal()) - else: - # Use https://en.wikipedia.org/wiki/Weighted_least_squares#Solution scheme - x = np.array(x) - y = np.array(y) - dt = np.array(dt) - X_mat_t = np.vander(dt, 2) - # x calculation - W_mat_x = np.diag(x_wt) - XTWX_mat_x = X_mat_t.T @ W_mat_x @ X_mat_t - pcov_x = np.linalg.inv(XTWX_mat_x) # Covariance Matrix - popt_x = pcov_x @ X_mat_t.T @ W_mat_x @ x # Linear Solution - perr_x = np.sqrt(np.diag(pcov_x)) # Uncertainty of Linear Solution - # y calculation - W_mat_y = np.diag(y_wt) - XTWX_mat_y = X_mat_t.T @ W_mat_y @ X_mat_t - pcov_y = np.linalg.inv(XTWX_mat_y) # Covariance Matrix - popt_y = pcov_y @ X_mat_t.T @ W_mat_y @ y # Linear Solution - perr_y = np.sqrt(np.diag(pcov_y)) # Uncertainty of Linear Solution - # prepare values to return - x0, vx = popt_x[1], popt_x[0] - y0, vy = popt_y[1], popt_y[0] - x0e, vxe = perr_x[1], perr_x[0] - y0e, vye = perr_y[1], perr_y[0] + else: + # Use https://en.wikipedia.org/wiki/Weighted_least_squares#Solution scheme + x = np.array(x) + y = np.array(y) + dt = np.array(dt) + X_mat_t = np.vander(dt, 2) + # x calculation + W_mat_x = np.diag(x_wt) + XTWX_mat_x = X_mat_t.T @ W_mat_x @ X_mat_t + pcov_x = np.linalg.inv(XTWX_mat_x) # Covariance Matrix + popt_x = pcov_x @ X_mat_t.T @ W_mat_x @ x # Linear Solution + perr_x = np.sqrt(np.diag(pcov_x)) # Uncertainty of Linear Solution + # y calculation + W_mat_y = np.diag(y_wt) + XTWX_mat_y = X_mat_t.T @ W_mat_y @ X_mat_t + pcov_y = np.linalg.inv(XTWX_mat_y) # Covariance Matrix + popt_y = pcov_y @ X_mat_t.T @ W_mat_y @ y # Linear Solution + perr_y = np.sqrt(np.diag(pcov_y)) # Uncertainty of Linear Solution + # prepare values to return + x0, vx = popt_x[1], popt_x[0] + y0, vy = popt_y[1], popt_y[0] + x0e, vxe = perr_x[1], perr_x[0] + y0e, vye = perr_y[1], perr_y[0] params = [x0, vx, y0, vy] param_errors = [x0e, vxe, y0e, vye] diff --git a/flystar/tests/test_align.py b/flystar/tests/test_align.py index d64ff9f..4d4cc73 100644 --- a/flystar/tests/test_align.py +++ b/flystar/tests/test_align.py @@ -19,7 +19,7 @@ def test_MosaicSelfRef(): msc = align.MosaicSelfRef(lists, ref_index=0, iters=2, dr_tol=[3, 3], dm_tol=[1, 1], trans_class=transforms.PolyTransform, - verbose=False, + verbose=2, trans_args={'order': 2}) msc.fit() @@ -168,7 +168,7 @@ def test_MosaicSelfRef_vel(): dr_tol=[5, 3, 3], dm_tol=[1, 1, 0.5], outlier_tol=None, trans_class=transforms.PolyTransform, trans_args={'order': 2}, default_motion_model='Linear', - verbose=False) + verbose=2) msc.fit() @@ -1418,3 +1418,6 @@ def test_masked_cols(): msc.fit() return + +if __name__ == "__main__": + test_MosaicSelfRef_vel() \ No newline at end of file From 85a1dbacbc391c9fda4d26762c083b3b609d29ef Mon Sep 17 00:00:00 2001 From: Lingfeng Wei Date: Thu, 6 Aug 2026 19:35:29 -0700 Subject: [PATCH 12/12] Added Parallax test data --- flystar/tests/test_data/list_of_starlists.pkl | Bin 0 -> 70070 bytes flystar/tests/test_data/my_gaia.pkl | Bin 0 -> 8805 bytes .../{ => test_data}/test_all_detected.fits | 0 flystar/tests/{ => test_data}/test_catalog.fits | 0 4 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 flystar/tests/test_data/list_of_starlists.pkl create mode 100644 flystar/tests/test_data/my_gaia.pkl rename flystar/tests/{ => test_data}/test_all_detected.fits (100%) rename flystar/tests/{ => test_data}/test_catalog.fits (100%) diff --git a/flystar/tests/test_data/list_of_starlists.pkl b/flystar/tests/test_data/list_of_starlists.pkl new file mode 100644 index 0000000000000000000000000000000000000000..3662f0f65f88a8a944aaac52225bb815feff41d4 GIT binary patch literal 70070 zcmeEvb$k`a*M5+q1&Xwk0!52Ua_E-58gZ>#g%+t#=J~meu@poJCcYIuY<)HZ3yo1J!964l=*RV08#o4?aJO&T( z7&5r&ut6F_XCWo~Xtz;A;$uG@?LBI&PmMuiJciU5V# z%J#f;FOo~=Y~|YOP}}A@;&l~!7qQbdq=S-a&qun`rg&{_V9!hEXd_;-kL~@*K6DSI zD^c1MFX^LnP`soU=~`Fpv?*S7p6;jq&P*=J%uFs_M{;fTB|a+K>!9RPeQ5I(yO-GA z#O@(>;wL+(`iMXE%@Hr%%YJjjtL$&fpBXRNi1e`EGUKIsl5JJ})isWIscyt?yM2h) zUJpmSW+9Mv+wDWVstpwXhj{I6{2^X@8$05q`x~_b@se-aZjR)Vu10M}x*D|u@fx)Q z@fx)Q@fp>Z_>AgHyhi1cZH>yMzG2@#9km7Vn#FFreTdgc8`Hf>rV+WsOZ76Muj2hD za*fzb$yMX7P)FM`vjg#xKN!`4c#X;>J|nzzZ|avue1~{x+%q!v5iixxNS{}{M(j(x zv>TO6GKt@I^ArUyvD^L*5fAluqje)*qxPkCBU{*Rj&z`QFw!qbU*a>$OS&2v&q%IO z-yz#lI~b9x+QEo@6))M%cFW9W#H;v?@RF{?Z@W2?OZPXbFYy|+1MwQ^Z=@&j8P(UQ z?~pw58zVlcc*!qpH%B%jxkmLh>c1q9e9}mtSG-301+^EouMs;CFYy_%gW@$}2jV4} zM*K|i8r6aL{)rAoZALPU>Oj0kb)fng8B>Xm`j8QSP`pOlnB*GOf%==#b|YSrZM%7j zf|uBBe}_n4qkTy68kK9b50PBrH6oXIjQ$1?FWJk8%@nVZc2m5hx9ygheHAaYhwbKw zm+o&=F7X=i9nyjNkI}jjuTi zKh9b;;#d&>B^^qa`nwh$U^2*q{Pg1FG;v&;EI!S_iJ?=`so%?>{XaPPR! zvv+SE!rurr9@YusLuvaXm>+IaJF4;QV6N5cw3Yjc>fIh_|8LAvD<;$bar44ADIR@{$6fhj%E>`-borLY>-$z0h6_138HsJIoJd zZ$ExUm@{CXKaTtRa$$`D_XY=YYL@KIk+|R40JKj}A6}9SIz0sU?H0+Y`9G`d%bQTe z!b3SV)AF=X-jNJ+I2`s43gy2Bh?YJX%G=TnzVYMRTZt>jOoDx@1o8R8O+1>&cK;}d zlbdu(4@5h)j)cCw!uhPPV@Ce4CIa%ehC=_{{+yaO|EzG{RGfbIFc9`G9Rj~t?uB+Z z7>@g~C{E_h(<_pXq6W$p&U2-d=r#Dd75B*#%00Z@=Ko@J80@;t z8}+RnhWiS?)cmZOKjiL@$_G9=cyjT$gjgVnYfZ?I6o)U z>siqUHmcyqSBb7LFn&Dx=8|C8Dl7o^%jb`_uM~x`$J>jO$+G?!j67bG><>HaJayp3cE2d--C_spbv6Crv;f z%NB{Tq_CCGog=o5-e~ip!MOiHUoJ+JF5jC2(C0r2Mqdc?g%37~04+2D*BuUneLoN8 zG$c)W>w__XPlPW|_QqHx#$xS$#m1vgtoG(*$gmTFFgD-uMLQn~ft}9CzVt8v^_VXE z=!P(isXYT>mvT~{dPJg+{1OR!dqzNClP~`&hww1r>zbd1 z1fag1teg(VRFQpZK_Ke2&kugW{V`TNkbU522=rMP4F4K0`()=}j8Vc*wBM5BVd#5B zL!nPXDBA6emD7;0GA06ke8>-d>wuMO{o!Fa+EMgH%`g3|7~Agn^Op3Za3U1`kr0Ib zQN<7aW7H&^&+Cip7esIxLhAoC0b@@7VASVT5c)y>N$7XABB}2>|E2GaGSGKFufoMB zSA9d@YMh3b+0Saix3V~+ZOhGocSfo6zRdO!$4AYeFBIXF{9KH$mS8Cj2feGodZk zo1otj6Mo+^Oq{AS`mDUq`S*0c7&fQgG2wUUt_eQ%*o5)(g$X|W(!`&MA*J81Cirv> zGcP56KDIV8qt6SfJx?$*o)6la@q9ACjJ6(Z#_vU>8TE`e!#5Jl`2AUC##pn;jPZVr z8TH&@##p%5%oD`pK>9H=+U1lye%p+CJ}_e}dtt^n|Jn>+`O}QC?h6aX_b)B5Qw0m` zVY0yfwJo?`eGA%;Tk!kY!h+xBwievCy#?*r!Gh79v0Z4zXiW%V=QPN ze|bF6g1!@JL7$JYz}|vt?K8uIF)-mh)t+lsSWy2}7WC~D3;OUz3#Z4iQ(G+PLp!A2 zhb`zkf@;?vvEcdalm&hCg4FA}1%2{{1%3981^)fL)a$WD`hmRPYYW=@4_S{~u4uPB zuBcC5SB#DMT;+43D}GA~yW%;nge&Y@+La#^k0b6buCSw^+V5~fSL(a<{?d27KHPV; zXSMN*@a*s}{nrlk?EbGG{IC0g=(|nr`)<~WLq6PjwHK9$aHgpC!#g^dX)mlPu}pih ziHp?FtFrbYIVqHBFT_DJnf8J(c{$Tw@CHs}+7HsAXnA}#Q~c_B?V)bYO361@m!rKv z-AiQ9V}?%Eb=r&5)`d)afy{H5X)ns9)0sAa{hrLUhtQ?-STXHMVu7wV?jNt82Y;H> zE1qdDXztS)^qQlWl^*mSEM6L=+%-%a_!h4ydB0^$dno*22h$$%KRLp* z7vRT9OnV8ao2K6f^*X{Pi1)9W*UvNUVeaJyX~#sJYA@2~x9aD$7p|%M^m=M9?lbo@ zO-%VaOndq8+avkU>UL1=tGy7`-y`k2LDw7SH!$smu;wm?AF$0#d)XPkn`tjgn>Xrq zPMvuC9Y$9%Lmo$NY|I6Afc`Nk7}9w+kL@_A<2ZerfMi*OTHz#|7jBn8T<5lseYq9eD}P}v5Tzr*eK;6VA`z9%Jocp1Fw2mIz_5vpBUY=B1xy+8bYsy$o;u2V^Il!L*r>)D3b_TcQtcu)`9m-)5%Gf?YkT zcXHKvv{_J(BeI^mbm!J)8M~a5dhcZPa<}05VJ6;XMAPhJ+ALZ~x^7oBn4&#bF(n`R zA7k2+W{$mbFuEYy@qq3ps`F#8$&iEU0lo8Uvza>&$iaQR92`zc`yFI7$#?aa-ArAN zADJz3FucvQH~yFt41-Yy!=QDLiO$>8zVl}D9PRPppI&X=ri~{*xAXfVC7+MuwHoZ| zJgb8{?_c@&Cx^TEVBVv*8{SDf4aB?0v;q9=nEp9W?H|g&Z8~hsi&BGm(gn-&Iv;iB zUgc;1IqkQ>yyYJMl6B6xaZ5oVu}~+zuV{rXAAR4OFDSTl*WOQu@@y0Hmao*WFK@NL zT4CRf_IS7L)e-YElUm|_Bi#5OLw>sJ=H|u|uE#}$9czbox)B3-g}WgOPhM$&cM0bX zyg~aj)B9}h%Ol@hy%^M?JLVvN7{EtVeLMU5#R2^2J(hfPUw^!#uj#?3c#kWdFQXCP z_FKf8tDgP%!hWK{&)f1Jw_S<;-L*63HA=Yg&J`1*Qb%><^~$yO4TOFT9;fF3`6g?F?mjqwf^(}|lVR4-G?uN8MQb#tx~)fny4 zs}_zAZp&McAukC#RIWGU_?~*a?Xo=0zw-W+?_H7Hp2u+>vAoAIQ})Wd*qxq5x7==w z`M+W=u+J|o_=DlXK>1qmX^B6(++W*-hxuh}{KmT$+G{`O)Jy@PP5AG^Oc#^NVy<~; zC4Q<{rz^$9@fIGx6nHqbDep3?%oOh|&G3%bwhr%8y79|FBkOZ=kFa0r@%6VK&fHK! z%zboN7eAzBRldUYWc~J?>hR}B_cvLSyCGlq6*Xv8KJR|XhKy5C(Dd|FydibgMU8r zd~Wx$75SQBH^xPeHuFja@@`r-ryAxMrd5Ui4X?(VU3&4R^-32$H&^k6{pL9H`R2uG zwI`Z6)n?VYD!ld=>Ag!gw(vE-_8E1sUs+ylh-iu0Uvann=c*35T$%Uiy?)P=el_Wx zcc6Xm%`W1~|EBxmm;Jf?h5swXkJ}6RBlkk*vp=rNe<#oWyf&@e=lNgcD?WRz9a3sZ?5+Jhk6K@h^K>J zVDUY8n4#~%CE>+R|CE;aFJAXg@AJRKU;KZ#K3^J?Z*SLss+axs|6BI^zj7a=c6Oxi z-(CNA_Wir_j_xZu`FQ(Io_**@Z$|+i@{=%mRg>cN8x!Su%4u}0&+nYd&for~(Znx* z(P(<$eU09Vxvf#2a1+NnUDT-c*aaN-2DQ%o4#y9l2c36Dqv>Brx-vs6JFPi}^N~kU z&V2~wX$L{~9{}9|>XayHtfY7L<9umQr=-0oeG{Pdhf7*d(nCo&K2g$+k~Wp} z3rQa&g7285%_Tjt5XUVGK(kBw$$XSA%>&Jlw1lLgb8%eJ>~kc40%(e)MI_xi8^?P} zI(8PygC%V%>8+VK&LvHWM>$#2K9aVTG%XJ2+e%tm(w;MLJZ(DYd`YKC>Mm(lNwZ68 zod*7vk|s<+IZ{$SnI5r2?YnPQf#S?OT1$4`-a}9JUmH1`9h~txWNPe4b|84~h=2m# z%(|e(wX(Gbv)eo7&8QLR&vL|!I5hROFB>x=snCH30j&IhB9)f$u|Dis z+g~#xhsCl>=ITvDhD>CGSnCxHw)ru?pT1vO;P=UF;-l7|ebpzJ&91Vw>fIU>m|OU= znO}|xVjG=8e|k7`8uRqaeLJQ7RQCB$ubD27;#q~S)0!s@_Ge8$d()|(w?FH)S&DSc~@W{yg@jmF2&``AzplVQl^4Dl^Y7NMvQo1VLGdO zt>Zk8lCf;$ufJWkJoabr`rb_GQDrKtbnnLO`%`1s%FZVruSuQ4T<&e@Iy`a<8yohz z&iLHZSe*;*F(p2VVdvGpw?07EU z^44$Xu~s25-nC9FV~6Ig9L~AcWJCwn-+J^ zi{qc>U)gx{{RHN-pnT~0{t0Z~wq2hzDzucf__;-&>8lsBm%XO-DD-YI8(Hk>FY_L+ zVb!8SUziSTU`uOcj-T?vnl@31eC?I~7xwAiFjow!KESatBUcoCM`_Ka^c z&NDoIpnEm}-`1irMPTM%`*m9pl3COaM6HEk+CZp3$@edW`u`$?nnLdTi`_-+>Ti?y zn{*{gTNklg#O^9~dtM56CSDaHO}xaH886*~c!|=cc&ThJSLs0arOi|9USf9>yNB3` zzn0jQj^a;!bHuA;>5SqJN4&LdndIM&cvanPFUgFT>`OY^Z<+CuOrvtiezcqI1ZO5! z@!H$sL%h_Mw%dnz?d#`=m)gy2D^&dbAzq{Ul1w8Z_{v`P?d(Vgy1&u75wB6XWLwh3 zc5|cy@fwv&x*D|u@fx)+@fx)Q@sZElZjS0k^)kv!d`4|XeaMKuA69jpUJSjmo9=B|U66Pf@V9gQ}O2wotrAY^HdP)QxzJ z_BY~HeZxrIl)oF{CEFUc1Llyke?ZqOT0$< zJk^c(jEpN)!T0~j?aD*YQ3j+LMMhHTrr+ee`{VfSunk=DbSp&E~MQ8YUBc!-{5N}S{w)p`F zM;K%U^vlf(n7|r;gav#mp>8Jv5wiS|gmhI12VcdBTnqXr9)u9_3LyvupBc)vppYM} zfX}=NfnK`<5SlF@ZrXP$0Zr4s7j3Wte(^F4*WHj1yyjMfXio}82>Vh$d=B3q1~^Go zAFf3No(M*$c=Z6-VT=#gzDr#V0wiFh72yZ-!x0V;8H)M|J+|0;@W-yffbp&fQhI7zYw@34nvxjTQ*BegJ}c0mB)Wk~Q`IX`H6+z-H! zp`i$HxfhDy|Jz{znyd|h$IP+9lViQm!o5NOWZD%9AW(ZxF1nAC(*SP--&Y+Ei}nwL z2WF3gC3=P+kVNzhjlX>Wg6~6vq2G65pig}PmU})CAcj8UasO`N04AOHgPcZQI9?=# zYhgjoz5q;J42Os2jzSOG9f81)O_SgeF%tmlIN%Qu#-UJv586dY(4{v=IM{fAIO_PK zhy3V;K=VbyKH@%_J+DRrq`(5v&i8#W0;Gh&UwT6`=Cc^aM}D-WJpal zz=D3L360;?gaxWOW9-T00*$-4;Ep3*0IS&Og0W$b3oL%ZMMBVApz#?Oz=Cs|U~vIA z(%``bP4M{Q66VvvA3of^){OU59VgHNOv>F zzOiP&O2(N1dzxm(Z^j(6glEa?&zUhcT`*%TzG%i6bXmfNZko9kc=*PQ(6nF8XopG` z$aS^=ekx!?+89{Zf*x7V0y_vekQSXS;6PgR?jQ-D@wA{POtfHZ3$~!=g{)b*=+&rN>FV?+iwBf>x=~+epbSmu39pK(YQDLphr@0TSxCj+L-b$ z{nrlv-Tok+BkjQOLCby&!RD!Ykh%(H2O?sz9;}S`pS?_jfaW?cgQ$1NxTsC?`c*O*dzBv4 z4#dMc9eJP-8VD|*tw;8&AY=rOFSFBkJJX<;+mA|p6LmW%M2kXrpx%r1`)lE=Wq0Vo z(OL|2*5f*YL?IcJo+?;e$x(<8g_J;$_ck3tf#BvertYHUM{{ z)0qYn?stW$U}+V6uk=&$RS>z-7s2_P^dNB{Ahzj2=s={LminHf&r0H}(tjWx+j!1cPP zE6ffe*3Ri30zKFN#S;|_0{7b~?~|anBw8X-_Y4I`zntGHn5aHixB0FO#r_lRErC*;9FtwkV2Ah7J!veV#~<->8S# zYj2Pfk4ioE=@=;GF~DlA(J@mR;N5k%j_E;rWk}E5u4BD$d@DmFM~*IM1E^` z{`uxSav72)1Gyi?lFysODPsPQQoVWbox=$`y1Mg44XT}=YaY)JG!a41?K<*+l6Agq zI&V0C*s{i-W!4SmXPrds`hWr4c}K+pc_ZBMTvclvf7Yu>qX#<%Z~++DzSX5Or-=Kw zc|&cY0ioGF zWC(c4PW<2>8+>2YX@ijQZXJ28)s1@e?$(=&!0CkzwhrLW9;eqWcfT8NA|mR4_+$Wo ze5ckI=^NWa&hy56xsM2rztfseoWHo-HowlimWVt*v+Wza^Umso=jslf09lyb6K`U_ zHs)5sFeKFDqnBja`B9-S_|N+WuNu&`G4Jxxk{R`~RKh#&$1Qkaise7pi02j&;FV*V z@^d$C<&DeTTt0Un5 zWZq57s^K|5rHSM>@t(&_2i6BP#&hs>0*WnW< z32=dciE99!dRGr5EQb@tH#+2lY*H=9@za z3t-ADG~IpV2R%}Ir*(QCZLXrR?Ws{d6-{k_y(&98ua2wxsOV!A6|K%I>PVhC?}$3O zP90ZK*(zFE9k-XKqQmWvt0-)BUPT|P<0@KP-A6?^tFnrAR`pg<&no&_MLR1wD(YFu zx2Huiuk=*MRayLM+ecfwkI&pA{V^Z;DO|wmpVH6% zhpzjt%Kh)szlf~Y-|3o zQ3 zCEaxd$9K!<;u)Y$2W6D;AW)}M2f#N-(x#x+JhEJCKhBo`wPx9e@?#l2JY3QoGJ1H| z4xDeZ4YazX{+m%gya{y3M$i@OLBl2Wmb9a!Z&Go7fQ%a6E~AFCuEO#5GAcMlMgvDn zI#|+5%fR>ZQqYejy}JbESV=D|MtPQ`Llu?LzbY#DXBqYDE~9>XN_uMnuJ1V?)FNr0 zd6`AKx=;C_M^3a<-`V46KbB`&i95GP_GfP&bi00MKo2&fo7J@E-C%aHQztL~CIi{4 z*SDIvATY+!cr+>ElwaT|sSw6SC zA9rlz$tr$%tn8)qiOkv6e^}S1foy+zocpbtUhH+lk=JJK9>pqktTttJLI7*I?nv9b z0RWb$!uwC%I2pZO=DYX^t$uiCyUvvDu;&ec{G{Dzj5=qd};>k6_;h+jh<6k zlI38&b~}PuSVFs;gBAs_eW$*h+RiP4^}Sf*VV4_$Y^VS6jbE*q#sYIySlA(UKI>7d z&X1vaMG^mzP| zM0Wq@G3hm*B(SXwYh>G0e>N-sx?}rwd6QVdN!N)Ej;eO|vAK9ALI)2PRVeyf>B zu~h}KHC@0K`fPf-ss2K?`$3+!X)ROO-j%iPw;z$j9(m8L9(FB}l|7L+uFA7yRxQr- zXNvW@CzSF?pL`~AB6>y0eE-Ggm8qc*ZhBuj?^`+6)%0ZroZVJ$t6mg;w63Tc}Xwg zrO(Z@5ijv&#!L4mUV9zXed)ORJ2M?fCh^$opm?c1w0Vl%OYCl9_YgbrlU%Ze{pN_5 z`0cmMc&XlIdqz=q#H-rM_L9tash(6XBjG@bSJhw1b|jbVNV09W5AoXD!Vxdo*Qmb4 zYgAw2C7o@z59wf3U*a>WFSWDsaivFQI;b`?D<(-l#B0P3O0JQ*5ijjVbs%0N_EmC? z>Ol7=-?ZHv*@1YC+JShD>OeAy-*$7PgVBCL^2l$D)J^f)k3E^`K)fW=h#eHK5jzks zwVhG9#B0RYNiNB>AJaU=F)y*(j=#i9HXzxynHbD+M!Y23 z$oQamjp(a*jo6HMX*bdq#Anoh6|Ygbq^ps3Be~SZWH;N*QCmSu)4Xulwt#7oq6bEJb&UnhCQOYLF1WhR&I zPwil&EflYjzNu_rR4%ozQGJQmh|QD^M*2MQlHVBBmw1ipOT6T3wwt52FxoFjU*a{| z7Q|~*U*aYCw%dpFH8M_;JZdwe_9b4jo9*UE2jV52ZMP5c8r7F{HmWc28uc^cGb)#O zjn<8LjrI$&ff2bRllq&zzT%^Y5m*3C-O1Gyv`i}X!VvPfGYY}T{lXB0+b{xQ-;ewd z!VnvYpx{@*2*zC(jL_6Gepqq;=U{}O77j(Yf7URB(&m;KOVYv-+&wDG$gT*J*H(@qvwiln%%7W8nzZSQ?7Z(t#mZv3`Wi5|TdxA-vmzVRs$^ zeRD)2*!trz)W5zj&Q}jXDD8nD*m1bb>~kUra+wuj{(EDfXZZ+(9E;$4ZE4H0GT6I# zB*OjQ1|yS!$P}Xm-Smq_i10}X57-_G7|Hz@wEvYD_{l5@YY^?Gg?W7!0+@nF2!g>` z3__thMnZn!NQ5Boj76yJlPIiQ|3u353I@#Jv@iU!ax~g`SSak15sP4V;a6I?S-WTi z`ELqF=&&gm?IyBxXyI#rL?J{sAqp$-i}u$R@J)|Is7v7(gzl~jK)b(-L^g|ck${Oj zllDsp$9+1+q8&>{qW_!=M*HRugrBsE!YY#^Vi43H9)^A+U?keY$$k;I{!I{l95Mo* zu<5y40N*!ZXz`no0GyNxhNV4%VX^U%2$p^&!{uK`A&5L(<|`4on6yAX_Xq$Bd8ALn}oBG<*AA;>jRBpe$@NXg$tO87BDFQH(EOuk~y+whQYH(`or_L z$q{a77=qe&$D*Cz_`^QdXjzXaT=#ttdSt5*j4Xv>Frxhw2ED(Hg~tyKfv3I*hrFTE zLz6-I5H3@}%zlg;z!}e&5>^0@cBd7*^YAWrpEd=elB?3LVV=#K|jwk@B zR?8e!!f&-zlEcEG=hbNV)5a+HqkjzA=VUbOA^b=SaQz_~azqX>ZFSzZ(HP0DMZ!OB z$-G70K?vI4BLQU}2cv!?L(vZ|1>$_wSb#Ae%G_q1xvW3_rskHLuXrxhbwRx3t? z+%On2EOja-z&ylV(DFc9}57?UC{hn9yPuOj0isJoc6eJy>M7(E>7mFv%=tCj16H zlh?hK@E8{}elr@F(SohaXu-B-jHLo@qWuQ+He+o5){HDcK{7kd6f;Jgxn_hJZkE|@ zwwlpXMOGUvKz5HABhg6-U&)Yg76D7qvf*7ZqsQJeOMN8_=DF14tr;Ui84G|}Rb`f% z8WxN^W|{4!mIXb8OE}8c7TmX;1!^-pa=}MNI;hb_FXAqGn*{nT5#}a3ztgMhnWz?~0MU zfGhm{Q&#{01=WJ-K6izF#a!`>@|CO1Cgh4|G>a>qVO(ADTin_ezX&5U&1#bw9LJ{T zQiyg@k44&V!N2rhJN)Yp|9|>}_{}u}Te9i73jCqqLkeu6z#{;6gvwkGRvm7n&MOKq z!+brLg#w2Fv@r88;E?|WZ7kAZLy8}IEY`su>OP>W^y3PAqToCVTmruLiEL%c&PT*p6zCz z(!n(f4u|$m*7b&+4=~6-r2Cn&w>IBUXg^blXW$1bBzWhdZg2R{c|8XUL6Fi z9huPQoXp*^MuKJb>+PptWZL{n?W0WP`cUvL1fri{8f3QDdpMrxxklh!Ha*uK)}<_` zj;pfj396?kSd@Z;sTNjz%F|U@jR?xq6|4)_FVW#uxUX2plQQ=_eh>Fjd@5^-(oeOd zN|~f!OzQp$$cU%%6>?-rlz_QyIt)!^u~9uv)mycI>iNopm4~aIkJL@yN%`qIEKGR{ zJnyjdjEg!9PU)rKWYGJ4>L+bsiLcK|`z+Oajy?3&QW$VP}9Oq{!>H z=p&7Sm1#3o0cUl+Q9dTGKd%EkF=A}gdyYyCrR)Lt*CD1nRl&qmdn>*1lzmFi`l7B^ z@HhY^&*=K9o`h^T+w|-_s%N5I4$FFPlnDh6$!s^5q&{bumXCMd0htZwte#azd%ACR zo&jEWNy6(6F+eqr%j=KIEJPVPP!qt-V|tby08+Oxj!Ty*L z7&s6aO1!$_{iRl0WOzB;ga0YeY^m=3x!@aV-S{zayX440@_o4;4}4c|?5U2;d9bGt zP-rm!vGq6MbryEw%_!x|wtD<^DUk!HMIBD)v-6$WAXCm)odD-6QkGW{326dLH|ARg z23~GhrWt08?p6dWFRm;9CGp#hk0x~H0$lj}+fUm9HrJ{t?3uSNr@SdIzpjKi?YGU5 zfu}$V%#<8z!M9D$@#3>v?U8BaaC030tRp9s-)}ja@bMkK{A_qbGhg3v#M}ufb@5*J zrwgyRFI?orX~_vKIQ;t-TzEqlubs`f$X6rqc{R~K&W(_nsCiw?f4r*=|6vWe*QQV3 z9bYIj_Oz%ExL$<{u-~&3{~!|Ev|QT=_3c(4Fv^e2 zfCakO;v#QSz`N~VaFI7?V&3c(_$b0XcCN}NR}!gkZZW*STdQ&5r4uI&DGmRPZNdpG z=rh5Dc6jKD1!-EA=2a<$jmV7CE;L_svh@x_(pKLGN0E?bg(+kM&W`s zJy*fu#LTSN9K#0f&#SUMbtF&msq^+!oyTwTGd)?3J)bHoepOcRJO$rV$L;ylaYs~L zudY*NNBruzg6XNUq6+q>sN%ENSHbzzaeEk_(nFnBWhGZp1>aL;#cxmT^-%Cg`{Syt zJQLDt>!?)Ny-QqB^c%kP0@b^iovWMNtJ`RM)BFN^f;s*;gG`*C~Fbr_xLD zspE=I>80f1neTn_9DBJ+o;t4h6~!~iL;B4%0;{B+D`1t8*_@opz(xNtE&Vswsq>23 z-^bB$dp=bzW$=EA|DRH&&)?nWpUSo8H+r9c%BT1h{kwWNI`8PXl51Zd`?5WsBYA&! zo&EW;PT16O`?5Oz-%ai7^LO(9DZi2@dae=psm*UL+2cr-uNxB$*iuEK?s3^QdTDw# zjrKeO7;0k7b54ATB_C-tJ@mdtOV_^*y5JUQc2H}cn_Af^+jWhm54eJI4p1I`Nh@1_ zyomFC&Ve33i{n!y47CoZb=67jxKmk4e>{Qma!_l8q=6EKy8S5R_=7qXk<{-nt}7<# zS_x0>0&2}K>0Q87ty3g@egN0qPY3;4(gG5Ox??wvFWd!MO46%4Q0^z;sIFU4uCoQy zX%p!6G|)In+eq4FBaWAobi)ReM@hPEJ<4Mw{X){EsW?6(1vFeHwP_^j2?=8jmUPxi z@YR!4!C0LooK?{bne--A(nLwCNP1-%v( z`x3bmZQ+tyC7mx}unHD?WC8fxCH+Lw6bXw>orm*Il4i_BdHx(wXGtH(x6>F%IPc(~ZJ(ZF=tG{huz(Hg709R(hD{gkwWk z-W%PkOlUcp)o|-}rbAEL1$uPj~j@;{H7IrU$bI%Z@PjPGi}-o_C9VzI_N=(J3&nYJ51eo*8mv!ma^q>B}Gb zZ+98R29-J#`my(LHszFgQraXddwz24{XUm{S<3qRg=^)EV%4j64N3TF2CFutW}6KW zVeI>)@$U+^v$CTtc*mxNrn1I+M=zb0*OyfdyO*)BWIQ{%waJMrzs9klcbA+gaA7(# zeU`CsOKb%D%(OUX$v%Fpd(gy=oxk^FWxcDMnmf$Of~zODnXzRCE7N~K;ckAx?83>r zuQsOpvM*6WMUf(LJYg2eDr=amc3|4oZwZ-GBli6?k-?&XFA=hQa^zu zAD&h3qtTmK)w`|R^t}_wZm(;V^|Ng;YK-1%4~M{O33hDeNvgvuC_Cxe@tXULfyy2Hcn;_i&wr8^UYG$rf=$;nn6h{ z&#+h9*L|PFcGcfLfBJ#NEZH~j<_aaZup`$K<|lPo#+>7tR)29Lk+oTM?fwz(Ijmc* zyuUw)j=x1vZ|-*eE7 z=*C;vnjdxyX|`(vTk@`A+tdPyY;38${@H(8&kkJuwA$;sdszOnXTI$?cOKhtzty%q z+g7uRXHuT_>ARe{b{^JY(YQS<=$FOeZSSS9vZS zR$$Jg;)#ovv)(<&Y`$M}1A7#>GSArSE7*a2b9eXIo5C6mYVdmaPb*mUX{A>l>a(7W zy;h!0Y`lye88CTX->XYmWI)O~mq8+PPQ7P)UVXQbx%JK1;XG}fV>q*sOf=&Sz?oYX z`8N|y4JU-Cb`jJpc4x8o7rTqtYuo%ye6_?*o8l!NN4#{7HsYm_XbPcLycF(Af73CN zOO!Un>ni@X=Ow*}x31V}Q@lzR-B10UnOu^YnOwS#>d{nmALCK~1(B>(2FR{Cc z-9zlePc~CJia+%&GhWiytaH5o!x1mZw%;=2CEaN!o7rw3;Pvje$Hkx8 z=15=SRee&jX6oxCk9gI%@F8B3M{Q=rW{TH{&4`!uFru&GHL8QsLx|V6%xp${X0aRb zN#Zr?XT)pNcNDJ?-yvS&BWk-js+&#7j0aqA&4~-~1C^BXuKrq^l91RJ^3O?dB;8UShZXeN((dZ8t}}ir*-&Q5{IG z5nm@hYGb2viPwlvl3b(yK)lpuL~XasY)1E{HZ#(0iq}Zph}Y;isd$OnZjSV&{$o^M z;x*FW=-yN>BfkN}OSU!AHx;ju-z?%)dWb*u&5<1xzs@NB_zY z0iIMR41mQh5rA~Jjt0cIWDxAJAO!HTCPAoQwlDxPhlB&*S4Kj6Uxng){wUZ_K!UXl z7E2?fyeP;I3*lO`TEhVJ`z;)+zn2OH@bEzd;7~JV{aVS4dAB2x?WkcWuD>qpiqTyRYY>iJV7V1>4s00cXSga(UDY1Jvk$4)CC!@ozuBeuuD!vkXgRjwY17X32@mMjwtpymZDJWljL z4XAipdQ8DkwAA`3fTRb-!qaz1sQS>!@PtMoSTtgGAbL)q?4eGwB}WCJUL&PF-o*mq zyf*?sRiVGu(kujeh!H^x?JjGDM@*lL<0~V;*D?q+OE@6vVv&bYl#pRw80xz=6eCu5 zS^rLQ1Qd%+Xrb;!g3!{ZC1iMvtoMot05(&C5b$Y{nd?G=;b9%E=(!DL&vlYLxlbWIQ;IwWGp@r9t1tb;uu=J z6N2#z5fp*;t}z+yw>22`YhoDsZPUXMcvuBJtVcEccD}BOG4hX^Xwl-%(CklV^z>#f zphH~H?9;BT8WtG7`2m0TL`Si5zYt07F(2 zJaDqi>?f#}V=-PrxkW}n1@)$+uM=0A0D|9P!U%KBgfaEFyx%FAnJ_~_!_S*gzZ)iu z)wfLu4Zr_B!=M&W`C3B6-X^|ITtc+lnGs^&PeQl*n*k*rY{p3C zZN}mlBIBTz6KRnd9iE&qdALYd77-F{)_8tmzn zgqmluVC*kxK|KWoTmu1#422rFMnJ#;-jL92lSM+5E%5x7GUH$?3)-`d1tY(JdTZIm zdPzvNn*|;`LPE?1)$(wRkq~o{(NG(S1=S$ylip`E)WAm5C10F`cF&Oz^d*0Vero~6 zJ7s3XeUeXPIMl!>M=gL-pLkET%r+-wCdP|0G7Kdza)Zcu)2?2VW+(BA>$fYrI0Ir;flP6ifj2)O1festn7;R zb9RLvG;&2xAMA<|e2gpfi8c%ox9Q;rOQ(7~7`I3pgZ`!e+TmY+_`lsB#Cy0AD7j4! zSHMjLtW>~E1&LOGQUx_vWtCA<@dHG#QHTC2U@7=lNq}jR4q;XRMFk{PkZYBBQ2{6c z^jxC@`4j|PLD8YdVO_q`O94Fa2t; zFOv}HcpciU@&-aLE7P(U-@GURohb|%Bri(9YN`Yz9+Z&uEAqY~TSovvFOM+oS!Td@ z`fLeIOUUzn9Uu$4Y|+~X_1PzZrrR0vEFRSJGphDL zJMPt?^$Pl~+6fstkLvdUAa);9dMUswpxE1WNVjT7_|G2Zq8Vt1gl=D83it{TNty(L zrpkPWha_P1g3Ob6K!?g}*^rl{O2B9u11RMjOBhWc`@K347Xa9F*-ocrrqyeDK1yvV zwFRg3B?FY6B{f5yluI8ZF$Gx%WIRLq(P7!%>-0R0YN-R*VU-SeReq}SOd=!dEvBIK z09hT@{Ymt2Bam~O9*!1^l`S2u=a^LxYxK})=@}w0RKp>?_iVIOumt-r*7@KW({;Ik zM@`Xth+515@{*<8P#q$z&j08A!}Xr1a?L7;GhjqXdTw6G4P^?-jIx*BGD<(?!SK|T zdOhqR^2io-tr`}-l6?^UL9=BmXS6jWdJG&SOB`MX2XWrovT4E5VCN3>nK9n^glBwj62 zfe|O2IcZzEzC-{Z`)T@O3JM~QjH3rwOYKH`^^B`(aSN4+7VU6Zju?ruhwPLi*d86) ztsv>D-f9GfpIy*1zQVuHOTReF62GA39QNv&P}L$4$fbNj`r!e6Bvm7w8X<6>I}8g| z?brJQ+Wn#&LC#Ak`vrMES%<=_5l8tw+GoEGQOD8^d-XE>;gFu|S@iIp_Bn$y$uZtHcImO~~so8#rt#rLV@A%v2o)>G6#hl=bvAn^q&a=vt9n4D}y~mP&8UyI` zPwsr4K04$Ilyx>&6e}x8@gIy32DQ-##MuFNuXaqnu>oU1pV;;=RBP5c6{b zc+o9qoNjDrh8ev%ZhY-#fhlj-p7Z3JyYpYG%JUJ1{*xY+MSyP0($h4 zdVpR#b>QwBA9-%t)sBDLbl8{|L#lEz=ia8B`9^`{?|8f||Dw~C;(mK8!ER&QU^v z%8ev+yfx+l>Txde>rT#*(URXSB3gP+W8PaN-o1Q+qrFFbjZC$}Tl2nzA+O5$;%b7ji_6jCq|{-e%P%gUn`oS!zr=E8$rnrHV2fxQZ3A@ykorj*Cig`>0T2{ z#8t1zYmtKmRR)y)V`sGgpt_j*E@9?3;shqVoP~d6DJWQL!XLkEAR+gRzw0Bw?7{Up zA~$ zRCoWio3pxdLf}8UT!&Bb9#=eHD(7k6{J8DP7xmyb>s>MLbgd5fS~cWbX;Ftr3uirg zuiq0}3H|O=O&(N%)&i)^&x@55?o4#WjO7x;8E5jbspG1ws5-9r?V-@>e(Jn}IIH6d zDy^W`O0J>`Dz5IM&AbhLCf_0+$={A=GOPD9`5W?tDfx=3>y$iIR{APFN9WaXMOA&2 z9971%!6Uhh%QKxS2)MdlU5DSe$GZNCPuWHB;WzM!UO#o6qU!$2jtYvd%1WNnPeIdF zS=m8RrH7KEu2=1+_*7Zxr{vn7SA6QcDl4k)qo~SsEPlg{K;5Z_3#dDN53dN<|1Z-D zcGugNzp^{8?q|=Zj{hAkFXgH09nrs&qpnkYs_aN##b-~A@;j2Rj@!$1R95#>*Oim@ zus^QKePY^7w>(w9|EgU3>;6t}`}2Qif6>E@ zK;rFw!ylYQM&_2qs%jAUTVa62JC_1=strhdM$Q+U9OBlE2O70jxv$a08@DuCIyb0w z5}@$bQ#Z8ZPL-}}lxJK5tqPjH<08tP&w+3JS*O+28Ltu^-H_=0qlBP1;b z>a<}u&R^XLx>eHJl1|%>a+&s`6-|eBqYA2r1>SiA))a*CB44_{27uaOS*VDj(3#w&@zebBmw=X#G^?aD7E3uvps|vkNJQCP(j1aLScKyZB>in6%8DM5<*f^FJXKOJ zNt;UARzl`gM(LgNB)_ER=c26W$FiJD(n=CapOyf=NJ)!GdTln2ACRk`&r?I?D#J=;E@r5js%{+pXau6wfBNB0KK@6wl5O~`#HprDlv%X&Kb z%=9rVcy&lZk$a2OTmUGW2cBuIU zm!z$}tW@V`2?>kFvxfb**2r=znhm>K^VF6;i&^uY#fO#-^<|T;+J;b7;e>^sF&J@-tYWv}p_qVc+ ziH+9BUYO5n&-&*4*lRP{{^D^x+07~JX~4)4z02%kl_u;dcePv;yLkO*@P_r_Z12Gl z-wv;{gw+Xta%4xyG?wtTns4#i(^-~>r}n*Fw}D-H+azvVj#&2M2hY+uk4$IZtZMv6 z!rC;}E-riQ;}tVl`#&#zn%5G;zIQuayinR=c6?IB*^&Da*zlPr%|#|nVMPNg<}CO$ ziIp7sph%A|XS1I!z1vf~@&dMgz)P?5t7o$tk2b70`b(mOwhMU^6 zPO0AKx_3`yX?x;-Z&GwKTitbgo%4m#SiL9LDwqFZ9UFDiqinw_$?S(}53R+Htzva& zG#IdAR%j4N<3H$Qesg-@oC9&~?Mid;s<1l+s?(F)m)mO14?mzFV zYFWZ6CXfHMf4>#%)i>)p@wv&YRt{e`_q7XImYigb(mgvT#l=0>|O4cWq$lTmCbD5sAgV^_=fPj&6Te&EobFS zmKPf*ejsD`cmo;3%Nl^#U;5}@h&7hGacx9(m~bDQ~d2Eb~mwmh@EsM8>#c+PknR5tM02aia#9jlFs&9X1ru0vpu6I zJL09f+Ha0{sh#b&%y`KL)DHIiAL6yIpCewgkW2nxyM2gP`K{u2#7k{qR4(x?r)=Sv4qg!ih(kyS0+YSGkt6C%D8l=O2O#|Mxj(`ddwC;Z)*}Gn zcx{3ZklW6R0Oa0Qgb~L3BOLH|e|*;16NqX4X(}W>9N~WXCnE6lsZ~x-2O!MwqaXzE z3dn;N3i7cJLiU$g5nfi?itEP)g9iBFbCS?oo1DuR06hozAUrND4EkIOLO5C(sl)`!pGhw>^v1km*iMM(F#Pz2>)2}9`c(hvlK23rxp*fInmz9YS$ zv5O}HC6|ZdgH#)DfFyE_M<8VTFqHR)p+%a6Ah%KWa0I3v^2POy#v?$mbvVkGCL+jw z=mZ2ZR`o();BlEZ?dM3ez`6(o^yT)&^&%UT_9?en1ai+zorDi2ExZxHd1gFW&`7o> z8w+X^)YmFQa$Gg|+V3@GNUbx1s7pA*M@zfFVqINiD5}iLw891O8L|4W_R)C1$?5B~8!OOv+3{--ARkOXLp9MD6Gw486RPs&iFH2HC7x1NOp9|u%u+-9l&UOm7>SkOpFkBtQp;xSQ%;)Ki3 zM$7e31r@TQLVr{!f}$$)M}<775D1msKm}i^>;JQ*Qe~*Z33>kudT^EsZ9&MxMm=Oi z1+C#T(r&5mVcim{#SnzHi6P6tIo$)$lJBzusGvg?#3n3gBxJ(If-95T^TO{Y@-V-Q zjh6-dPg|Y;O;;{Hz=s{}IFJ{16YJlbhhw#LuZ~!qeSb%uu~B>r_UOd>h{fXvwQYsa zhfcMz+WwW+SnYgtJN{sJ)wi=7b^zaRU3ut2@eSU)4Hs+vmu{TYgnMV5o2&THCit{n zxE$|5cf8sHzpYKnBQ*R(ZGP=ykY{s|Rj&BLewB}>*5!H2i;wNO>+zf-f}-Qd1_;f6 zS(~RnS=gXj*Sfs&@lOtK|G6GQKpHpX;)C|Uz}yYESbM)c-)OX+&w4j zc0ZSedhtN}ad6}Z$3g8WMy#f;sQT3Dh^kMcimFeS>XV)QadlpOdbBUA>l9Vz)h9ni z)p`4}I<7wDDXKmRs!w|A)1RUkr5-X>R#aHf$R|JYLGj6to(BuSbpKLm2Vp@YpHysm zFumf|qt3p+a4O$g@smaeTz#w2(pm3lH2v~*&=S`)>Xh#S%BRn3lxLOn z>{+dBHGw+S$N+yqxoYrHd?K;-I1Ii=_@rU=k~DoU_&0!BXUnpW{Djd%uI`&6KV{sL zpDv;$J+}kbhfA8c4dp?Sn&hXAsLeROZWHz5Q2Smy<%7LAz4N5oj{*iWm*lw7H_r59 zZ%&q;_Q&W*mT&RQb{*gNvKq&CWIY(|$vS-1`?q#qPGuF7x^}M-J%t@DAN$j%JtwnY zvczZ0(=eV*yYf}taz8F*^>e(Poxfll8xe9~S=}MincIcWf0#CC8oTFRcZnR zax|}ZA#KV^$E$9ZH83Su-N2OK&Yb^F3DOi7z5aI*yIJfsb)lx?=oMc5O?>KgUh&c^ zyE<>rt6t}cm*m^>(knb&N3ZcD*CKY>h?oAR<4UeQFP)=J$tBtLymT+ysXHZC$)fY> zZ>58+T+$&kxnvuXtEMGA#ou0HcN4pZ*hvqvk&-3;)Hg@Gs_r_Y_`?yevj5-lQr(Q| zK)gnEpn9uz{W~3qkL+bHJF~hGuiA~sC0=SD+wDVgjZTG`g+Q{0?dC`>@hU%1{Em3Z z2KMdoAzq_8kPVF5f%uHtf&7eYY`ZzqH#2|#5U){xAYF~tEi-@5Ob2QM@<}6Yp?J-< zMi75Hl1uh7DwlYTv@yw}-AKC;AKBKZTxvHX_Eq;c;tz_~sJ^7D{g~#+zQk+a-^3%9 z<7wZSNd?6vVtYSxM9&0b{?}bR?;}%n2U3%2lK_A!B`T+Ed;Znt%DE-{x%4+kA;IUcX!JQn6X&+LYvp9 z?}xdyHnq0}V=nHu2*8xP zTLDL|xXVyP1?yfxl~YwxzA<+Pcn0ijsTriK@0 zSG_#3bjLTrs7cWKeKmUWiUcW2DJw^wBs zxm3XWd1|di@5m)2esDp}p1N>tHsFN|zM1^$0%%7a6P8wBV!nkYtQ(uqBK{`%jmCr} z52l;s(iJ9n%s#mU#34EVe8PlpH_uG)i03AJQTWk>Wr(ty(U%LF(GrEs65?;hA}^ns zu_Q%tGg`2g8Dm2Sv#g(4dX8MuA>0f*Ofh3Ve}PHnA=M9UbNU$<$Q|AcNw>(E^_YpP~*6%-Y-zB6?I z8?IlY&vz?;A!f!xCFE6Qd7dS6*rv+)`3(|6yhES2$BfxgJx8trJSs@C0vxK@X#fh3 z%I^UybqKWr7-~zLRo|!2dSm|IszY(rJime*t666NDmUo!^qBYGAob7C^;eKw%t&5j z04$!7*{2Wdv)`g+jLh5Hw2b;vp{UA!tia>Ybczm7Z$q zrCLEL)t46a-2vYrw(5{ulp|&S?F^>!xZ}DRI;a~S`<>nr3WBP>Ux=15@;$(&W#~J= zH?=(hshjPCdHRYYIVG{ZHf0>Iva)c-fl5=5P1R!*ccHYjoxL~#(mFS4$Hi>_^|)`V z@doYBOz$(kCttd2Z_x@HyKs>kw|_#y2rNzXNnK~EXo_tY*YUk%(Xn^}wY=*48tD6D}`Fnf5Lm>EGZ0~_5 zjM6;-(aYW**>8QCoN@o&}GQMBIXeBI$cmbL(Dt_Qp|wcE_mh{VDjkm4zT$Iba8Y(Y}_9v4%-3*(+`^` zfXTzQEWzYqG;IC>oe#-7@8Da$V0@T8Fg|R{6gmx)huI6$ht7v>iGpoGg1H+e4x6Wd z&1=Ba!T7+KA!hyp6f?m23s8SX3KH}KMoSE^{pkLerjdJqF+r>{lx*~zylT#+F+UrRPECIij{h z)H!d1(1uX@?^X!^JCrVi(zTl*;-@x3Xsh)Q+6YQpLAPootbvHXf^Nk~f^NN-0o^*W zX)#28*CGflwh%%eSO6-`h?#FeikUBwhs|Cr$+9=zu2|&Hm}39s;+1(z^^5H550~l- zXy)3-ut{C)g|IwshOhm}$Rv+Xopv z+gbKJ>yXu)d wZO>hMu_Q8Tjr~>a4chOIY_V57(h;Vszs!E~?z226qBr86$I%7Nj!p)1wo|>E&}eINhfX7LS1|I?Abj|&&->5 z-@W(Ud*8cv-o5Gk*kv?hAii~einwZ0Wr(? zMH()(^Q6Gv?{f28d~WEap3g_B_}nNfU(Eyi{HXk)RUnhd9a=kE?iFO?Vx}rFeql`r zDTa%Bf;XGZ`6+;pu)KBQ-Ko(C;;|}~3WZ7)QdfL>p>I*@EK0dTf&xfsQEC8$N+69w zt5c&IZBC9-i>j2ELV-#kRH@ORYPAAI6%r^XN2NkFs18OX5USKFFttX3!9bW7V5zkl zWsa^L1SD8jEEP=NaYn5IRp>CCS|t$zd!1UT%z^nLL26|fDR3r#*swL(W>|?x)1CWIvzIfES`lw zD;Hg+;ORuVzp(S@0Uh=oJ&21gLOX6VgRsr+{|PSf5!!6IpO>e%B2kB$0uNl6l+jc= z*SP|LrF~Xg(asJ@&CHS|;Hb?2R1_m$BSYB~J7Hw*| zu=nHC$FCY0^KG$@#w42t{pn3h@|?lJIin}(uDy_8+IQ;jNZ^O$;H^1|Lgex=)3~JH zp5JjKEl52Q{4CqpFNiEu;pb1a22OpnAjXk=o5|dN>%r=^DZ!mFhs|5o4-J~jcKvR| zQ;DYAlPl*%k4z7O)kbi)NS^C=J6NCS>O2=a4LYGArSw8=qM)Sqwwk zA#%XeLArTL35)Fl-o1k0U`+r^j8~4jNB#H(822jEiLX z3Lf4qc8|?s;kk!G{BarSF7`V~M&_bOpP!*zprg!FBx7kx#*h#pWmbxvOHpnS;-uWg zEG_f6flNj+wgl2`mx0*;)$Z{7c)wMc&+C~>F@m7^Q6z}%vyl*aCmNT2Hle9tHO3`d z><~N!15$R2jkM7ePiy>E`fgEsJ5>zY!;_&Xo}%->mJgEj1ZWCK`d+~haZa#l8P6rm ztC&4+rp>|F1kJ#%Hjf=*cJ$k?@_zb$P8H(Oy@txMYMnyRv&ips$`}tTsOkp!&0JD} zqnKu8q|ZjVeF9MQDp(i44b~+E#FkTYO6TN>GV@-d9IT8hh4P%XSygox?2ZDkm@Hl?|CEX zz&Rt5yT2+aJ>E3`HEq)bPeKq;B>rXBOUDCou{r19UMfSUEV*fKFwrWw*Rf7 z?^A>Gm+9`v+tS~J5aimH?dc|D$KkIl3Wo%dC$~JlZQMVN!8iTrk<1j6DR8c$H7&y= z8yl>Bs99!0jz&HIQS4!(>?K9Q_q8&UX_V4=+qrbpqYZrzygzGX&{Ulk(h3uuu$hgRrS-51ZYui8uBSE73DYHsDiJ zYO%kroQS`m+<>0@;0mUy9*E7Tmf_ppsX;HO)36!7X2ZU(mg}$i&*G0%Xt9=5bdKrg{nkhajczQY)b7KBB z(M^W~0~xRa4&Az{5K}q`JjMFJp>KbGrL(}AL(5cSRuNqUZjc1#VqBz;C&g&8CL|z) z5(pBRE{D_U5JFiGGMKp$Q(1_1Bx9Gs7V$7b@L+BfDN=n9ROBMAJ{&ZDUZ}W7n4oj0 zh%fK`)?9_~sQ1A*{k`X^uv?~2#XljRM~~~i!5WRpSZ3K!yiZgm8htzkJ6iaTVe3Lx zUw^g%-@9lGwm$F)%4KiIjx}I-+>u7SDQOkPY{|e^J)^=-sJ5a@_&w;Rvl+Pj;BY)G zb125GYrz_7^MtLrnSR>L4eM$-@Yn7gdSF+p0q(zNwhkU0U0QYnZ|-mi>fW0_Kfkt4 zo!&qMq%*YoE9XP06Q&rYzO4SwrlfurWB&Q8uPuGK!YJSOdVy(`H9WpWF(fnLjnRQ) z|1ghU-f9fGSRzTV}( z`CzS){k}A@rcZ%!_7}@vXU2Ut$)kYh3@;tl9&I{4-d@UNlahYSQs9MB?@*q zIF^?8%4QeWeR2MSY+G2@sexmCZmoGd+@D^teBnq!Yt&Yq{%&)|!^REgXUW%;`NQL< z$HuMORX!;2PGe)%mFk&cxMTyVoiaCHfxkip`mI0u>W9{tAoe@4vJYIf|x1RYCOUqYP&M&qw8r^No7=pod3bD{Vy*cu(Yyo**ep~GySVg;LS5&P@+RHEaJp?68JiQU%##%tHNN}g$Fyd?DBS;R3s(Cf z(rSHu9zSD)e8qm<0MNNYL8k+Yi26jQ^i+vyIw!G0+6xZX#A&wSo&HT(K!5TVo}^z#0qU4 zu{h}h(a0l2b=4T6E^0Ng`gd;;Ylf8*t4@t3melVcYDc9J%Qo7Hm5MB)q3%aw#WTAI z-n55UTxTL`KN(6aymFOj{Bt(3Vx6C;B@?94zf}*koeyflIYW@`}|6i_J;;pbk3nIy!Y=0C_g> zZqugI!94KEO5tjtCzp%@!9sb0aO^G)9lJ{~QQUSRF)>Ycdq6DGDFZyHi}>vnF3VQQ zwK6Zo$b?!-w+fV~fmewJ9^6e2ki}tf%-q0<*;b(>0mTe=JI_!C-+ejmiL9Taj5uSzhKs3J9d9o^rmG5a3VOm7t1qZ7ShQKY2rEC=VtVxJp`hz2Jt#>V33{<&lhh!!7Odoh7V4EZ2I)Mnx`fk zrH{V4`EYPBH!?gP)tb>#osbhAH*MOzdi)Tju|eH2Jeti!z2y%s zzMvZy@8h<3NrrU6kC)f`LD&)b@$QmqeGZ>+@q+(B5U*d_kEE0MNv$+4wGzCiR>F^L YBk82^NH`d8Nmt|f#owEr<*o960VeOH%>V!Z literal 0 HcmV?d00001 diff --git a/flystar/tests/test_all_detected.fits b/flystar/tests/test_data/test_all_detected.fits similarity index 100% rename from flystar/tests/test_all_detected.fits rename to flystar/tests/test_data/test_all_detected.fits diff --git a/flystar/tests/test_catalog.fits b/flystar/tests/test_data/test_catalog.fits similarity index 100% rename from flystar/tests/test_catalog.fits rename to flystar/tests/test_data/test_catalog.fits