Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
*.png
*random*.fits
*.DS_Store
*.pyc
.ipynb_checkpoints/

Expand Down
90 changes: 58 additions & 32 deletions flystar/align.py
Original file line number Diff line number Diff line change
Expand Up @@ -1127,15 +1127,15 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
else:
idx_good = np.arange(0, len(ref_table), 1)
idx_ref = np.where(ref_table['use_in_trans'] == True)

# Initialize output arrays
x_trans_arr = np.ones((len(ref_table['x']), n_boot, n_epochs)) * -999
y_trans_arr = np.ones((len(ref_table['x']), n_boot, n_epochs)) * -999
m_trans_arr = np.ones((len(ref_table['x']), n_boot, n_epochs)) * -999
xe_trans_arr = np.ones((len(ref_table['x']), n_boot, n_epochs)) * -999
ye_trans_arr = np.ones((len(ref_table['x']), n_boot, n_epochs)) * -999
me_trans_arr = np.ones((len(ref_table['x']), n_boot, n_epochs)) * -999

# Initialize sums for output
x_boot_sum = np.zeros((len(ref_table['x']), n_epochs))
x2_boot_sum = np.zeros((len(ref_table['x']), n_epochs))
y_boot_sum = np.zeros((len(ref_table['x']), n_epochs))
y2_boot_sum = np.zeros((len(ref_table['x']), n_epochs))
m_boot_sum = np.zeros((len(ref_table['x']), n_epochs))
m2_boot_sum = np.zeros((len(ref_table['x']), n_epochs))
# Set up motion model parameters
motion_model_list = ['Fixed', self.default_motion_model]
if 'motion_model_used' in ref_table.keys():
Expand All @@ -1144,9 +1144,13 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
motion_model_list += ref_table['motion_model_input'].tolist()
motion_col_list = motion_model.get_list_motion_model_param_names(np.unique(motion_model_list).tolist(), with_errors=False, with_fixed=False)
if calc_vel_in_bootstrap:
motion_data = {}
motion_boot_sum = {}
motion2_boot_sum = {}
for col in motion_col_list:
motion_data[col] = np.ones((len(ref_table['x']), n_boot)) * -999
motion_boot_sum[col] = np.zeros((len(ref_table['x'])))
motion2_boot_sum[col] = np.zeros((len(ref_table['x'])))
motion_boot_min_epochs = np.max([self.motion_model_dict[mod].n_pts_req
for mod in np.unique(motion_model_list)])

### IF MEMORY PROBLEMS HERE:
### DEFINE MEAN, STD VARIABLES AND BUILD THEM RATHER THAN SAVING FULL ARRAY
Expand All @@ -1158,6 +1162,15 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
# reference stars. Use a loop for each epoch here, so we
# can handle case where different reference stars are used
# in different epochs

# Initialize data arrays
x_trans_arr = np.ones((len(ref_table['x']), n_epochs)) * -999
y_trans_arr = np.ones((len(ref_table['x']), n_epochs)) * -999
m_trans_arr = np.ones((len(ref_table['x']), n_epochs)) * -999
xe_trans_arr = np.ones((len(ref_table['x']), n_epochs)) * -999
ye_trans_arr = np.ones((len(ref_table['x']), n_epochs)) * -999
me_trans_arr = np.ones((len(ref_table['x']), n_epochs)) * -999

for jj in range(n_epochs):
# Extract bootstrap sample of matched reference stars
good = np.where(~np.isnan(ref_table['x_orig'][idx_ref][:,jj]))
Expand Down Expand Up @@ -1236,13 +1249,21 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
starlist_T.transform_xy(trans)

# Add output to pos arrays
x_trans_arr[:,ii,jj] = starlist_T['x']
y_trans_arr[:,ii,jj] = starlist_T['y']
m_trans_arr[:,ii,jj] = starlist_T['m']
xe_trans_arr[:,ii,jj] = starlist_T['xe']
ye_trans_arr[:,ii,jj] = starlist_T['ye']
me_trans_arr[:,ii,jj] = starlist_T['me']

x_trans_arr[:,jj] = starlist_T['x']
y_trans_arr[:,jj] = starlist_T['y']
m_trans_arr[:,jj] = starlist_T['m']
xe_trans_arr[:,jj] = starlist_T['xe']
ye_trans_arr[:,jj] = starlist_T['ye']
me_trans_arr[:,jj] = starlist_T['me']

x_boot_sum += x_trans_arr
x2_boot_sum += x_trans_arr**2
y_boot_sum += y_trans_arr
y2_boot_sum += y_trans_arr**2
if self.mag_trans:
m_boot_sum += m_trans_arr
m2_boot_sum += m_trans_arr**2

t2 = time.time()
#print('=================================================')
#print('Time to do {0} epochs: {1}s'.format(n_epochs, t2-t1))
Expand All @@ -1252,18 +1273,18 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
# for each star, if desired. Draw a full-sample bootstrap over the epochs
# for each star, and then run it through the startable fit_velocities machinery
if calc_vel_in_bootstrap:
# TODO: consider confirming we reach some threshold of unique time values here?
# TODO: Like, grab n_pts needed for the default motion model maybe
boot_idx = np.random.choice(np.arange(0, n_epochs, 1), size=n_epochs)
while len(np.unique(boot_idx)) < motion_boot_min_epochs:
boot_idx = np.random.choice(np.arange(0, n_epochs, 1), size=n_epochs)
t_boot = t_arr[boot_idx]

star_table = StarTable(name=ref_table['name'],
x=x_trans_arr[:,ii,boot_idx],
y=y_trans_arr[:,ii,boot_idx],
m=m_trans_arr[:,ii,boot_idx],
xe=xe_trans_arr[:,ii,boot_idx],
ye=ye_trans_arr[:,ii,boot_idx],
me=me_trans_arr[:,ii,boot_idx],
x=x_trans_arr[:,boot_idx],
y=y_trans_arr[:,boot_idx],
m=m_trans_arr[:,boot_idx],
xe=xe_trans_arr[:,boot_idx],
ye=ye_trans_arr[:,boot_idx],
me=me_trans_arr[:,boot_idx],
t=np.tile(t_boot, (len(ref_table),1)))

# Now, do proper motion calculation, making sure to fix t0 to the
Expand All @@ -1278,7 +1299,8 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot

# Save proper motion fit results to output arrays
for col in motion_col_list:
motion_data[col][:,ii] = star_table[col]
motion_boot_sum[col] += star_table[col]
motion2_boot_sum[col] += star_table[col]**2

# Quick check to make sure bootstrap calc was valid: output t0 should be
# same as input t0_arr, since we used fixed_t0 option
Expand All @@ -1290,15 +1312,19 @@ def calc_bootstrap_errors(self, n_boot=100, boot_epochs_min=-1, calc_vel_in_boot
#print('=================================================')

# Calculate the bootstrap error values.
x_err_b = np.std(x_trans_arr, ddof=1, axis=1)
y_err_b = np.std(y_trans_arr, ddof=1, axis=1)
m_err_b = np.std(m_trans_arr, ddof=1, axis=1)
#pdb.set_trace()
x_boot_mean = x_boot_sum/n_boot
x_err_b = np.sqrt((x2_boot_sum - 2*x_boot_mean*x_boot_sum + n_boot*x_boot_mean**2)/n_boot)
y_boot_mean = y_boot_sum/n_boot
y_err_b = np.sqrt((y2_boot_sum - 2*y_boot_mean*y_boot_sum + n_boot*y_boot_mean**2)/n_boot)
m_boot_mean = m_boot_sum/n_boot
m_err_b = np.sqrt((m2_boot_sum - 2*m_boot_mean*m_boot_sum + n_boot*m_boot_mean**2)/n_boot)

motion_data_err = {}
if calc_vel_in_bootstrap:
for col in motion_col_list:
motion_data_err[col] = np.nanstd(motion_data[col], ddof=1,axis=1)
mot_boot_mean = motion_boot_sum[col]/n_boot
motion_data_err[col] = np.sqrt((motion2_boot_sum[col] -
2*mot_boot_mean*motion_boot_sum[col] + n_boot*mot_boot_mean**2)/n_boot)
else:
for col in motion_col_list:
motion_data_err[col] = np.nan
Expand Down
1 change: 1 addition & 0 deletions flystar/tests/test_align.py
Original file line number Diff line number Diff line change
Expand Up @@ -1141,6 +1141,7 @@ def test_bootstrap():
assert np.sum(np.isnan(match1.ref_table['ye_boot'])) == 0
assert np.sum(np.isnan(match1.ref_table['vx_err_boot'])) == 0
assert np.sum(np.isnan(match1.ref_table['vy_err_boot'])) == 0
#pdb.set_trace()

# Test 2: make sure boot_epochs_min is working
# Eliminate some rows to list2, so some stars are only in 1 epoch.
Expand Down