From 07a525f2220056dffd87d85f3a6cf6b705fdff17 Mon Sep 17 00:00:00 2001 From: Zstone19 Date: Wed, 19 Aug 2026 17:15:13 -0500 Subject: [PATCH 1/3] Allow user to specify number of threads. Put prior on delta_i when delay_dist=True of [0, baseline*5], where baseline is the maximum of the reference+line baselines. Improves performance significantly when delta_i cannot blow up. --- PyROA/PyROA.py | 44 ++++++++++++++++++++++++++++---------------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/PyROA/PyROA.py b/PyROA/PyROA.py index 6c989f8..d99b3ab 100644 --- a/PyROA/PyROA.py +++ b/PyROA/PyROA.py @@ -1021,9 +1021,12 @@ def log_prior(params, priors, add_var, data, delay_dist, AccDisc, wavelengths, i if (add_var == True): V = params_chunks[i][-1]/(init_params_chunks[i][-1]*5) - - if (delay_dist == True and i>0): - if (params_chunks[i][3]>=0.0): + + if (delay_dist == True and i>0): #########ADDED BY ZS - STOPS DELTA_i FROM BLOWING UP (08/19/2026)######## + baseline_0 = np.nanmax(data[0][:,0]) - np.nanmin(data[0][:,0]) + baseline_i = np.nanmax(data[i][:,0]) - np.nanmin(data[i][:,0]) + maxval = max(baseline_0, baseline_i)*5 + if (params_chunks[i][3]>=0.0) and (params_chunks[i][3] <= maxval): tau_rms = params_chunks[i][3] #pr.append(2.0*np.log((1.0/np.sqrt(2.0*np.pi*(rms_prior_width**2)))*np.exp(-0.5*(tau_rms/rms_prior_width)**2))) check.append(0.0) @@ -1138,7 +1141,7 @@ def Slow(t, S0, dS, t0): def FullFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, Nburnin, include_slow_comp, slow_comp_delta, calc_P, delay_dist, psi_types, pos_ref, AccDisc, wavelengths, filters, use_backend, - resume_progress, plot_corner,memfunction, gridsize): + resume_progress, plot_corner,memfunction, gridsize, nthread=1): Nchunk = 2 @@ -1384,7 +1387,8 @@ def FullFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, #print('int(Npar - param_delete)', int(Npar - param_delete)) - pos = [pos_min + psize*np.random.rand(int(Npar - param_delete)) for i in range(2*((Npar-param_delete)))] + # pos = [pos_min + psize*np.random.rand(int(Npar - param_delete)) for i in range(2*((Npar-param_delete)))] + pos = 0.2*np.array(pos)* np.random.randn(int(2.0*Npar), int(Npar - param_delete)) + np.array(pos) + 1e-4* np.random.randn(int(2.0*Npar), int(Npar - param_delete)) pos = np.array(pos) #print(np.array(pos)) @@ -1400,7 +1404,8 @@ def FullFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, else: backend = None - with Pool() as pool: + ########ADDED BY ZS - ALLOWS THE USER TO CHOOSE THE NUMBER OF THREADS USED######## + with Pool(nthread) as pool: sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=[data, priors, add_var, size,sig_level, include_slow_comp, slow_comp_delta, P_func, slow_comps, P_slow, init_delta, delay_dist, psi_types, @@ -1667,9 +1672,12 @@ def FullFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, mx=max(merged_mjd) mn=min(merged_mjd) length = abs(mx-mn) - t = np.arange(mn, mx, length/(gridsize)) - + if gridsize is not None: + t = np.arange(mn, mx, length/(gridsize)) + else: + t = np.arange(mn, mx, length/1000) + ts, Xs, errss = RunningOptimalAverageOutConv(t, merged_mjd, merged_flux, merged_err, factors, conv, prev, x, delta_new) @@ -1814,13 +1822,15 @@ def __init__(self, datadir, objName, filters, priors, delay_ref = None, init_tau delay_dist=False , psi_types = None, add_var=True, sig_level = 4.0, Nsamples=10000, Nburnin=0, include_slow_comp=False, slow_comp_delta=30.0, calc_P=False, AccDisc=False, wavelengths=None, - use_backend = False, resume_progress = False, plot_corner=False,memfunction='gaussian', gridsize = None): + use_backend = False, resume_progress = False, plot_corner=False,memfunction='gaussian', gridsize = None, + nthread=1): if datadir[-1] != '/': datadir += '/' #Add forward slash in case it isn't there self.datadir=datadir self.objName=objName self.filters=filters self.gridsize = gridsize + self.nthread = nthread data=[] for i in range(len(filters)): file = datadir + str(self.objName) +"_"+ str(self.filters[i]) + ".dat" @@ -1888,7 +1898,8 @@ def __init__(self, datadir, objName, filters, priors, delay_ref = None, init_tau self.sig_level, self.Nsamples, self.Nburnin, self.include_slow_comp, self.slow_comp_delta, self.calc_P, self.delay_dist, self.psi_types, self.delay_ref_pos, self.AccDisc, self.wavelengths, self.filters, - self.use_backend, self.resume_progress,plot_corner,memfunction, self.gridsize) + self.use_backend, self.resume_progress,plot_corner,memfunction, self.gridsize, + nthread=self.nthread) self.samples = run[0] self.samples_flat = run[1] @@ -2356,7 +2367,7 @@ def log_probability2(params, data, priors, sig_level, init_params_chunks,memfunc -def InterCalib(data, priors, init_delta, sig_level, Nsamples, Nburnin, filter,plot_corner,memfunction, gridsize): +def InterCalib(data, priors, init_delta, sig_level, Nsamples, Nburnin, filter,plot_corner,memfunction, gridsize, nthread=1): ######################################################################################## #Run MCMC to fit to data @@ -2426,7 +2437,7 @@ def InterCalib(data, priors, init_delta, sig_level, Nsamples, Nburnin, filter,pl np.savetxt('test_initial_points.txt',pos.T) nwalkers, ndim = pos.shape - with Pool() as pool: + with Pool(nthread) as pool: sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability2, args=(data, priors, sig_level, init_params_chunks,memfunction, gridsize), pool=pool) @@ -2564,11 +2575,12 @@ def InterCalib(data, priors, init_delta, sig_level, Nsamples, Nburnin, filter,pl class InterCalibrate(): def __init__(self, datadir, objName, filter, scopes, priors, init_delta=1.0, sig_level = 3.0, - Nsamples=15000, Nburnin=10000,plot_corner=False,memfunction='gaussian'): + Nsamples=15000, Nburnin=10000,plot_corner=False,memfunction='gaussian', nthread=1): self.datadir=datadir self.objName=objName self.filter=filter self.scopes=scopes + self.nthread = nthread scopes_array = [] data=[] for i in range(len(scopes)): @@ -2590,7 +2602,7 @@ def __init__(self, datadir, objName, filter, scopes, priors, init_delta=1.0, sig run = InterCalib(data, self.priors, self.init_delta, self.sig_level, self.Nsamples, - self.Nburnin, self.filter,plot_corner,memfunction) + self.Nburnin, self.filter,plot_corner,memfunction, nthread=self.nthread) self.samples = run[0] self.samples_flat = run[1] @@ -2958,7 +2970,7 @@ def log_probability3(params, data, priors, add_var, size, sig_level): -def LensFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, Nburnin, image, file): +def LensFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, Nburnin, image, file, nthread=1): if (add_var == True): Npar = 7*len(data) + 3 @@ -3085,7 +3097,7 @@ def LensFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, Nb nwalkers, ndim = pos.shape print("Nwalkers = ", nwalkers) - with Pool() as pool: + with Pool(nthread) as pool: sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability3, args=(data, priors, add_var, size, sig_level), pool=pool) sampler.run_mcmc(pos, Nsamples, progress=True); From ab848223127d59782d0aaf70102a814b2446bab5 Mon Sep 17 00:00:00 2001 From: Zstone19 Date: Tue, 25 Aug 2026 13:18:01 -0500 Subject: [PATCH 2/3] Changed starting position back to normal --- PyROA/PyROA.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/PyROA/PyROA.py b/PyROA/PyROA.py index d99b3ab..7a1eb7e 100644 --- a/PyROA/PyROA.py +++ b/PyROA/PyROA.py @@ -1022,7 +1022,7 @@ def log_prior(params, priors, add_var, data, delay_dist, AccDisc, wavelengths, i V = params_chunks[i][-1]/(init_params_chunks[i][-1]*5) - if (delay_dist == True and i>0): #########ADDED BY ZS - STOPS DELTA_i FROM BLOWING UP (08/19/2026)######## + if (delay_dist == True and i>0): #Upper limit on delta_i stops it from blowing up in MCMC baseline_0 = np.nanmax(data[0][:,0]) - np.nanmin(data[0][:,0]) baseline_i = np.nanmax(data[i][:,0]) - np.nanmin(data[i][:,0]) maxval = max(baseline_0, baseline_i)*5 @@ -1387,8 +1387,7 @@ def FullFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, #print('int(Npar - param_delete)', int(Npar - param_delete)) - # pos = [pos_min + psize*np.random.rand(int(Npar - param_delete)) for i in range(2*((Npar-param_delete)))] - pos = 0.2*np.array(pos)* np.random.randn(int(2.0*Npar), int(Npar - param_delete)) + np.array(pos) + 1e-4* np.random.randn(int(2.0*Npar), int(Npar - param_delete)) + pos = [pos_min + psize*np.random.rand(int(Npar - param_delete)) for i in range(2*((Npar-param_delete)))] pos = np.array(pos) #print(np.array(pos)) From eda4ab97004be3b3d0c26a1e2e9cec8cd072b96b Mon Sep 17 00:00:00 2001 From: Zstone19 Date: Tue, 25 Aug 2026 13:23:39 -0500 Subject: [PATCH 3/3] Fixed comments --- PyROA/PyROA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PyROA/PyROA.py b/PyROA/PyROA.py index 7a1eb7e..53b6477 100644 --- a/PyROA/PyROA.py +++ b/PyROA/PyROA.py @@ -1403,7 +1403,7 @@ def FullFit(data, priors, init_tau, init_delta, add_var, sig_level, Nsamples, else: backend = None - ########ADDED BY ZS - ALLOWS THE USER TO CHOOSE THE NUMBER OF THREADS USED######## + with Pool(nthread) as pool: sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=[data, priors, add_var, size,sig_level, include_slow_comp,