From a2822bf4c2388e0c0f605cbcc6c16752c45ccbd4 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:12:36 -0700 Subject: [PATCH 01/24] Refactored run_search to prepare for run_region_search --- analysis/analysis_utils.py | 260 +++++++++++++++++++++++++++++++++++-- analysis/run_search.py | 239 ++++++---------------------------- 2 files changed, 287 insertions(+), 212 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index edd77cbb1..ab47573ef 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -19,17 +19,17 @@ def kalman_filter(obs, var): xhatminus = np.zeros(len(obs)) Pminus = np.zeros(len(obs)) K = np.zeros(len(obs)) - + Q = 1. R = np.copy(var) - + xhat[0] = obs[0] P[0] = R[0] for k in range(1,len(obs)): xhatminus[k] = xhat[k-1] Pminus[k] = P[k-1] + Q - + K[k] = Pminus[k] / (Pminus[k] + R[k]) xhat[k] = xhatminus[k] + K[k]*(obs[k]-xhatminus[k]) P[k] = (1-K[k])*Pminus[k] @@ -47,14 +47,14 @@ def return_indices(psi_values, phi_values, val_on): inv_flux = np.array(phi_values[flux_idx]) inv_flux[inv_flux < -999.] = 9999999. f_var = (1./inv_flux) - + ## 1st pass #f_var = #var_curve[flux_idx]#np.var(fluxes)*np.ones(len(fluxes)) kalman_flux, kalman_error = kalman_filter(fluxes, f_var) if np.min(kalman_error) < 0.: return ([], [-1], []) deviations = np.abs(kalman_flux - fluxes) / kalman_error**.5 - + #print(deviations, fluxes) # keep_idx = np.where(deviations < 500.)[0] keep_idx = np.where(deviations < 5.)[0] @@ -91,9 +91,9 @@ def stamp_filter_parallel(stamps): if len(peak_2) > 1: peak_2 = np.max(np.abs(peak_2-10.)) - if ((mom_list[0] < 35.5) & (mom_list[1] < 35.5) & - (np.abs(mom_list[2]) < 1.) & - (np.abs(mom_list[3]) < .25) & (np.abs(mom_list[4]) < .25) & + if ((mom_list[0] < 35.5) & (mom_list[1] < 35.5) & + (np.abs(mom_list[2]) < 1.) & + (np.abs(mom_list[3]) < .25) & (np.abs(mom_list[4]) < .25) & (np.abs(peak_1 - 10.) < 2.) & (np.abs(peak_2 - 10.) < 2.)): keep_stamps = 1 else: @@ -108,14 +108,252 @@ def __init__(self): return - def cluster_results(self, results, x_size, y_size, + + def return_filename(self, visit_num): + + hits_file = 'v%i-fg.fits' % visit_num + + return hits_file + + def get_folder_visits(self, folder_visits): + + patch_visit_ids = np.array([int(visit_name[1:7]) for visit_name in folder_visits]) + + return patch_visit_ids + + def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): + """ + This function loads images and ingests them into a search object + + Input + --------- + + im_filepath : string + Image file path from which to load images + + time_file : string + File name containing image times + + Output + --------- + + search : kb.stack_search object + """ + + visit_nums, visit_times = np.genfromtxt(time_file, unpack=True) + image_time_dict = OrderedDict() + for visit_num, visit_time in zip(visit_nums, visit_times): + image_time_dict[str(int(visit_num))] = visit_time + + chunk_size = 100000 + + start = time.time() + + patch_visits = sorted(os.listdir(im_filepath)) + patch_visit_ids = self.get_folder_visits(patch_visits) + patch_visit_times = np.array([image_time_dict[str(visit_id)] for visit_id in patch_visit_ids]) + + if mjd_lims is None: + use_images = patch_visit_ids + else: + visit_only = np.where(((patch_visit_times > mjd_lims[0]) + & (patch_visit_times < mjd_lims[1])))[0] + print(visit_only) + use_images = patch_visit_ids[visit_only] + + image_mjd = np.array([image_time_dict[str(visit_id)] for visit_id in use_images]) + times = image_mjd - image_mjd[0] + + flags = ~0 # mask pixels with any flags + flag_exceptions = [32,39] # unless it has one of these special combinations of flags + master_flags = int('100111', 2) # mask any pixels which have any of + # these flags in more than two images + + hdulist = fits.open('%s/%s' % (im_filepath, self.return_filename(use_images[0]))) + w = WCS(hdulist[1].header) + ec_angle = self.calc_ecliptic_angle(w) + ec_angle = np.pi + 1.25 + del(hdulist) + + images = [kb.layered_image('%s/%s' % (im_filepath, self.return_filename(f))) for f in np.sort(use_images)] + + print('Images Loaded') + + p = kb.psf(psf_val) + stack = kb.image_stack(images) + + # Apply masks + stack.apply_mask_flags(flags, flag_exceptions) + stack.apply_master_mask(master_flags, 2) + + #stack.grow_mask() + #stack.grow_mask() + + #stack.apply_mask_threshold(120.) + + stack.set_times(times) + print("Times set") + + x_size = stack.get_width() + y_size = stack.get_width() + + search = kb.stack_search(stack, p) + + return(search) + + def process_results(search): + + keep = {'stamps': [], 'new_lh': [], 'results': [], 'times': [], + 'lc': [], 'final_results': []} + likelihood_limit = False + res_num = 0 + chunk_size = 500000 + print('---------------------------------------') + print("Processing Results") + print('---------------------------------------') + while likelihood_limit is False: + pool = mp.Pool(processes=16) + results = search.get_results(res_num,chunk_size) + chunk_headers = ("Chunk Start", "Chunk Size", "Chunk Max Likelihood", + "Chunk Min. Likelihood") + chunk_values = (res_num, len(keep['results']), results[0].lh, results[-1].lh) + for header, val, in zip(chunk_headers, chunk_values): + if type(val) == np.int: + print('%s = %i' % (header, val)) + else: + print('%s = %.2f' % (header, val)) + print('---------------------------------------') + psi_curves = [] + phi_curves = [] + print(search.get_results(0,10)) + for line in results: + psi_curve, phi_curve = search.lightcurve(line) + psi_curves.append(np.array(psi_curve).flatten()) + phi_curve = np.array(phi_curve).flatten() + phi_curve[phi_curve == 0.] = 99999999. + phi_curves.append(phi_curve) + if line.lh < likelihood_level: + likelihood_limit = True + break + keep_idx_results = pool.starmap_async(return_indices, + zip(psi_curves, phi_curves, + [j for j in range(len(psi_curves))])) + pool.close() + pool.join() + keep_idx_results = keep_idx_results.get() + + if len(keep_idx_results[0]) < 3: + keep_idx_results = [(0, [-1], 0.)] + + for result_on in range(len(psi_curves)): + + if keep_idx_results[result_on][1][0] == -1: + continue + elif len(keep_idx_results[result_on][1]) < 3: + continue + elif keep_idx_results[result_on][2] < likelihood_level: + continue + else: + keep_idx = keep_idx_results[result_on][1] + new_likelihood = keep_idx_results[result_on][2] + keep['results'].append(results[result_on]) + keep['new_lh'].append(new_likelihood) + stamps = search.sci_stamps(results[result_on], 10) + stamp_arr = np.array([np.array(stamps[s_idx]) for s_idx in keep_idx]) + keep['stamps'].append(np.sum(stamp_arr, axis=0)) + keep['lc'].append((psi_curves[result_on]/phi_curves[result_on])[keep_idx]) + keep['times'].append(image_mjd[keep_idx]) + print(len(keep['results'])) + + if len(keep['results']) > 500000: + with open('%s/memory_error_tr_%i_patch_%s.txt' % + (res_filepath, tract, patch), 'w') as f: + f.write('In %i total results, %i were kept. Needs manual look.' % + (res_num + chunk_size, len(keep['results']))) + memory_error = True + likelihood_limit = True + + if res_num+chunk_size >= 4000000: + likelihood_level = 20. + with open('%s/overload_error_tr_%i_patch_%s.txt' % + (res_filepath, tract, patch), 'w') as f: + f.write('In %i total results, %i were kept. Likelihood level down to %f.' % + (res_num + chunk_size, len(keep['results']), line.lh)) + + res_num += chunk_size + + return(keep) + + def cluster_results(self,keep): + lh_sorted_idx = np.argsort(np.array(keep['new_lh']))[::-1] + + print(len(lh_sorted_idx)) + if len(lh_sorted_idx) > 0: + print("Stamp filtering %i results" % len(lh_sorted_idx)) + pool = mp.Pool(processes=16) + stamp_filt_pool = pool.map_async(stamp_filter_parallel, + np.array(keep['stamps'])[lh_sorted_idx]) + pool.close() + pool.join() + stamp_filt_results = stamp_filt_pool.get() + stamp_filt_idx = lh_sorted_idx[np.where(np.array(stamp_filt_results) == 1)] + if len(stamp_filt_idx) > 0: + print("Clustering %i results" % len(stamp_filt_idx)) + cluster_idx = self.cluster_results(np.array(keep['results'])[stamp_filt_idx], + x_size, y_size, [vel_min, vel_max], + [ang_min, ang_max]) + keep['final_results'] = stamp_filt_idx[cluster_idx] + else: + cluster_idx = [] + keep['final_results'] = [] + del(cluster_idx) + del(stamp_filt_results) + del(stamp_filt_idx) + del(stamp_filt_pool) + else: + final_results = lh_sorted_idx + + print('Keeping %i results' % len(keep['final_results'])) + return(keep) + + def save_results(self, res_filepath, out_suffix, keep): + """ + Save results from a given search method + (either region search or grid search) + + Input + -------- + + res_filepath : string + + out_suffix : string + Suffix to append to the output file name + + keep : dictionary + Dictionary that contains the values to keep and print to filtering + """ + + np.savetxt('%s/results_%s.txt' % (res_filepath, out_suffix), + np.array(keep['results'])[keep['final_results']], fmt='%s') + with open('%s/lc_%s.txt' % (res_filepath, out_suffix), 'w') as f: + writer = csv.writer(f) + writer.writerows(np.array(keep['lc'])[keep['final_results']]) + with open('%s/times_%s.txt' % (res_filepath, out_suffix), 'w') as f: + writer = csv.writer(f) + writer.writerows(np.array(keep['times'])[keep['final_results']]) + np.savetxt('%s/filtered_likes_%s.txt' % (res_filepath, out_suffix), + np.array(keep['new_lh'])[keep['final_results']], fmt='%.4f') + np.savetxt('%s/ps_%s.txt' % (res_filepath, out_suffix), + np.array(keep['stamps']).reshape(len(keep['stamps']), 441)[keep['final_results']], fmt='%.4f') + + def cluster_results(self, results, x_size, y_size, v_lim, ang_lim, dbscan_args=None): default_dbscan_args = dict(eps=0.03, min_samples=-1, n_jobs=-1) if dbscan_args is not None: default_dbscan_args.update(dbscan_args) - dbscan_args = default_dbscan_args + dbscan_args = default_dbscan_args x_arr = [] y_arr = [] @@ -168,7 +406,7 @@ def calc_ecliptic_angle(self, test_wcs, angle_to_ecliptic=0.): vel_par_arr = [vel_par_arr] if type(vel_perp_arr) is not np.ndarray: vel_perp_arr = [vel_perp_arr] - for start_loc, vel_par, vel_perp in zip(pixel_start, + for start_loc, vel_par, vel_perp in zip(pixel_start, vel_par_arr, vel_perp_arr): start_coord = astroCoords.SkyCoord.from_pixel(start_loc[0], diff --git a/analysis/run_search.py b/analysis/run_search.py index b7b29f55c..f14ebba78 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -16,6 +16,30 @@ return_indices, stamp_filter_parallel from collections import OrderedDict +class run_region_search(analysis_utils): + + def __init(self,v_guess,radius): + + """ + Input + -------- + + v_guess : float + + Initial object velocity guess. Algorithm will search velocities + within 'radius' of 'v_guess' + + radius : float + + radius in velocity space to search, centered around 'v_guess' + """ + + return + + def run_search(self, im_filepath, res_filepath, out_suffix, time_file, + likelihood_level=10., mjd_lims=None): + # Load images to search + return class run_search(analysis_utils): @@ -39,88 +63,24 @@ def __init__(self, v_list, ang_list, num_obs): Number of images a trajectory must be unmasked. """ - + self.v_arr = np.array(v_list) self.ang_arr = np.array(ang_list) self.num_obs = num_obs return - def return_filename(self, visit_num): - - hits_file = 'v%i-fg.fits' % visit_num - - return hits_file - - def get_folder_visits(self, folder_visits): - - patch_visit_ids = np.array([int(visit_name[1:7]) for visit_name in folder_visits]) - - return patch_visit_ids - def run_search(self, im_filepath, res_filepath, out_suffix, time_file, likelihood_level=10., mjd_lims=None): - visit_nums, visit_times = np.genfromtxt(time_file, unpack=True) - image_time_dict = OrderedDict() - for visit_num, visit_time in zip(visit_nums, visit_times): - image_time_dict[str(int(visit_num))] = visit_time - - chunk_size = 100000 - - start = time.time() + # Initialize some values - patch_visits = sorted(os.listdir(im_filepath)) - patch_visit_ids = self.get_folder_visits(patch_visits) - patch_visit_times = np.array([image_time_dict[str(visit_id)] for visit_id in patch_visit_ids]) - - if mjd_lims is None: - use_images = patch_visit_ids - else: - visit_only = np.where(((patch_visit_times > mjd_lims[0]) - & (patch_visit_times < mjd_lims[1])))[0] - print(visit_only) - use_images = patch_visit_ids[visit_only] - - image_mjd = np.array([image_time_dict[str(visit_id)] for visit_id in use_images]) - times = image_mjd - image_mjd[0] - - flags = ~0 # mask pixels with any flags - flag_exceptions = [32,39] # unless it has one of these special combinations of flags - master_flags = int('100111', 2) # mask any pixels which have any of - # these flags in more than two images - - hdulist = fits.open('%s/%s' % (im_filepath, self.return_filename(use_images[0]))) - w = WCS(hdulist[1].header) - ec_angle = self.calc_ecliptic_angle(w) - ec_angle = np.pi + 1.25 - del(hdulist) - - images = [kb.layered_image('%s/%s' % (im_filepath, self.return_filename(f))) for f in np.sort(use_images)] -# for im in images: -# im.set_variance(kb.raw_image(np.median(im.variance())*np.ones(np.shape(im.variance())))) -# print(np.median(im.variance())) - print('Images Loaded') - - p = kb.psf(1.4) - stack = kb.image_stack(images) - del(images) - stack.apply_mask_flags(flags, flag_exceptions) - stack.apply_master_mask(master_flags, 2) - - #stack.grow_mask() - #stack.grow_mask() - - #stack.apply_mask_threshold(120.) + memory_error = False - stack.set_times(times) - print("Times set") + # Load images to search + search = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) - x_size = stack.get_width() - y_size = stack.get_width() - - search = kb.stack_search(stack, p) - del(stack) + # Run the grid search ang_min = ec_angle - self.ang_arr[0] ang_max = ec_angle + self.ang_arr[1] vel_min = self.v_arr[0] @@ -135,141 +95,18 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, search.gpu(int(self.ang_arr[2]),int(self.v_arr[2]),ang_min,ang_max, vel_min,vel_max,int(self.num_obs)) - keep_stamps = [] - keep_new_lh = [] - keep_results = [] - keep_times = [] - memory_error = False - keep_lc = [] - - likelihood_limit = False - res_num = 0 - chunk_size = 500000 - print('---------------------------------------') - print("Processing Results") - print('---------------------------------------') - while likelihood_limit is False: - pool = mp.Pool(processes=16) - results = search.get_results(res_num,chunk_size) - chunk_headers = ("Chunk Start", "Chunk Size", "Chunk Max Likelihood", - "Chunk Min. Likelihood") - chunk_values = (res_num, len(keep_results), results[0].lh, results[-1].lh) - for header, val, in zip(chunk_headers, chunk_values): - if type(val) == np.int: - print('%s = %i' % (header, val)) - else: - print('%s = %.2f' % (header, val)) - print('---------------------------------------') - psi_curves = [] - phi_curves = [] - print(search.get_results(0,10)) - for line in results: - psi_curve, phi_curve = search.lightcurve(line) - psi_curves.append(np.array(psi_curve).flatten()) - phi_curve = np.array(phi_curve).flatten() - phi_curve[phi_curve == 0.] = 99999999. - phi_curves.append(phi_curve) - if line.lh < likelihood_level: - likelihood_limit = True - break - keep_idx_results = pool.starmap_async(return_indices, - zip(psi_curves, phi_curves, - [j for j in range(len(psi_curves))])) - pool.close() - pool.join() - keep_idx_results = keep_idx_results.get() - - if len(keep_idx_results[0]) < 3: - keep_idx_results = [(0, [-1], 0.)] - - for result_on in range(len(psi_curves)): - - if keep_idx_results[result_on][1][0] == -1: - continue - elif len(keep_idx_results[result_on][1]) < 3: - continue - elif keep_idx_results[result_on][2] < likelihood_level: - continue - else: - keep_idx = keep_idx_results[result_on][1] - new_likelihood = keep_idx_results[result_on][2] - keep_results.append(results[result_on]) - keep_new_lh.append(new_likelihood) - stamps = search.sci_stamps(results[result_on], 10) - stamp_arr = np.array([np.array(stamps[s_idx]) for s_idx in keep_idx]) - keep_stamps.append(np.sum(stamp_arr, axis=0)) - keep_lc.append((psi_curves[result_on]/phi_curves[result_on])[keep_idx]) - keep_times.append(image_mjd[keep_idx]) - print(len(keep_results)) - - if len(keep_results) > 500000: - with open('%s/memory_error_tr_%i_patch_%s.txt' % - (res_filepath, tract, patch), 'w') as f: - f.write('In %i total results, %i were kept. Needs manual look.' % - (res_num + chunk_size, len(keep_results))) - memory_error = True - likelihood_limit = True - - if res_num+chunk_size >= 4000000: - likelihood_level = 20. - with open('%s/overload_error_tr_%i_patch_%s.txt' % - (res_filepath, tract, patch), 'w') as f: - f.write('In %i total results, %i were kept. Likelihood level down to %f.' % - (res_num + chunk_size, len(keep_results), line.lh)) - - res_num += chunk_size - + # Process the search results + keep = self.process_results(search) del(search) - lh_sorted_idx = np.argsort(np.array(keep_new_lh))[::-1] + # Cluster the results + keep = self.cluster_results(keep) - print(len(lh_sorted_idx)) - if len(lh_sorted_idx) > 0: - print("Stamp filtering %i results" % len(lh_sorted_idx)) - pool = mp.Pool(processes=16) - stamp_filt_pool = pool.map_async(stamp_filter_parallel, - np.array(keep_stamps)[lh_sorted_idx]) - pool.close() - pool.join() - stamp_filt_results = stamp_filt_pool.get() - stamp_filt_idx = lh_sorted_idx[np.where(np.array(stamp_filt_results) == 1)] - if len(stamp_filt_idx) > 0: - print("Clustering %i results" % len(stamp_filt_idx)) - cluster_idx = self.cluster_results(np.array(keep_results)[stamp_filt_idx], - x_size, y_size, [vel_min, vel_max], - [ang_min, ang_max]) - final_results = stamp_filt_idx[cluster_idx] - else: - cluster_idx = [] - final_results = [] - del(cluster_idx) - del(stamp_filt_results) - del(stamp_filt_idx) - del(stamp_filt_pool) - else: - final_results = lh_sorted_idx - - print('Keeping %i results' % len(final_results)) - - np.savetxt('%s/results_%s.txt' % (res_filepath, out_suffix), - np.array(keep_results)[final_results], fmt='%s') - with open('%s/lc_%s.txt' % (res_filepath, out_suffix), 'w') as f: - writer = csv.writer(f) - writer.writerows(np.array(keep_lc)[final_results]) - with open('%s/times_%s.txt' % (res_filepath, out_suffix), 'w') as f: - writer = csv.writer(f) - writer.writerows(np.array(keep_times)[final_results]) - np.savetxt('%s/filtered_likes_%s.txt' % (res_filepath, out_suffix), - np.array(keep_new_lh)[final_results], fmt='%.4f') - np.savetxt('%s/ps_%s.txt' % (res_filepath, out_suffix), - np.array(keep_stamps).reshape(len(keep_stamps), 441)[final_results], fmt='%.4f') + # Save the results + self.save_results(res_filepath, out_suffix, keep) end = time.time() - del(keep_stamps) - del(keep_times) - del(keep_results) - del(keep_new_lh) - del(keep_lc) - + del(keep) + print("Time taken for patch: ", end-start) From dc061a1fa3759aefee35b848d2c4b7aabba0a1a2 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:39:38 -0700 Subject: [PATCH 02/24] Added OrderedDict into analysis_utils.py --- analysis/analysis_utils.py | 2 ++ analysis/run_search.py | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index ab47573ef..dee3ab926 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -11,6 +11,8 @@ from astropy.wcs import WCS from sklearn.cluster import DBSCAN from skimage import measure +from collections import OrderedDict + def kalman_filter(obs, var): diff --git a/analysis/run_search.py b/analysis/run_search.py index f14ebba78..7bb406240 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -14,7 +14,6 @@ from skimage import measure from analysis_utils import analysis_utils, \ return_indices, stamp_filter_parallel -from collections import OrderedDict class run_region_search(analysis_utils): From 00d3f76252689a576c50381f6fb3745c841ece1b Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:43:58 -0700 Subject: [PATCH 03/24] Return ec_angle from load_images --- analysis/analysis_utils.py | 4 +++- analysis/run_search.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index dee3ab926..abe0aaae9 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -140,6 +140,8 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): --------- search : kb.stack_search object + + ec_angle : ecliptic angle """ visit_nums, visit_times = np.genfromtxt(time_file, unpack=True) @@ -201,7 +203,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): search = kb.stack_search(stack, p) - return(search) + return(search,ec_angle) def process_results(search): diff --git a/analysis/run_search.py b/analysis/run_search.py index 7bb406240..95bbf41a9 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -77,7 +77,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, memory_error = False # Load images to search - search = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) + search,ec_angle = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) # Run the grid search ang_min = ec_angle - self.ang_arr[0] From d1b025a1834c057145a66fc09c46965f67827f0c Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:46:15 -0700 Subject: [PATCH 04/24] Added self in process_results() --- analysis/analysis_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index abe0aaae9..b01378a35 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -205,7 +205,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): return(search,ec_angle) - def process_results(search): + def process_results(self,search): keep = {'stamps': [], 'new_lh': [], 'results': [], 'times': [], 'lc': [], 'final_results': []} From 93aed591c5cc36be527aecc634b3aade3246442e Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:49:54 -0700 Subject: [PATCH 05/24] Added likelihood_level to process_results() --- analysis/analysis_utils.py | 2 +- analysis/run_search.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index b01378a35..7ad3ef9a3 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -205,7 +205,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): return(search,ec_angle) - def process_results(self,search): + def process_results(self,search,likelihood_level): keep = {'stamps': [], 'new_lh': [], 'results': [], 'times': [], 'lc': [], 'final_results': []} diff --git a/analysis/run_search.py b/analysis/run_search.py index 95bbf41a9..1ed5bae51 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -95,7 +95,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, vel_min,vel_max,int(self.num_obs)) # Process the search results - keep = self.process_results(search) + keep = self.process_results(search,likelihood_level=likelihood_level) del(search) # Cluster the results From d198452db8bf84c11037e4031583b30b4eec12dd Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:50:58 -0700 Subject: [PATCH 06/24] Typo fix in run_search --- analysis/run_search.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/run_search.py b/analysis/run_search.py index 1ed5bae51..0a5bce226 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -95,7 +95,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, vel_min,vel_max,int(self.num_obs)) # Process the search results - keep = self.process_results(search,likelihood_level=likelihood_level) + keep = self.process_results(search,likelihood_level) del(search) # Cluster the results From a2604e0f6b14fc6621f7de362057f6f2479897d8 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:54:59 -0700 Subject: [PATCH 07/24] Renamed cluster_results to avoid duplicates --- analysis/analysis_utils.py | 2 +- analysis/run_search.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 7ad3ef9a3..d20dbb2c9 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -288,7 +288,7 @@ def process_results(self,search,likelihood_level): return(keep) - def cluster_results(self,keep): + def filter_results(self,keep): lh_sorted_idx = np.argsort(np.array(keep['new_lh']))[::-1] print(len(lh_sorted_idx)) diff --git a/analysis/run_search.py b/analysis/run_search.py index 0a5bce226..ccbc5bd6a 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -99,7 +99,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, del(search) # Cluster the results - keep = self.cluster_results(keep) + keep = self.filter_results(keep) # Save the results self.save_results(res_filepath, out_suffix, keep) From e5115a6c10955e4d4cdfa6a2226e4e7d37dea078 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 14:57:16 -0700 Subject: [PATCH 08/24] Import csv in analysis_utils.py --- analysis/analysis_utils.py | 1 + analysis/run_search.py | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index d20dbb2c9..e858cb393 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -4,6 +4,7 @@ import numpy as np import time import multiprocessing as mp +import csv import astropy.coordinates as astroCoords import astropy.units as u from kbmodpy import kbmod as kb diff --git a/analysis/run_search.py b/analysis/run_search.py index ccbc5bd6a..a50ff9da2 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -6,7 +6,6 @@ import multiprocessing as mp import astropy.coordinates as astroCoords import astropy.units as u -import csv from kbmodpy import kbmod as kb from astropy.io import fits from astropy.wcs import WCS From f94aa8de504aac8d630637649285af8c8c8e75b1 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 15:00:01 -0700 Subject: [PATCH 09/24] Moved start time back into run_search() --- analysis/analysis_utils.py | 4 +--- analysis/run_search.py | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index e858cb393..6f6bf1a04 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -151,9 +151,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): image_time_dict[str(int(visit_num))] = visit_time chunk_size = 100000 - - start = time.time() - + patch_visits = sorted(os.listdir(im_filepath)) patch_visit_ids = self.get_folder_visits(patch_visits) patch_visit_times = np.array([image_time_dict[str(visit_id)] for visit_id in patch_visit_ids]) diff --git a/analysis/run_search.py b/analysis/run_search.py index a50ff9da2..7bc3ec9b6 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -72,9 +72,9 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, likelihood_level=10., mjd_lims=None): # Initialize some values + start = time.time() memory_error = False - # Load images to search search,ec_angle = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) From 38aae4f5ef011b1e665fe9a1276931c1cdd2d101 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 15:03:10 -0700 Subject: [PATCH 10/24] moved final_results into dictionary keep --- analysis/analysis_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 6f6bf1a04..3de6f9362 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -151,7 +151,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): image_time_dict[str(int(visit_num))] = visit_time chunk_size = 100000 - + patch_visits = sorted(os.listdir(im_filepath)) patch_visit_ids = self.get_folder_visits(patch_visits) patch_visit_times = np.array([image_time_dict[str(visit_id)] for visit_id in patch_visit_ids]) @@ -314,7 +314,7 @@ def filter_results(self,keep): del(stamp_filt_idx) del(stamp_filt_pool) else: - final_results = lh_sorted_idx + keep['final_results'] = lh_sorted_idx print('Keeping %i results' % len(keep['final_results'])) return(keep) From b6b56793b5a253e47249ddb9e24f58e68d4c494d Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 15:16:57 -0700 Subject: [PATCH 11/24] Added image_mjd as an input to process_results() --- analysis/analysis_utils.py | 4 ++-- analysis/run_search.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 3de6f9362..233c1471f 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -202,9 +202,9 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): search = kb.stack_search(stack, p) - return(search,ec_angle) + return(search,ec_angle,image_mjd) - def process_results(self,search,likelihood_level): + def process_results(self,search,likelihood_level,image_mjd): keep = {'stamps': [], 'new_lh': [], 'results': [], 'times': [], 'lc': [], 'final_results': []} diff --git a/analysis/run_search.py b/analysis/run_search.py index 7bc3ec9b6..890cb90de 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -76,7 +76,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, memory_error = False # Load images to search - search,ec_angle = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) + search,ec_angle,image_mjd = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) # Run the grid search ang_min = ec_angle - self.ang_arr[0] @@ -94,7 +94,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, vel_min,vel_max,int(self.num_obs)) # Process the search results - keep = self.process_results(search,likelihood_level) + keep = self.process_results(search,likelihood_level,image_mjd) del(search) # Cluster the results From 497f0b4afa74d8a8775a43a7a66ba0af21ca4410 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:29:57 -0700 Subject: [PATCH 12/24] Created image_params dictionary --- analysis/analysis_utils.py | 31 ++++++++++++++++++------------- analysis/run_search.py | 12 +++++++++--- 2 files changed, 27 insertions(+), 16 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 233c1471f..afd742b71 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -142,9 +142,14 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): search : kb.stack_search object - ec_angle : ecliptic angle + image_params : dictionary + Contains image parameters such as ecliptic angle and mean Julian day """ + # Empty for now. Will contain x_size, y_size, ec_angle, and mjd before being returned. + image_params = {} + + visit_nums, visit_times = np.genfromtxt(time_file, unpack=True) image_time_dict = OrderedDict() for visit_num, visit_time in zip(visit_nums, visit_times): @@ -164,8 +169,8 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): print(visit_only) use_images = patch_visit_ids[visit_only] - image_mjd = np.array([image_time_dict[str(visit_id)] for visit_id in use_images]) - times = image_mjd - image_mjd[0] + image_params['mjd'] = np.array([image_time_dict[str(visit_id)] for visit_id in use_images]) + times = image_params['mjd'] - image_params['mjd'][0] flags = ~0 # mask pixels with any flags flag_exceptions = [32,39] # unless it has one of these special combinations of flags @@ -174,8 +179,8 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): hdulist = fits.open('%s/%s' % (im_filepath, self.return_filename(use_images[0]))) w = WCS(hdulist[1].header) - ec_angle = self.calc_ecliptic_angle(w) - ec_angle = np.pi + 1.25 + image_params['ec_angle'] = self.calc_ecliptic_angle(w) + image_params['ec_angle'] = np.pi + 1.25 del(hdulist) images = [kb.layered_image('%s/%s' % (im_filepath, self.return_filename(f))) for f in np.sort(use_images)] @@ -197,14 +202,14 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): stack.set_times(times) print("Times set") - x_size = stack.get_width() - y_size = stack.get_width() + image_params['x_size'] = stack.get_width() + image_params['y_size'] = stack.get_width() search = kb.stack_search(stack, p) - return(search,ec_angle,image_mjd) + return(search,image_params) - def process_results(self,search,likelihood_level,image_mjd): + def process_results(self,search,image_params,likelihood_level): keep = {'stamps': [], 'new_lh': [], 'results': [], 'times': [], 'lc': [], 'final_results': []} @@ -265,7 +270,7 @@ def process_results(self,search,likelihood_level,image_mjd): stamp_arr = np.array([np.array(stamps[s_idx]) for s_idx in keep_idx]) keep['stamps'].append(np.sum(stamp_arr, axis=0)) keep['lc'].append((psi_curves[result_on]/phi_curves[result_on])[keep_idx]) - keep['times'].append(image_mjd[keep_idx]) + keep['times'].append(image_params['mjd'][keep_idx]) print(len(keep['results'])) if len(keep['results']) > 500000: @@ -287,7 +292,7 @@ def process_results(self,search,likelihood_level,image_mjd): return(keep) - def filter_results(self,keep): + def filter_results(self,keep,image_params): lh_sorted_idx = np.argsort(np.array(keep['new_lh']))[::-1] print(len(lh_sorted_idx)) @@ -303,8 +308,8 @@ def filter_results(self,keep): if len(stamp_filt_idx) > 0: print("Clustering %i results" % len(stamp_filt_idx)) cluster_idx = self.cluster_results(np.array(keep['results'])[stamp_filt_idx], - x_size, y_size, [vel_min, vel_max], - [ang_min, ang_max]) + image_params['x_size'], image_params['y_size'], + image_params['v_lim'], image_params['ang_lim']) keep['final_results'] = stamp_filt_idx[cluster_idx] else: cluster_idx = [] diff --git a/analysis/run_search.py b/analysis/run_search.py index 890cb90de..a928770f4 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -76,13 +76,19 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, memory_error = False # Load images to search - search,ec_angle,image_mjd = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) + search,ec_angle,image_params = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) # Run the grid search + + # Set min and max values for angle and velocity ang_min = ec_angle - self.ang_arr[0] ang_max = ec_angle + self.ang_arr[1] vel_min = self.v_arr[0] vel_max = self.v_arr[1] + # Save values in image_params for use in filter_results + image_params['ang_lims'] = [ang_min, ang_max] + image_params['vel_lims'] = [vel_min, vel_max] + print("Starting Search") print('---------------------------------------') param_headers = ("Ecliptic Angle", "Min. Search Angle", "Max Search Angle", @@ -94,11 +100,11 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, vel_min,vel_max,int(self.num_obs)) # Process the search results - keep = self.process_results(search,likelihood_level,image_mjd) + keep = self.process_results(search,image_params,likelihood_level) del(search) # Cluster the results - keep = self.filter_results(keep) + keep = self.filter_results(keep,image_params) # Save the results self.save_results(res_filepath, out_suffix, keep) From edd48db8a825886a2a0d9c106fb14b008cdb36e2 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:37:54 -0700 Subject: [PATCH 13/24] Typo fix in call to load_images --- analysis/run_search.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/run_search.py b/analysis/run_search.py index a928770f4..5d3ed45ca 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -76,7 +76,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, memory_error = False # Load images to search - search,ec_angle,image_params = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) + search,image_params = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) # Run the grid search @@ -88,7 +88,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, # Save values in image_params for use in filter_results image_params['ang_lims'] = [ang_min, ang_max] image_params['vel_lims'] = [vel_min, vel_max] - + print("Starting Search") print('---------------------------------------') param_headers = ("Ecliptic Angle", "Min. Search Angle", "Max Search Angle", From 6eb27dc28d55d8db813252331bc7a1010525f721 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:41:31 -0700 Subject: [PATCH 14/24] Fixed call to ec_angle (now image_params['ec_angle'] --- analysis/run_search.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/run_search.py b/analysis/run_search.py index 5d3ed45ca..82d9ad507 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -81,8 +81,8 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, # Run the grid search # Set min and max values for angle and velocity - ang_min = ec_angle - self.ang_arr[0] - ang_max = ec_angle + self.ang_arr[1] + ang_min = image_params['ec_angle'] - self.ang_arr[0] + ang_max = image_params['ec_angle'] + self.ang_arr[1] vel_min = self.v_arr[0] vel_max = self.v_arr[1] # Save values in image_params for use in filter_results From 26f09ec46fec3b4ff439754c09406065770b1210 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:46:29 -0700 Subject: [PATCH 15/24] Incorporate image_params as function inputs --- analysis/run_search.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/analysis/run_search.py b/analysis/run_search.py index 82d9ad507..b94b26b3f 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -93,11 +93,11 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, print('---------------------------------------') param_headers = ("Ecliptic Angle", "Min. Search Angle", "Max Search Angle", "Min Velocity", "Max Velocity") - param_values = (ec_angle, ang_min, ang_max, vel_min, vel_max) + param_values = (image_params['ec_angle'], *image_params['ang_lims'], *image_params['vel_lims']) for header, val in zip(param_headers, param_values): print('%s = %.4f' % (header, val)) - search.gpu(int(self.ang_arr[2]),int(self.v_arr[2]),ang_min,ang_max, - vel_min,vel_max,int(self.num_obs)) + search.gpu(int(self.ang_arr[2]),int(self.v_arr[2]),*image_params['ang_lim'], + *image_params['vel_lims'],int(self.num_obs)) # Process the search results keep = self.process_results(search,image_params,likelihood_level) From 01448f3f7e558836693bd9146517cb61a7c2c1ef Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:48:12 -0700 Subject: [PATCH 16/24] Typo: ang_lim -> ang_lims --- analysis/run_search.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/run_search.py b/analysis/run_search.py index b94b26b3f..7eb8ce3ab 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -96,7 +96,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, param_values = (image_params['ec_angle'], *image_params['ang_lims'], *image_params['vel_lims']) for header, val in zip(param_headers, param_values): print('%s = %.4f' % (header, val)) - search.gpu(int(self.ang_arr[2]),int(self.v_arr[2]),*image_params['ang_lim'], + search.gpu(int(self.ang_arr[2]),int(self.v_arr[2]),*image_params['ang_lims'], *image_params['vel_lims'],int(self.num_obs)) # Process the search results From 45fe62a5375521b39177b08fc977473355900732 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:51:38 -0700 Subject: [PATCH 17/24] Typo: vel_lim -> vel_lims --- analysis/analysis_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index afd742b71..f466a5e2e 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -309,7 +309,7 @@ def filter_results(self,keep,image_params): print("Clustering %i results" % len(stamp_filt_idx)) cluster_idx = self.cluster_results(np.array(keep['results'])[stamp_filt_idx], image_params['x_size'], image_params['y_size'], - image_params['v_lim'], image_params['ang_lim']) + image_params['v_lims'], image_params['ang_lims']) keep['final_results'] = stamp_filt_idx[cluster_idx] else: cluster_idx = [] From 0e8cab74e3908d1d2e4bffa89d24c392cd04d897 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 16:55:22 -0700 Subject: [PATCH 18/24] Typo: v_lims->vel_lims --- analysis/analysis_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index f466a5e2e..8ada4aafb 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -309,7 +309,7 @@ def filter_results(self,keep,image_params): print("Clustering %i results" % len(stamp_filt_idx)) cluster_idx = self.cluster_results(np.array(keep['results'])[stamp_filt_idx], image_params['x_size'], image_params['y_size'], - image_params['v_lims'], image_params['ang_lims']) + image_params['vel_lims'], image_params['ang_lims']) keep['final_results'] = stamp_filt_idx[cluster_idx] else: cluster_idx = [] From 5f0c3cc67eede96b6be2a5e15036d615ff145a81 Mon Sep 17 00:00:00 2001 From: smotherh Date: Fri, 11 May 2018 17:11:15 -0700 Subject: [PATCH 19/24] Removed hard-coded ec_angle value --- analysis/analysis_utils.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 8ada4aafb..76a3495d7 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -180,7 +180,6 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): hdulist = fits.open('%s/%s' % (im_filepath, self.return_filename(use_images[0]))) w = WCS(hdulist[1].header) image_params['ec_angle'] = self.calc_ecliptic_angle(w) - image_params['ec_angle'] = np.pi + 1.25 del(hdulist) images = [kb.layered_image('%s/%s' % (im_filepath, self.return_filename(f))) for f in np.sort(use_images)] @@ -210,6 +209,9 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): return(search,image_params) def process_results(self,search,image_params,likelihood_level): + """ + Processes results that are output by the gpu search. + """ keep = {'stamps': [], 'new_lh': [], 'results': [], 'times': [], 'lc': [], 'final_results': []} @@ -293,6 +295,21 @@ def process_results(self,search,image_params,likelihood_level): return(keep) def filter_results(self,keep,image_params): + """ + Filters results from a given search and clusters duplicate detections + + Input + --------- + + keep : Dictionary + Contains the values of which results were kept from the search algorithm + + image_params : dictionary + Contains values concerning the image and search initial settings + + Output + --------- + """ lh_sorted_idx = np.argsort(np.array(keep['new_lh']))[::-1] print(len(lh_sorted_idx)) From 3e62d3567fb955c699e6254c5648e2ddb29bc69d Mon Sep 17 00:00:00 2001 From: smotherh Date: Mon, 14 May 2018 12:26:49 -0700 Subject: [PATCH 20/24] Uncommented out mask commands --- analysis/analysis_utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 76a3495d7..feba785c1 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -193,10 +193,10 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): stack.apply_mask_flags(flags, flag_exceptions) stack.apply_master_mask(master_flags, 2) - #stack.grow_mask() - #stack.grow_mask() + stack.grow_mask() + stack.grow_mask() - #stack.apply_mask_threshold(120.) + stack.apply_mask_threshold(120.) stack.set_times(times) print("Times set") From e4160b89521865e1475b82c7ba5f743d025030da Mon Sep 17 00:00:00 2001 From: smotherh Date: Mon, 14 May 2018 12:34:54 -0700 Subject: [PATCH 21/24] Output number of images loaded to console --- analysis/analysis_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index feba785c1..398802a21 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -184,7 +184,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): images = [kb.layered_image('%s/%s' % (im_filepath, self.return_filename(f))) for f in np.sort(use_images)] - print('Images Loaded') + print('Loaded %i images' %(len(images))) p = kb.psf(psf_val) stack = kb.image_stack(images) From de9e4878cfc8fab3cec94ab7d53a5fb8e7e92161 Mon Sep 17 00:00:00 2001 From: smotherh Date: Wed, 16 May 2018 11:09:29 -0700 Subject: [PATCH 22/24] Fixed indentation error and stared region_search class --- analysis/analysis_utils.py | 5 +++- analysis/run_search.py | 48 ++++++++++++++++++++++++++++++++++---- 2 files changed, 47 insertions(+), 6 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 398802a21..97ae988fd 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -292,7 +292,7 @@ def process_results(self,search,image_params,likelihood_level): res_num += chunk_size - return(keep) + return(keep) def filter_results(self,keep,image_params): """ @@ -309,6 +309,9 @@ def filter_results(self,keep,image_params): Output --------- + + keep : Dictionary + Contains the values of which results were kept from the search algorithm """ lh_sorted_idx = np.argsort(np.array(keep['new_lh']))[::-1] diff --git a/analysis/run_search.py b/analysis/run_search.py index 7eb8ce3ab..60f7509f5 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -14,29 +14,67 @@ from analysis_utils import analysis_utils, \ return_indices, stamp_filter_parallel -class run_region_search(analysis_utils): +class region_search(analysis_utils): - def __init(self,v_guess,radius): + def __init__(self,v_guess,radius,num_obs): """ Input -------- - v_guess : float + v_guess : float array - Initial object velocity guess. Algorithm will search velocities - within 'radius' of 'v_guess' + Initial object velocity guess. Given as an array or tuple. + Algorithm will search velocities within 'radius' of 'v_guess' radius : float radius in velocity space to search, centered around 'v_guess' + + num_obs : int + + The minimum number of observations required to keep the object """ + self.v_guess = v_guess + self.radius = radius + self.num_obs = num_obs + return def run_search(self, im_filepath, res_filepath, out_suffix, time_file, likelihood_level=10., mjd_lims=None): + # Initialize some values + start = time.time() + + memory_error = False # Load images to search + search,image_params = self.load_images(im_filepath,time_file,mjd_lims=mjd_lims) + + # Run the region search + # Save values in image_params for use in filter_results + + print("Starting Search") + print('---------------------------------------') + param_headers = ("Central Velocity Guess","Radius in velocity space") + param_values = (*self.v_guess,self.radius) + for header, val in zip(param_headers, param_values): + print('%s = %.4f' % (header, val)) + search.region_search(*self.v_guess, self.radius, likelihood_level, int(self.num_obs)) + + # Process the search results + keep = self.process_results(search,image_params,likelihood_level) + del(search) + + # Cluster the results + #keep = self.filter_results(keep,image_params) + + # Save the results + self.save_results(res_filepath, out_suffix, keep) + + end = time.time() + + del(keep) return class run_search(analysis_utils): From 200dcc9eba48125a6ca6da713809ed42b954dfe2 Mon Sep 17 00:00:00 2001 From: smotherh Date: Wed, 16 May 2018 11:59:59 -0700 Subject: [PATCH 23/24] Removed uninitialized variables for file names --- analysis/analysis_utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 97ae988fd..6164f3852 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -208,7 +208,7 @@ def load_images(self, im_filepath, time_file,psf_val=1.4, mjd_lims=None): return(search,image_params) - def process_results(self,search,image_params,likelihood_level): + def process_results(self,search,image_params,res_filepath,likelihood_level): """ Processes results that are output by the gpu search. """ @@ -276,8 +276,8 @@ def process_results(self,search,image_params,likelihood_level): print(len(keep['results'])) if len(keep['results']) > 500000: - with open('%s/memory_error_tr_%i_patch_%s.txt' % - (res_filepath, tract, patch), 'w') as f: + with open('%s/memory_error.txt' % + (res_filepath), 'w') as f: f.write('In %i total results, %i were kept. Needs manual look.' % (res_num + chunk_size, len(keep['results']))) memory_error = True @@ -285,8 +285,8 @@ def process_results(self,search,image_params,likelihood_level): if res_num+chunk_size >= 4000000: likelihood_level = 20. - with open('%s/overload_error_tr_%i_patch_%s.txt' % - (res_filepath, tract, patch), 'w') as f: + with open('%s/overload_error.txt' % + (res_filepath), 'w') as f: f.write('In %i total results, %i were kept. Likelihood level down to %f.' % (res_num + chunk_size, len(keep['results']), line.lh)) From 33cefeaa46b469446920551abfa175bf325b413f Mon Sep 17 00:00:00 2001 From: smotherh Date: Wed, 16 May 2018 12:01:27 -0700 Subject: [PATCH 24/24] Added res_path to call to process_results() --- analysis/run_search.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/run_search.py b/analysis/run_search.py index 60f7509f5..17e1b0f2e 100644 --- a/analysis/run_search.py +++ b/analysis/run_search.py @@ -138,7 +138,7 @@ def run_search(self, im_filepath, res_filepath, out_suffix, time_file, *image_params['vel_lims'],int(self.num_obs)) # Process the search results - keep = self.process_results(search,image_params,likelihood_level) + keep = self.process_results(search,image_params,res_filepath,likelihood_level) del(search) # Cluster the results