diff --git a/popsycle/lightcurves.py b/popsycle/lightcurves.py index f36ddd6..17618af 100644 --- a/popsycle/lightcurves.py +++ b/popsycle/lightcurves.py @@ -238,7 +238,7 @@ def coords_and_prop_motion(event): return raL, decL, muL, muS -def get_pspl_lightcurve_parameters(events, photometric_system, filter_name, event_id = None): +def get_pspl_lightcurve_parameters(events, filter_dict, event_id = None): """ Find the parameters for PSPL_PhotAstrom_Par_Param1 from event_table. @@ -296,8 +296,12 @@ def get_pspl_lightcurve_parameters(events, photometric_system, filter_name, even beta = event['u0'] * event['theta_E'] # mas dL = event['rad_L'] * 10 ** 3 # Distance to lens dS = event['rad_S'] * 10 ** 3 # Distance to source - mag_src = [event['%s_%s_app_S' % (photometric_system, filter_name)]] - b_sff = [event['f_blend_%s' % filter_name]] + mag_src = [] + b_sff = [] + for photometric_system in filter_dict: + for filter_name in filter_dict[photometric_system]: + mag_src += [event['%s_%s_app_S' % (photometric_system, filter_name)]] + b_sff += [event['f_blend_%s' % filter_name]] parameter_dict = {'raL': raL, 'decL': decL, 'mL': mL, 't0': t0, 'beta': beta, 'dL': dL, 'dL_dS': dL / dS, @@ -307,7 +311,7 @@ def get_pspl_lightcurve_parameters(events, photometric_system, filter_name, even return parameter_dict, obj_id_L, obj_id_S, model_name -def get_psbl_lightcurve_parameters(events, companions, comp_idx_L, photometric_system, filter_name, event_id = None): +def get_psbl_lightcurve_parameters(events, companions, comp_idx_L, filter_dict, event_id = None): """ Find the parameters for PSBL_PhotAstrom_Par_EllOrbs_Param7 from event_table and comp_table. @@ -379,16 +383,27 @@ def get_psbl_lightcurve_parameters(events, companions, comp_idx_L, photometric_s beta_p = event['u0'] * event['theta_E'] # 5.0 dL = event['rad_L'] * 10 ** 3 # Distance to lens dS = event['rad_S'] * 10 ** 3 # Distance to source - mag_src = [event['%s_%s_app_S' % (photometric_system, filter_name)]] - b_sff = [event['f_blend_%s' % filter_name]] # ASSUMES ALL BINARY LENSES ARE BLENDED +# mag_src = [event['%s_%s_app_S' % (photometric_system, filter_name)]] +# b_sff = [event['f_blend_%s' % filter_name]] # ASSUMES ALL BINARY LENSES ARE BLENDED omega = companions['omega'][comp_idx_L] big_omega = companions['Omega'][comp_idx_L] i = companions['i'][comp_idx_L] e = companions['e'][comp_idx_L] tp = companions['tp'][comp_idx_L] a = 10**(companions['log_a'][comp_idx_L]) - dmag_Lp_Ls = [event['%s_%s_L' % (photometric_system, filter_name)] - companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_L]] +# dmag_Lp_Ls = [event['%s_%s_L' % (photometric_system, filter_name)] - companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_L]] + #FIX HERE + mag_src = [] + b_sff = [] + dmag_Lp_Ls = [] + for photometric_system in filter_dict: + for filter_name in filter_dict[photometric_system]: + filt = photometric_system+'_'+filter_name + mag_src += [event['%s_%s_app_S' % (photometric_system, filter_name)]] + b_sff += [event['f_blend_%s' % filter_name]] + dmag_Lp_Ls += [event['%s_%s_L' % (photometric_system, filter_name)] - companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_L]] + parameter_dict = {'raL': raL, 'decL': decL, 'mLp': mLp, 'mLs': mLs, 't0_p': t0_p, @@ -401,7 +416,7 @@ def get_psbl_lightcurve_parameters(events, companions, comp_idx_L, photometric_s return parameter_dict, obj_id_L, obj_id_S, model_name -def get_bspl_lightcurve_parameters(events, companions, comp_idx_S, photometric_system, filter_name, red_law, event_id = None): +def get_bspl_lightcurve_parameters(events, companions, comp_idx_S, filter_dict, red_law, event_id = None): """ Find the parameters for BSPL_PhotAstrom_Par_EllOrbs_Param4 from event_table and comp_table. @@ -467,9 +482,9 @@ def get_bspl_lightcurve_parameters(events, companions, comp_idx_S, photometric_s model_name = 'BSPL_PhotAstrom_Par_EllOrbs_Param4' - filt_dict = phot_utils.make_filt_dict() - f_i = filt_dict[photometric_system + '_' + filter_name][red_law] - abs_mag_sec = companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_S] +# filt_dict = phot_utils.make_filt_dict() +# f_i = filt_dict[photometric_system + '_' + filter_name][red_law] +# abs_mag_sec = companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_S] mL = event['mass_L'] # msun (Lens current mass) t0_p = event['t0'] # mjd @@ -477,10 +492,10 @@ def get_bspl_lightcurve_parameters(events, companions, comp_idx_S, photometric_s dL = event['rad_L'] * 10 ** 3 # Distance to lens dL_dS = dL / (event['rad_S'] * 10 ** 3) # Distance to lens/Distance to source xS0 = np.array([0, 0]) # arbitrary offset (arcsec) - mag_src_sec = synthetic.calc_app_mag(event['rad_S'], abs_mag_sec, event['exbv_S'], f_i) - mag_src_pri = binary_utils.subtract_magnitudes( - event['%s_%s_app_S' % (photometric_system, filter_name)], mag_src_sec) - b_sff = event['f_blend_%s' % filter_name] # ASSUMES THAT SOURCE BINARIES ARE BLENDED +# mag_src_sec = synthetic.calc_app_mag(event['rad_S'], abs_mag_sec, event['exbv_S'], f_i) +# mag_src_pri = binary_utils.subtract_magnitudes( +# event['%s_%s_app_S' % (photometric_system, filter_name)], mag_src_sec) +# b_sff = event['f_blend_%s' % filter_name] # ASSUMES THAT SOURCE BINARIES ARE BLENDED omega = companions['omega'][comp_idx_S] big_omega = companions['Omega'][comp_idx_S] i = companions['i'][comp_idx_S] @@ -490,16 +505,33 @@ def get_bspl_lightcurve_parameters(events, companions, comp_idx_S, photometric_s mass_source_p = event['mass_S'] mass_source_s = companions['mass'][comp_idx_S] + #FIX HERE + b_sff = [] + mag_src_sec = [] + mag_src_pri = [] + filt_dict = phot_utils.make_filt_dict() + for photometric_system in filter_dict: + for filter_name in filter_dict[photometric_system]: + filt = photometric_system+'_'+filter_name + f_i = filt_dict[photometric_system + '_' + filter_name][red_law] + abs_mag_sec = companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_S] + mag_src_sec_tmp = synthetic.calc_app_mag(event['rad_S'], abs_mag_sec, event['exbv_S'], f_i) + + mag_src_sec += [mag_src_sec_tmp] + mag_src_pri += [binary_utils.subtract_magnitudes( + event['%s_%s_app_S' % (photometric_system, filter_name)], mag_src_sec_tmp)] + b_sff += [event['f_blend_%s' % filter_name]] + parameter_dict = {'raL': raL, 'decL': decL, 'mL': mL, 't0': t0_p, 'beta': beta_p, 'dL': dL, 'dL_dS': dL_dS, 'xS0_E': xS0[0], 'xS0_N': xS0[1], 'muL_E': muL[0], 'muL_N': muL[1], 'muS_E': muS[0], 'muS_N': muS[1], - 'mag_src_pri': [mag_src_pri], 'mag_src_sec': [mag_src_sec], 'b_sff': [b_sff], 'omega_pri': omega, 'big_omega_sec': big_omega, + 'mag_src_pri': mag_src_pri, 'mag_src_sec': mag_src_sec, 'b_sff': b_sff, 'omega_pri': omega, 'big_omega_sec': big_omega, 'i': i, 'e': e, 'log_a': log_a, 'mass_source_p': mass_source_p, 'mass_source_s': mass_source_s, 'tp': tp} return parameter_dict, obj_id_L, obj_id_S, model_name -def get_bsbl_lightcurve_parameters(events, companions, comp_idx_L, comp_idx_S, photometric_system, filter_name, red_law, event_id = None): +def get_bsbl_lightcurve_parameters(events, companions, comp_idx_L, comp_idx_S, filter_dict, red_law, event_id = None): """ Find the parameters for BSBL_PhotAstrom_Par_EllOrbs_Param3 from event_table and comp_table. @@ -569,8 +601,8 @@ def get_bsbl_lightcurve_parameters(events, companions, comp_idx_L, comp_idx_S, p model_name = 'BSBL_PhotAstrom_Par_EllOrbs_Param3' - f_i = synthetic.filt_dict[photometric_system + '_' + filter_name][red_law] - abs_mag_sec = companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_S] +# f_i = synthetic.filt_dict[photometric_system + '_' + filter_name][red_law] +# abs_mag_sec = companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_S] mLp = event['mass_L'] # msun (Lens current mass) mLs = companions['mass'][comp_idx_L] # msun (Companion lens current mass) @@ -580,10 +612,10 @@ def get_bsbl_lightcurve_parameters(events, companions, comp_idx_L, comp_idx_S, p dS = event['rad_S'] * 10 ** 3 # Distance to source xS0_E = 0.0 # arbitrary offset (arcsec) xS0_N = 0.0 # arbitrary offset (arcsec) - mag_src_sec = synthetic.calc_app_mag(event['rad_S'], abs_mag_sec, event['exbv_S'], f_i) - mag_src_pri = binary_utils.subtract_magnitudes( - event['%s_%s_app_S' % (photometric_system, filter_name)], mag_src_sec) - b_sff = event['f_blend_%s' % filter_name] # ASSUMES THAT SOURCE BINARIES ARE BLENDED +# mag_src_sec = synthetic.calc_app_mag(event['rad_S'], abs_mag_sec, event['exbv_S'], f_i) +# mag_src_pri = binary_utils.subtract_magnitudes( +# event['%s_%s_app_S' % (photometric_system, filter_name)], mag_src_sec) +# b_sff = event['f_blend_%s' % filter_name] # ASSUMES THAT SOURCE BINARIES ARE BLENDED omegaL = companions['omega'][comp_idx_L] big_omegaL = companions['Omega'][comp_idx_L] iL = companions['i'][comp_idx_L] @@ -596,15 +628,35 @@ def get_bsbl_lightcurve_parameters(events, companions, comp_idx_L, comp_idx_S, p eS = companions['e'][comp_idx_S] tpS = companions['tp'][comp_idx_S] aS = 10**(companions['log_a'][comp_idx_S]) - dmag_Lp_Ls = [event['%s_%s_L' % (photometric_system, filter_name)] - companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_L]] +# dmag_Lp_Ls = [event['%s_%s_L' % (photometric_system, filter_name)] - companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_L]] mass_source_p = event['mass_S'] mass_source_s = companions['mass'][comp_idx_S] + + #FIX HERE + b_sff = [] + mag_src_sec = [] + mag_src_pri = [] + dmag_Lp_Ls = [] + filt_dict = phot_utils.make_filt_dict() + for photometric_system in filter_dict: + for filter_name in filter_dict[photometric_system]: + filt = photometric_system+'_'+filter_name + f_i = filt_dict[photometric_system + '_' + filter_name][red_law] + abs_mag_sec = companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_S] + + mag_src_sec_tmp = synthetic.calc_app_mag(event['rad_S'], abs_mag_sec, event['exbv_S'], f_i) + mag_src_sec += [mag_src_sec_tmp] + mag_src_pri += [binary_utils.subtract_magnitudes( + event['%s_%s_app_S' % (photometric_system, filter_name)], mag_src_sec_tmp)] + b_sff += [event['f_blend_%s' % filter_name]] + dmag_Lp_Ls += [event['%s_%s_L' % (photometric_system, filter_name)] - companions['m_%s_%s' % (photometric_system, filter_name)][comp_idx_L]] + parameter_dict = {'raL': raL, 'decL': decL, 'mLp': mLp, 'mLs': mLs, 't0_p': t0_p, 'xS0_E': xS0_E, 'xS0_N': xS0_N, 'beta_p': beta_p, 'muL_E': muL[0], 'muL_N': muL[1], 'muS_E': muS[0], 'muS_N': muS[1], 'dL': dL, 'dS': dS, - 'mag_src_pri': [mag_src_pri], 'mag_src_sec': [mag_src_sec], 'b_sff': [b_sff], + 'mag_src_pri': mag_src_pri, 'mag_src_sec': mag_src_sec, 'b_sff': b_sff, 'omegaL_pri': omegaL, 'big_omegaL_sec': big_omegaL, 'iL': iL, 'eL': eL, 'tpL': tpL, 'aL': aL, 'omegaS_pri': omegaS, 'big_omegaS_sec': big_omegaS, 'iS': iS, 'eS': eS, 'tpS': tpS, 'aS': aS, 'dmag_Lp_Ls': dmag_Lp_Ls, 'mass_source_p': mass_source_p, 'mass_source_s': mass_source_s} diff --git a/popsycle/run.py b/popsycle/run.py index 0ce5816..d24ec67 100755 --- a/popsycle/run.py +++ b/popsycle/run.py @@ -745,8 +745,7 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, filter_name = popsycle_config['filter_dict'][photometric_system][0] _check_refine_binary_events(events=filename_dict['refined_events_filename'], companions=refined_events_comp_filename, - filter_name=filter_name, - photometric_system=photometric_system, + filter_dict=popsycle_config['filter_dict'], n_proc=n_cores_refine_binary_events, overwrite=overwrite, output_file='default', save_phot=True, @@ -1131,8 +1130,7 @@ def run(output_root='root0', filter_name = popsycle_config['filter_dict'][photometric_system][0] _check_refine_binary_events(events=filename_dict['refined_events_filename'], companions=refined_events_comp_filename, - filter_name=filter_name, - photometric_system=photometric_system, + filter_dict=popsycle_config['filter_dict'], n_proc=n_cores_refine_binary_events, overwrite=overwrite, output_file='default', save_phot=True, @@ -1270,8 +1268,7 @@ def run(output_root='root0', phot_dir = '%s_bin_phot' % output_root synthetic.refine_binary_events(events=filename_dict['refined_events_filename'], companions=refined_events_comp_filename, - filter_name=filter_name, - photometric_system=photometric_system, + filter_dict=popsycle_config['filter_dict'], n_proc=n_cores_refine_binary_events, overwrite=overwrite, output_file='default', save_phot=True, @@ -1328,7 +1325,7 @@ def main(): 'Default is --n-cores=1 or serial processing.', default=1) optional.add_argument('--multi-proc-refine-binary-events', type=bool, - help='Controls multi processing for refine bianry events ' + help='Controls multi processing for refine bianary events ' 'even if n-cores=1', default=True) optional.add_argument('--seed', type=int, diff --git a/popsycle/synthetic.py b/popsycle/synthetic.py index 2a6dbbf..ddcc359 100755 --- a/popsycle/synthetic.py +++ b/popsycle/synthetic.py @@ -5226,7 +5226,7 @@ def _add_multiples_parameters(companion_table, event_table): def _check_refine_binary_events(events, companions, - photometric_system, filter_name, + filter_dict, overwrite, output_file, save_phot, phot_dir, n_proc, multi_proc): """ @@ -5280,11 +5280,16 @@ def _check_refine_binary_events(events, companions, if not isinstance(companions, str): raise Exception('companions (%s) must be a string.' % str(companions)) - if not isinstance(filter_name, str): - raise Exception('filter_name (%s) must be a string.' % str(filter_name)) + if filter_dict is not None: + if not isinstance(filter_dict, dict): + raise Exception('filter_dict (%s) must be a dictionary.' % str(filter_dict)) + if not all(isinstance(key,str) for key in filter_dict): + raise Exception('All filter_dict keys must be strings.') + if not all(isinstance(filt,str) for key,val in filter_dict.items() for filt in val): + raise Exception('All filter_dict vaues must be lists of strings.') - if not isinstance(photometric_system, str): - raise Exception('photometric_system (%s) must be a string.' % str(photometric_system)) +# if not isinstance(photometric_system, str): +# raise Exception('photometric_system (%s) must be a string.' % str(photometric_system)) if not isinstance(output_file, str): raise Exception('output_file (%s) must be a string.' % str(output_file)) @@ -5308,26 +5313,30 @@ def _check_refine_binary_events(events, companions, if n_proc > 1 and multi_proc == False: raise Exception('if multi_proc is False, n_proc must = 1') - # Check to see that the filter name, photometric system, red_law are valid - if photometric_system not in photometric_system_dict: - exception_str = 'photometric_system must be a key in ' \ - 'photometric_system_dict. \n' \ - 'Acceptable values are : ' - for photometric_system in photometric_system_dict: - exception_str += '%s, ' % photometric_system - exception_str = exception_str[:-2] - raise Exception(exception_str) + # Check to see that the filter name and photometric system in filter_dict and red_law are valid + for system in filter_dict: + if system not in photometric_system_dict: + exception_str = 'photometric_system must be a key in ' \ + 'photometric_system_dict. \n' \ + 'Acceptable values are : ' + for photometric_system in photometric_system_dict: + exception_str += '%s, ' % photometric_system + exception_str = exception_str[:-2] + raise Exception(exception_str) - if filter_name not in photometric_system_dict[photometric_system]: - exception_str = 'filter_name must be a value in ' \ - 'photometric_system_dict[%s]. \n' \ - 'Acceptable values are : ' % photometric_system - for filter_name in photometric_system_dict[photometric_system]: - exception_str += '%s, ' % filter_name - exception_str = exception_str[:-2] - raise Exception(exception_str) + for filt in filter_dict[system]: + if filt not in photometric_system_dict[system]: + exception_str = 'filter_name must be a value in ' \ + 'photometric_system_dict[%s]. \n' \ + 'Acceptable values are : ' % system + for filter_name in photometric_system_dict[system]: + exception_str += '%s, ' % filter_name + exception_str = exception_str[:-2] + raise Exception(exception_str) + + key = system + '_' + filt -def refine_binary_events(events, companions, photometric_system, filter_name, +def refine_binary_events(events, companions, filter_dict, red_law = 'Damineli16', overwrite = False, output_file = 'default', save_phot = False, phot_dir = None, @@ -5399,6 +5408,9 @@ def refine_binary_events(events, companions, photometric_system, filter_name, """ start_time = time.time() + photometric_system = list(filter_dict.keys())[0] + filter_name = filter_dict[photometric_system][0] + if not overwrite and os.path.isfile(output_file): raise Exception('That refined_events.fits file name is taken! ' 'Either delete the .fits file, or pick a new name.') @@ -5408,7 +5420,7 @@ def refine_binary_events(events, companions, photometric_system, filter_name, # Error handling/complaining if input types are not right. _check_refine_binary_events(events, companions, - photometric_system, filter_name, + filter_dict, overwrite, output_file, save_phot, phot_dir, n_proc, multi_proc) @@ -5477,7 +5489,7 @@ def refine_binary_events(events, companions, photometric_system, filter_name, obj_id_L_S = (obj_id_L, obj_id_S) event_table_df['companion_idx_list'].loc[obj_id_L_S] = list(grouped_comps.groups[i]['companion_idx']) inputs[i] = [[event_table_df.loc[obj_id_L_S]], grouped_comps.groups[i].to_pandas(), obj_id_L, obj_id_S, - photometric_system, filter_name, red_law, save_phot, phot_dir, overwrite] + filter_dict, red_law, save_phot, phot_dir, overwrite] if multi_proc: results = pool.starmap(one_lightcurve_analysis, inputs) @@ -5596,7 +5608,7 @@ def refine_binary_events(events, companions, photometric_system, filter_name, return def one_lightcurve_analysis(event_table_row, comp_table_rows, obj_id_L, obj_id_S, - photometric_system, filter_name, red_law, + filter_dict, red_law, save_phot = False, phot_dir = None, overwrite = False): """ Generate BAGLE model, photometry, and generate binary lightcurve parameters. @@ -5683,8 +5695,7 @@ def one_lightcurve_analysis(event_table_row, comp_table_rows, obj_id_L, obj_id_S for comp_idx_L in comp_idxs_L: global_comp_idx_L = comp_table_rows['companion_idx'][comp_idx_L] name = "L_{}_S_{}".format(obj_id_L, obj_id_S) + "compL_{}_compS_{}".format(global_comp_idx_L, global_comp_idx_S) - model_parameter_dict, _, _, model_name = lightcurves.get_bsbl_lightcurve_parameters(event_table_row, comp_table_rows, int(comp_idx_L), int(comp_idx_S), - photometric_system, filter_name, red_law, event_id = 0) + model_parameter_dict, _, _, model_name = lightcurves.get_bsbl_lightcurve_parameters(event_table_row, comp_table_rows, int(comp_idx_L), int(comp_idx_S), filter_dict, red_law, event_id = 0) mod_class = getattr(model, model_name) try: mod = mod_class(**model_parameter_dict) @@ -5708,7 +5719,7 @@ def one_lightcurve_analysis(event_table_row, comp_table_rows, obj_id_L, obj_id_S global_comp_idx = comp_table_rows['companion_idx'][comp_idx] name = "L_{}_S_{}".format(obj_id_L, obj_id_S) + "compL_{}".format(global_comp_idx) model_parameter_dict, _, _, model_name = lightcurves.get_psbl_lightcurve_parameters(event_table_row, comp_table_rows, comp_idx, - photometric_system, filter_name, event_id = 0) + filter_dict, event_id = 0) mod_class = getattr(model, model_name) try: mod = mod_class(**model_parameter_dict) @@ -5730,7 +5741,7 @@ def one_lightcurve_analysis(event_table_row, comp_table_rows, obj_id_L, obj_id_S global_comp_idx = comp_table_rows['companion_idx'][comp_idx] name = "L_{}_S_{}".format(obj_id_L, obj_id_S) + "compS_{}".format(global_comp_idx) model_parameter_dict, _, _, model_name = lightcurves.get_bspl_lightcurve_parameters(event_table_row, comp_table_rows, comp_idx, - photometric_system, filter_name, red_law, event_id = 0) + filter_dict, red_law, event_id = 0) mod_class = getattr(model, model_name) try: mod = mod_class(**model_parameter_dict) @@ -5789,13 +5800,13 @@ def model_param_dict2fits_header(model_parameter_dict, phot_dir, name): """ fits_file = phot_dir + '/' + name + '_phot.fits' + + param_table = Table([model_parameter_dict]) + param_table_hdu = fits.table_to_hdu(param_table) + with fits.open(fits_file, 'update', memmap=False) as f: - hdr = f[0].header - for key in list(model_parameter_dict.keys()): - try: - hdr[key] = model_parameter_dict[key] - except ValueError: - hdr[key] = str(model_parameter_dict[key]) + f.append(param_table_hdu) + f.writeto(fits_file, overwrite=True) return @@ -5864,7 +5875,7 @@ def lightcurve_parameter_gen(model, model_parameter_dict, comp_idxs, obj_id_L, o # Handles the case of modeling a triple source as a binary source, # But it's CO + CO + star and you're modeling the CO + CO pair (so no flux) if 'mag_src_sec' in list(model_parameter_dict.keys()): - if np.isnan(model_parameter_dict['mag_src_sec']) and np.isnan(model_parameter_dict['mag_src_pri']): + if np.isnan(model_parameter_dict['mag_src_sec'][0]) and np.isnan(model_parameter_dict['mag_src_pri'][0]): return param_dict # These covers extreme cases that will never