From 1dc412bda44d6d9e05b8f1544b04fff5fddf6763 Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Mon, 2 Feb 2026 12:38:12 -0800 Subject: [PATCH 1/8] new function for monte carlo event estimate for popclass --- popsycle/analysis.py | 68 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index ac3871f..d3b4bc9 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -127,3 +127,71 @@ def count_stars_hdf5(hdf5_file, filt='ubv_I', mag_threshold=21): N_stars = len(gdx) return N_stars + +def events_for_popclass(h5_file, use_stars_per_bin=1e3): + """ + Draw microlensing events for random lens, source pairs from a + PopSyCLE singles-only catalog. + + Parameters + ---------- + hdf5_file : str + Filename of an hdf5 file. + + use_stars_per_bin : str + Order of magnitude number of stars per bin to use for the + calculation. Prevents excessive memory/computation use. + Default is 1e3. + + Returns + ------- + thetaEs : np.array + array of Einstein ring radii for events (mas) + piEs : np.array + array of microlensing parallaxes for events (unitless) + tEs : np.array + array of timescales for events (days) + weights : np.array + array of relative weights for events (mu_rel * thetaE) + """ + hf = h5py.File(h5_file, 'r') + thetaEs = [] + piEs = [] + tEs = [] + weights = [] + for k in list(hf.keys()): + if '_' not in k: + print('running', k) + dat = hf[k] + + if dat.shape[0] > 0: + ds = int(np.floor(dat.shape[0]/use_stars_per_bin)) + patch = dat[::ds] + dists = patch['rad'] + mul = patch['mu_lcosb'] + mub = patch['mu_b'] + masses = patch['mass'] + idx = np.arange(len(masses)) + del patch + + src_idxs, lens_idxs = np.meshgrid(idx, idx) + src_idxs, lens_idxs = src_idxs.ravel(), lens_idxs.ravel() + + dist_comp = (dists[src_idxs] > dists[lens_idxs]) #source further than lens + use_srcs = src_idxs[dist_comp] + use_lens = lens_idxs[dist_comp] + + # Microlensing math + pi_rel = (1/dists[use_lens] - 1/dists[use_srcs]) + c, G, mSun, pctom = 299792458, 6.6743e-11, 1.98840987e+30, 3.08567758e+16 + theta_e = np.sqrt(4*G*mSun*masses[use_lens]*pi_rel/(1000*pctom*c**2)) * 180/np.pi * 60**2 * 1000 + pi_e = pi_rel / theta_e + mu_rel = np.sqrt((mul[use_lens]-mul[use_srcs])**2 + (mub[use_lens]-mub[use_srcs])**2) + t_e = theta_e/mu_rel * 365.25 # years -> days + thetamu = theta_e*mu_rel + thetaEs.append(theta_e) + piEs.append(pi_e) + tEs.append(t_e) + weights.append(thetamu) + print(f'Drew {len(np.concatenate(thetaEs))} events total') + return np.concatenate(thetaEs), np.concatenate(piEs), np.concatenate(tEs), np.concatenate(weights) From 48d8cfee56dcbc2f3c0d7ca10a1506277065a2dc Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Mon, 2 Feb 2026 13:43:42 -0800 Subject: [PATCH 2/8] make bins scale appropriately in monte carlo selection --- popsycle/analysis.py | 50 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 3 deletions(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index d3b4bc9..46d676d 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -128,7 +128,49 @@ def count_stars_hdf5(hdf5_file, filt='ubv_I', mag_threshold=21): return N_stars -def events_for_popclass(h5_file, use_stars_per_bin=1e3): +def count_stars_all(h5_file): + """ + Finds the number of stars or systems in the field. + + Parameters + ---------- + hdf5_file : str + Filename of an hdf5 file. + + Returns + ------- + n_stars : int + Number of stars. + """ + hf = h5py.File(h5_file, 'r') + n_stars = 0 + for k in list(hf.keys()): + if '_' not in k: + n_stars += hf[k].shape[0] + return n_stars + +def count_stars_all_per_bin(h5_file): + """ + Finds the number of stars or systems per bin. + + Parameters + ---------- + hdf5_file : str + Filename of an hdf5 file. + + Returns + ------- + n_stars : list + Number of stars per bin. + """ + hf = h5py.File(h5_file, 'r') + n_stars = [] + for k in list(hf.keys()): + if '_' not in k: + n_stars.append(hf[k].shape[0]) + return n_stars + +def events_for_popclass(h5_file, max_stars_per_bin=3e3): """ Draw microlensing events for random lens, source pairs from a PopSyCLE singles-only catalog. @@ -159,14 +201,15 @@ def events_for_popclass(h5_file, use_stars_per_bin=1e3): piEs = [] tEs = [] weights = [] + stars_per_bin = count_stars_all_per_bin(h5_file) + scale_stars_used = np.maximum(1,int(np.floor(np.max(stars_per_bin)/max_stars_per_bin))) for k in list(hf.keys()): if '_' not in k: print('running', k) dat = hf[k] if dat.shape[0] > 0: - ds = int(np.floor(dat.shape[0]/use_stars_per_bin)) - patch = dat[::ds] + patch = dat[::scale_stars_used] dists = patch['rad'] mul = patch['mu_lcosb'] mub = patch['mu_b'] @@ -195,3 +238,4 @@ def events_for_popclass(h5_file, use_stars_per_bin=1e3): weights.append(thetamu) print(f'Drew {len(np.concatenate(thetaEs))} events total') return np.concatenate(thetaEs), np.concatenate(piEs), np.concatenate(tEs), np.concatenate(weights) + From 6d43cdc32da37b5b901b82a0c90c1ff30ec32e8c Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Mon, 2 Feb 2026 13:50:04 -0800 Subject: [PATCH 3/8] update comments --- popsycle/analysis.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index 46d676d..477436c 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -180,10 +180,10 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): hdf5_file : str Filename of an hdf5 file. - use_stars_per_bin : str - Order of magnitude number of stars per bin to use for the + max_stars_per_bin : str + Maxinum number of stars per bin to use for the calculation. Prevents excessive memory/computation use. - Default is 1e3. + Default is 3e3. Returns ------- From 412787e74fc24d54ecb4c950ad6b571a131d789b Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Fri, 6 Feb 2026 15:56:59 -0800 Subject: [PATCH 4/8] actually include lens type in popclass shortcut function --- popsycle/analysis.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index 477436c..a0b6e39 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -187,6 +187,8 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): Returns ------- + rem_ids : np.array + array of lens types thetaEs : np.array array of Einstein ring radii for events (mas) piEs : np.array @@ -197,6 +199,7 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): array of relative weights for events (mu_rel * thetaE) """ hf = h5py.File(h5_file, 'r') + rem_ids = [] thetaEs = [] piEs = [] tEs = [] @@ -214,6 +217,7 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): mul = patch['mu_lcosb'] mub = patch['mu_b'] masses = patch['mass'] + rem_id_catalog = patch['rem_id'] idx = np.arange(len(masses)) del patch @@ -223,6 +227,7 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): dist_comp = (dists[src_idxs] > dists[lens_idxs]) #source further than lens use_srcs = src_idxs[dist_comp] use_lens = lens_idxs[dist_comp] + rem_id = rem_id_catalog[use_lens] # Microlensing math pi_rel = (1/dists[use_lens] - 1/dists[use_srcs]) @@ -232,10 +237,11 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): mu_rel = np.sqrt((mul[use_lens]-mul[use_srcs])**2 + (mub[use_lens]-mub[use_srcs])**2) t_e = theta_e/mu_rel * 365.25 # years -> days thetamu = theta_e*mu_rel + rem_ids.append(rem_id) thetaEs.append(theta_e) piEs.append(pi_e) tEs.append(t_e) weights.append(thetamu) print(f'Drew {len(np.concatenate(thetaEs))} events total') - return np.concatenate(thetaEs), np.concatenate(piEs), np.concatenate(tEs), np.concatenate(weights) + return np.concatenate(rem_ids), np.concatenate(thetaEs), np.concatenate(piEs), np.concatenate(tEs), np.concatenate(weights) From efcd100fb86949c7088dd905c7d507a1798ca443 Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Mon, 6 Apr 2026 13:58:36 -0700 Subject: [PATCH 5/8] new blending function in analysis.py --- popsycle/analysis.py | 156 ++++++++++++++++++++++++++++++++++++++++-- popsycle/synthetic.py | 6 +- 2 files changed, 152 insertions(+), 10 deletions(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index a0b6e39..6bfd15e 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -2,11 +2,15 @@ import numpy as np import h5py import pylab as plt - -from popsycle import utils -from popsycle import synthetic +from scipy.spatial import KDTree +from popsycle import utils, synthetic, phot_utils import time import os +import pdb +import warnings +from itertools import zip_longest + +filt_dict = phot_utils.make_filt_dict() def get_star_system_pos_mag(hdf5_file, filt='ubv_I', recalc=True): """ @@ -57,7 +61,7 @@ def get_star_system_pos_mag(hdf5_file, filt='ubv_I', recalc=True): list_of_df = [] # Loop through the fields and aggregate the stars. - print('Countint stars in patches: ', end="") + print('Counting stars in patches: ', end="") for k in field_keys: print(f'{k}, ', end="") patch = np.array(hf[k]) @@ -145,7 +149,7 @@ def count_stars_all(h5_file): hf = h5py.File(h5_file, 'r') n_stars = 0 for k in list(hf.keys()): - if '_' not in k: + if (k.startswith('l')) & ('_' not in k): n_stars += hf[k].shape[0] return n_stars @@ -166,7 +170,7 @@ def count_stars_all_per_bin(h5_file): hf = h5py.File(h5_file, 'r') n_stars = [] for k in list(hf.keys()): - if '_' not in k: + if (k.startswith('l')) & ('_' not in k): n_stars.append(hf[k].shape[0]) return n_stars @@ -207,7 +211,7 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): stars_per_bin = count_stars_all_per_bin(h5_file) scale_stars_used = np.maximum(1,int(np.floor(np.max(stars_per_bin)/max_stars_per_bin))) for k in list(hf.keys()): - if '_' not in k: + if (k.startswith('l')) & ('_' not in k): print('running', k) dat = hf[k] @@ -245,3 +249,141 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): print(f'Drew {len(np.concatenate(thetaEs))} events total') return np.concatenate(rem_ids), np.concatenate(thetaEs), np.concatenate(piEs), np.concatenate(tEs), np.concatenate(weights) + +# For the blend calculation, sum the mags +def blend_magsum(magss): + return -2.5*np.log10(np.sum(10**(-0.4*magss),axis=1)) + +# For the blend calculation, flux-weight average the astrometric values +def blend_mag_weighted_sum(mags, vals): + fluxes = 10**(-0.4*mags) + sum_vals = (np.sum(fluxes.T*vals.T, axis=1)/np.sum(fluxes,axis=1)).T + return sum_vals + +# For the blend calculation, use a KDTree to assign stars to blend "clusters" +def calc_blend_clusters(points, radius, star_idxs): + # Set up KDTree and array to note whether each point has been assigned a cluster + tree = KDTree(points) + n_points = len(points) + assigned = np.zeros(n_points, dtype=bool) + clusters = [] + prim_idxs = [] + for i in range(n_points): + # Find unassigned neighbors and save cluster + if not assigned[i]: + neighbor_indices = tree.query_ball_point(points[i], r=radius) + current_cluster = [idx for idx in neighbor_indices if not assigned[idx]] + assigned[current_cluster] = True + clusters.append(current_cluster) + prim_idxs.append(star_idxs[i]) + return clusters, np.array(prim_idxs) + +""" +Take a resolved star catalog and blend it, summing magnitudes +and calculating flux-weighted positions, proper +motions, and parallaxes. +""" +def calc_blends_bin(star_dat, blend_rad, filters, + primary_filter=None, ext_law='Damineli16'): + # Clean up table and sort by magnitude + if primary_filter is None: + warnings.warn(f"No primary_filter provided. Using {filters[0]}"+ + " to sort stars and flux-weight additional parameters.") + primary_filter = filters[0] + for filt in filters: + star_dat.loc[:,filt] = synthetic.calc_app_mag(star_dat['rad'].to_numpy(), + star_dat[filt].to_numpy(), star_dat['exbv'].to_numpy(), + filt_dict[filt][ext_law]) + star_dat = star_dat[~np.isnan(star_dat[primary_filter])].copy() + star_dat.sort_values(by=primary_filter, inplace=True) + star_dat.reset_index(drop=True,inplace=True) + + # Get coordinates into simplified delta_l_cosb and delta_b frame + l_mean = np.mean(star_dat['glon']); b_mean = np.mean(star_dat['glat']) + all_l, all_b = star_dat['glon'].to_numpy(), star_dat['glat'].to_numpy() + delta_l_cosb = ((all_l-360*(all_l>180))-l_mean)*np.cos(all_b*np.pi/180) + delta_b = all_b-b_mean + all_pts = np.transpose([all_l,all_b]) + star_idxs = star_dat['obj_id'].to_numpy() + + # Set up arrays with necessary data + mags = star_dat[filters].to_numpy() + mags_prim = star_dat[primary_filter].to_numpy() + star_dat.loc[:,'plx'] = 1/star_dat['rad'] + astrom_vals = star_dat[['glon','glat','mu_lcosb', + 'mu_b','plx','exbv']].to_numpy() + blend_vals = [] + + # Run the cluster calculation on the sorted/cleaned points + clusters0, prim_idxs = calc_blend_clusters(all_pts, blend_rad, star_idxs) + # Turn this into a masked array for easy mag + param sums + clusters_arr = np.array(list(zip_longest(*clusters0, fillvalue=-1))).T + clusters = np.ma.masked_values(clusters_arr, -1) + + # Compute the magnitude sums & other values if desired + mags = np.append(mags, [np.repeat(np.inf, len(filters))], axis=0) + bmags = blend_magsum(mags[clusters,:]) + out = {} + for i,f in enumerate(filters): + out[f] = bmags[:,i] + mags_prim = np.append(mags_prim, [np.inf]) + astrom_vals = np.append(astrom_vals, [np.repeat(0, 6)], axis=0) + bvals = blend_mag_weighted_sum(mags_prim[clusters], astrom_vals[clusters,:]) + for i, col in enumerate(['glon','glat','mu_lcosb', + 'mu_b','plx', 'exbv']): + out[col] = bvals[:,i] + out['obj_id'] = prim_idxs + print(out.keys()) + return out + +def calc_blends(hdf5_file, blend_rad, filters, + primary_filter=None, recalc=False): + """ + Artificially blend the catalog based on a given radius & + save the blend catalog to a file. + + Parameters + ---------- + hdf5_file : str + Filename of an hdf5 file. + blend_radius : float + blend radius in degrees + filters : list of str + column headers for photometry to blend + primary_filter : str = None + filter to use for astrometric parameter flux-weighting. + if None, use the first filter in the list. + recalc : boolean = False + if True, recalculate and replace existing file + """ + outfile_name = hdf5_file.replace('.h5', '_blended_catalog.h5') + + if os.path.exists(outfile_name) and recalc == False: + print(f'{outfile_name} exists. Set recalc=True to recalculate.') + return + + hf = h5py.File(hdf5_file, 'r') + outfile = h5py.File(outfile_name, 'a') + + for dset_name in list(hf.keys()): + if (dset_name.startswith('l')) & ('_' not in dset_name): + print(dset_name) + blend_bin = calc_blends_bin(pd.DataFrame(hf[dset_name][:]), + blend_rad, filters, + primary_filter=primary_filter) + compound_dtype = synthetic._generate_compound_dtype(blend_bin) + save_data = np.empty(len(blend_bin['plx']), dtype=compound_dtype) + for colname in blend_bin: + save_data[colname] = blend_bin[colname] + if dset_name not in outfile: + dataset = outfile.create_dataset(dset_name, shape=(0,), + chunks=(1e4,), + maxshape=(None,), + dtype=compound_dtype) + else: + dataset = outfile[dset_name] + new_size = len(blend_bin['plx']) + dataset.resize((new_size, )) + dataset[:] = save_data + hf.close() + outfile.close() diff --git a/popsycle/synthetic.py b/popsycle/synthetic.py index 2bcd0fa..2862ea1 100755 --- a/popsycle/synthetic.py +++ b/popsycle/synthetic.py @@ -6406,11 +6406,11 @@ def calc_app_mag(r, M, E, f): Parameters ---------- - M : float or array - Absolute magnitude of star - r : float or array Distance of star from sun (in kpc) + + M : float or array + Absolute magnitude of star E : float or array Extinction law From 64c731a458f59f4166bd2da3d4cd4c086d4c744f Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Mon, 6 Apr 2026 14:35:52 -0700 Subject: [PATCH 6/8] clean up + improve flexibility of calc_blends --- popsycle/analysis.py | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index 6bfd15e..35e6012 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -337,7 +337,8 @@ def calc_blends_bin(star_dat, blend_rad, filters, return out def calc_blends(hdf5_file, blend_rad, filters, - primary_filter=None, recalc=False): + primary_filter=None, ext_law='Damineli16', + recalc=False, combine_bins=False): """ Artificially blend the catalog based on a given radius & save the blend catalog to a file. @@ -353,8 +354,12 @@ def calc_blends(hdf5_file, blend_rad, filters, primary_filter : str = None filter to use for astrometric parameter flux-weighting. if None, use the first filter in the list. + ext_law : str = 'Damineli16' + extinction law recalc : boolean = False if True, recalculate and replace existing file + combine_bins : boolean = False + if True, save a single dataset instead of binned data """ outfile_name = hdf5_file.replace('.h5', '_blended_catalog.h5') @@ -370,20 +375,31 @@ def calc_blends(hdf5_file, blend_rad, filters, print(dset_name) blend_bin = calc_blends_bin(pd.DataFrame(hf[dset_name][:]), blend_rad, filters, - primary_filter=primary_filter) + primary_filter=primary_filter, ext_law=ext_law) compound_dtype = synthetic._generate_compound_dtype(blend_bin) save_data = np.empty(len(blend_bin['plx']), dtype=compound_dtype) for colname in blend_bin: save_data[colname] = blend_bin[colname] - if dset_name not in outfile: - dataset = outfile.create_dataset(dset_name, shape=(0,), - chunks=(1e4,), - maxshape=(None,), - dtype=compound_dtype) + if combine_bins: + if 'data' not in outfile: + dataset = outfile.create_dataset('data', shape=(0,), + chunks=(1e4,), + maxshape=(None,), + dtype=compound_dtype) + else: + dataset = outfile['data'] + old_size = dataset.shape[0] else: - dataset = outfile[dset_name] - new_size = len(blend_bin['plx']) + if dset_name not in outfile: + dataset = outfile.create_dataset(dset_name, shape=(0,), + chunks=(1e4,), + maxshape=(None,), + dtype=compound_dtype) + else: + dataset = outfile[dset_name] + old_size = 0 + new_size = old_size + len(blend_bin['plx']) dataset.resize((new_size, )) - dataset[:] = save_data + dataset[old_size:new_size] = save_data hf.close() outfile.close() From d716975d1c596c1158af53d46405a9b31633d337 Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Mon, 6 Apr 2026 14:46:50 -0700 Subject: [PATCH 7/8] clean up + document calc_blends --- popsycle/analysis.py | 42 ++++++++++++++++-------------------------- 1 file changed, 16 insertions(+), 26 deletions(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index 35e6012..cf82ccb 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -250,18 +250,10 @@ def events_for_popclass(h5_file, max_stars_per_bin=3e3): return np.concatenate(rem_ids), np.concatenate(thetaEs), np.concatenate(piEs), np.concatenate(tEs), np.concatenate(weights) -# For the blend calculation, sum the mags -def blend_magsum(magss): - return -2.5*np.log10(np.sum(10**(-0.4*magss),axis=1)) - -# For the blend calculation, flux-weight average the astrometric values -def blend_mag_weighted_sum(mags, vals): - fluxes = 10**(-0.4*mags) - sum_vals = (np.sum(fluxes.T*vals.T, axis=1)/np.sum(fluxes,axis=1)).T - return sum_vals - -# For the blend calculation, use a KDTree to assign stars to blend "clusters" -def calc_blend_clusters(points, radius, star_idxs): +def _calc_blend_clusters(points, radius, star_idxs): + """ + Helper function for calc_blends to run the KDTree clustering + """ # Set up KDTree and array to note whether each point has been assigned a cluster tree = KDTree(points) n_points = len(points) @@ -278,13 +270,11 @@ def calc_blend_clusters(points, radius, star_idxs): prim_idxs.append(star_idxs[i]) return clusters, np.array(prim_idxs) -""" -Take a resolved star catalog and blend it, summing magnitudes -and calculating flux-weighted positions, proper -motions, and parallaxes. -""" -def calc_blends_bin(star_dat, blend_rad, filters, +def _calc_blends_bin(star_dat, blend_rad, filters, primary_filter=None, ext_law='Damineli16'): + """ + Helper function for calc_blends to process 1 bin + """ # Clean up table and sort by magnitude if primary_filter is None: warnings.warn(f"No primary_filter provided. Using {filters[0]}"+ @@ -303,7 +293,7 @@ def calc_blends_bin(star_dat, blend_rad, filters, all_l, all_b = star_dat['glon'].to_numpy(), star_dat['glat'].to_numpy() delta_l_cosb = ((all_l-360*(all_l>180))-l_mean)*np.cos(all_b*np.pi/180) delta_b = all_b-b_mean - all_pts = np.transpose([all_l,all_b]) + all_pts = np.transpose([delta_l_cosb, delta_b]) star_idxs = star_dat['obj_id'].to_numpy() # Set up arrays with necessary data @@ -312,28 +302,28 @@ def calc_blends_bin(star_dat, blend_rad, filters, star_dat.loc[:,'plx'] = 1/star_dat['rad'] astrom_vals = star_dat[['glon','glat','mu_lcosb', 'mu_b','plx','exbv']].to_numpy() - blend_vals = [] # Run the cluster calculation on the sorted/cleaned points - clusters0, prim_idxs = calc_blend_clusters(all_pts, blend_rad, star_idxs) + clusters0, prim_idxs = _calc_blend_clusters(all_pts, blend_rad, star_idxs) # Turn this into a masked array for easy mag + param sums clusters_arr = np.array(list(zip_longest(*clusters0, fillvalue=-1))).T clusters = np.ma.masked_values(clusters_arr, -1) - # Compute the magnitude sums & other values if desired + # Compute the magnitude sums & weighted astrometry parameters mags = np.append(mags, [np.repeat(np.inf, len(filters))], axis=0) - bmags = blend_magsum(mags[clusters,:]) + bmags = -2.5*np.log10(np.sum(10**(-0.4*mags[clusters,:]),axis=1)) out = {} for i,f in enumerate(filters): out[f] = bmags[:,i] mags_prim = np.append(mags_prim, [np.inf]) astrom_vals = np.append(astrom_vals, [np.repeat(0, 6)], axis=0) - bvals = blend_mag_weighted_sum(mags_prim[clusters], astrom_vals[clusters,:]) + fluxes = 10**(-0.4*mags_prim[clusters]) + bvals = (np.sum(fluxes.T*astrom_vals[clusters,:].T, axis=1) + / np.sum(fluxes,axis=1)).T for i, col in enumerate(['glon','glat','mu_lcosb', 'mu_b','plx', 'exbv']): out[col] = bvals[:,i] out['obj_id'] = prim_idxs - print(out.keys()) return out def calc_blends(hdf5_file, blend_rad, filters, @@ -373,7 +363,7 @@ def calc_blends(hdf5_file, blend_rad, filters, for dset_name in list(hf.keys()): if (dset_name.startswith('l')) & ('_' not in dset_name): print(dset_name) - blend_bin = calc_blends_bin(pd.DataFrame(hf[dset_name][:]), + blend_bin = _calc_blends_bin(pd.DataFrame(hf[dset_name][:]), blend_rad, filters, primary_filter=primary_filter, ext_law=ext_law) compound_dtype = synthetic._generate_compound_dtype(blend_bin) From f6ce9c96a1ce5f5c1d1dd228ad5244163c5d8b54 Mon Sep 17 00:00:00 2001 From: Macy Huston Date: Tue, 7 Apr 2026 15:19:19 -0700 Subject: [PATCH 8/8] improve documentation --- popsycle/analysis.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/popsycle/analysis.py b/popsycle/analysis.py index cf82ccb..2d4d695 100644 --- a/popsycle/analysis.py +++ b/popsycle/analysis.py @@ -331,7 +331,10 @@ def calc_blends(hdf5_file, blend_rad, filters, recalc=False, combine_bins=False): """ Artificially blend the catalog based on a given radius & - save the blend catalog to a file. + save the blend catalog to a file. Outputs are: summed magnitudes + in selected filters, primary filter flux-weighted positions, + proper motions, parallaxes, and reddening, and the obj_id in + the original table of the brightest star in the blend. Parameters ----------