diff --git a/ares/analysis/GalaxyPopulation.py b/ares/analysis/GalaxyPopulation.py index b9ad9b758..c04d025e3 100644 --- a/ares/analysis/GalaxyPopulation.py +++ b/ares/analysis/GalaxyPopulation.py @@ -41,7 +41,8 @@ datasets_lf = ('bouwens2015', 'finkelstein2015', 'bowler2020', 'stefanon2019', 'mclure2013', 'parsa2016', 'atek2015', 'alavi2016', 'reddy2009', 'weisz2014', 'bouwens2017', 'oesch2018', 'oesch2013', - 'oesch2014', 'vanderburg2010', 'morishita2018', 'rojasruiz2020') + 'oesch2014', 'vanderburg2010', 'morishita2018', 'rojasruiz2020', + 'harikane2022') datasets_smf = ('song2016', 'stefanon2017', 'duncan2014', 'tomczak2014', 'moustakas2013', 'mortlock2011', 'marchesini2009_10', 'perez2008') datasets_mzr = ('sanders2015',) @@ -53,7 +54,7 @@ 'dropouts': ('parsa2016', 'bouwens2015', 'finkelstein2015', 'bowler2020','stefanon2019', 'mclure2013', 'vanderburg2010', 'reddy2009', 'oesch2018', 'oesch2013', 'oesch2014', - 'morishita2018', 'rojasruiz2020'), + 'morishita2018', 'rojasruiz2020', 'harikane2022'), 'lensing': ('alavi2016', 'atek2015', 'bouwens2017'), 'local': ('weisz2014,'), 'all': datasets_lf, diff --git a/ares/physics/HaloMassFunction.py b/ares/physics/HaloMassFunction.py index 0b593b247..f08692ce3 100755 --- a/ares/physics/HaloMassFunction.py +++ b/ares/physics/HaloMassFunction.py @@ -1421,7 +1421,7 @@ def tab_prefix_hmf(self, with_size=False): def SaveHMF(self, fn=None, clobber=False, destination=None, fmt='hdf5', save_MAR=True): - self.save(fn=fn, clobber=clobber, destination=destination, fm=fmt, + self.save(fn=fn, clobber=clobber, destination=destination, fmt=fmt, save_MAR=save_MAR) def save(self, fn=None, clobber=False, destination=None, fmt='hdf5', diff --git a/ares/populations/GalaxyEnsemble.py b/ares/populations/GalaxyEnsemble.py index 9b0363965..0aa53c012 100755 --- a/ares/populations/GalaxyEnsemble.py +++ b/ares/populations/GalaxyEnsemble.py @@ -272,6 +272,72 @@ def _cache_halos(self): def _cache_halos(self, value): self._cache_halos_ = value + def _get_target_halos(self, raw): + Nz = raw['z'].size + tvol = self.pf['pop_target_volume'] + tred = self.pf['pop_target_redshift'] + tseed = self.pf['pop_target_seed'] + + ## + # Modify raw histories if volume provided. + assert tred is not None + + # Here, each element corresponds to a distinct halo. + # Will add Poisson noise only. + + # Loop over bins, draw halos, then add noise etc. + iztred = np.argmin(np.abs(tred - raw['z'])) + + # Mean number of halos in each mass bin at target redshift + Nlam = raw['nh'][:,iztred] * tvol + + np.random.seed(seed=tseed) + + # Actual number in target volume + Nact = np.random.poisson(lam=Nlam, size=Nlam.size) + + Mmin = self.guide.Mmin(raw['z']) + + # Loop over mass bins, create halo histories. + # Need to create dictionary containing 'Mh', 'MAR', 'nh', etc., + # where each is 2-D array in (halo number, time) + ct = 0 + for i, Nbin in enumerate(Nact): + + # Skip Mh < Mmin for efficiency. + if np.all(raw['Mh'][i,iztred:] < Mmin[iztred:]): + continue + + if Nbin == 0: + continue + + _Mh = raw['Mh'][i,:] + _MAR = raw['MAR'][i,:] + + _MhT = np.tile(_Mh, (Nbin, 1)) + _MART = np.tile(_MAR, (Nbin, 1)) + + if ct == 0: + print("# Generating {} individual halos from {} mass bins...".format( + Nact[i:].sum(), Nact.size)) + + Mh = _MhT + MAR = _MART + else: + Mh = np.vstack((Mh, _MhT)) + MAR = np.vstack((MAR, _MART)) + + ct += 1 + + nh = np.ones(Mh.shape) + + out = raw.copy() + out['Mh'] = Mh + out['nh'] = nh + out['MAR'] = MAR + + return out + def _gen_halo_histories(self): """ From a set of smooth halo assembly histories, build a bigger set @@ -281,9 +347,13 @@ def _gen_halo_histories(self): if hasattr(self, '_cache_halos_'): return self._cache_halos + thin = self.pf['pop_thin_hist'] + raw = self.load() - thin = self.pf['pop_thin_hist'] + if self.pf['pop_target_volume'] is not None: + assert thin in [None, 0, 1, False] + raw = self._get_target_halos(raw) if thin > 1: assert self.pf['pop_histories'] is None, \ @@ -1031,10 +1101,16 @@ def _gen_prescribed_galaxy_histories(self, zstop=0): if 'SFR' in halos: SFR = halos['SFR'][:,-1::-1] else: - iz = np.argmin(np.abs(6. - z)) - SFR = self.guide.get_fstar(z=z2d, Mh=Mh) - np.multiply(SFR, MAR, out=SFR) - SFR *= fb + if self.pf['pop_sfr'] is not None: + SFR = np.ones_like(MAR) * self.pf['pop_sfr'] + sigma_mar = self.pf['pop_scatter_mar'] + np.random.seed(self.pf['pop_scatter_mar_seed']) + noise = self.get_noise_lognormal(SFR, sigma_mar) + SFR += noise + else: + SFR = self.guide.get_fstar(z=z2d, Mh=Mh) + np.multiply(SFR, MAR, out=SFR) + SFR *= fb ## # Duty cycle effects @@ -1105,12 +1181,21 @@ def _gen_prescribed_galaxy_histories(self, zstop=0): Mmin = self.halos.VirialMass(z2d, self.pf['pop_Tmin']) above_Mmin = Mh >= Mmin + # Cut out Mh < Mmin galaxies + if self.pf['pop_Mmax'] is not None: + above_Mmax = Mh >= self.pf['pop_Mmax'] + else: + above_Mmax = np.zeros_like(Mh) + # Bye bye guys SFR[above_Mmin==False] = 0 + SFR[above_Mmax==True] = 0 ## # Introduce some by-hand quenching. - if self.pf['pop_quench'] is not None: + qmethod = self.pf['pop_quench_method'] + #print('HELLO', qmethod, self.pf['pop_quench']) + if (self.pf['pop_quench'] is not None) and (qmethod == 'zreion'): zreion = self.pf['pop_quench'] if type(zreion) in [np.ndarray, np.ma.core.MaskedArray]: assert zreion.size == Mh.size, \ @@ -1166,6 +1251,23 @@ def _gen_prescribed_galaxy_histories(self, zstop=0): if sum(SFR[ok==1] == 0) < sum(is_quenched[ok==1]): raise ValueError(err) + elif (self.pf['pop_quench'] is not None) and \ + (qmethod in ['Mh', 'mass', 'sfr']): + + is_quenched = np.zeros_like(Mh) + + if qmethod in ['mass', 'Mh']: + is_quenched[Mh > self.pf['pop_quench']] = 1 + elif qmethod == 'sfr': + is_quenched[SFR > self.pf['pop_quench']] = 1 + else: + raise NotImplemented('help') + + # Bye bye guys + SFR[is_quenched==True] = 0 + + print('quenched', is_quenched.sum() / float(is_quenched.size)) + # Stellar mass should have zeros padded at the 0th time index Ms = np.hstack((zeros_like_Mh, @@ -1181,19 +1283,24 @@ def _gen_prescribed_galaxy_histories(self, zstop=0): ## if have_dust: - delay = self.pf['pop_dust_yield_delay'] + if self.pf['pop_dust_yield_delay'] not in [0, None]: + have_delay = True + else: + have_delay = False + + delay = self.guide.dust_yield_delay(z=z2d, Mh=Mh) if np.all(fg == 0): - if type(fd) in [int, float, np.float64] and delay == 0: + if type(fd) in [int, float, np.float64] and (not have_delay): Md = fd * fZy * Ms else: - if delay > 0: + if have_delay: assert np.allclose(np.diff(dt_myr), 0.0, rtol=1e-5, atol=1e-5) - shift = int(delay // dt_myr[0]) + shift = np.array(delay // dt_myr[0], dtype=int) # Need to fix so Mh-dep fd can still work. assert type(fd) in [int, float, np.float64] @@ -1361,6 +1468,11 @@ def _gen_prescribed_galaxy_histories(self, zstop=0): else: pos = None + # Can null dust reddening by hand if we want. + if np.isfinite(self.pf['pop_dust_kill_redshift']): + Sd[np.argwhere(z2d > self.pf['pop_dust_kill_redshift'])] = 0.0 + Md[np.argwhere(z2d > self.pf['pop_dust_kill_redshift'])] = 0.0 + del z2d # Pack up @@ -1578,6 +1690,36 @@ def get_field(self, z, field): iz = np.argmin(np.abs(z - self.histories['z'])) return self.histories[field][:,iz] + def get_ages(self, z, mass_frac=0.5): + """ + Get age of galaxies (in Myr), assumed to be time since they assembled + `mass_frac` times their mass at redshift `z`. + """ + + iz = np.argmin(np.abs(z - self.histories['z'])) + + Ms = self.get_field(z, 'Ms') + ages = np.zeros_like(Ms) + + tasc = self.histories['t'][-1::-1] + zasc = self.histories['z'][-1::-1] + for i, _z_ in enumerate(zasc): + if _z_ <= z: + continue + + Ms_z_ = self.get_field(_z_, 'Ms') + is_double = Ms_z_ < mass_frac * Ms + is_new = ages == 0 + + ok = is_double * is_new + dt = self.histories['t'][iz] - tasc[i] + ages[ok==1] = dt + + if np.all(ages > 0): + break + + return ages + def StellarMassFunction(self, z, bins=None, units='dex'): return self.get_smf(z, bins=bins, units=units) @@ -2800,11 +2942,20 @@ def get_lf(self, z, bins=None, use_mags=True, wave=1600., assert y[yok==1].min() < x.min() assert y[yok==1].max() > x.max() - hist, bin_histedges = np.histogram(y[yok==1], - weights=w[yok==1], bins=bin_c2e(x), density=True) + #if self.pf['pop_fobsc']: + #fobsc = (1. - self.guide.fobsc(z=z, Mh=self.halos.tab_M)) - N = np.sum(w[yok==1]) - phi = hist * N + if np.all(w[yok==1] == 1): + hist, bin_histedges = np.histogram(y[yok==1], + bins=bin_c2e(x), density=False) + dbin = bin_histedges[1] - bin_histedges[0] + phi = hist / self.pf['pop_target_volume'] / dbin + else: + hist, bin_histedges = np.histogram(y[yok==1], + weights=w[yok==1], bins=bin_c2e(x), density=True) + + N = np.sum(w[yok==1]) + phi = hist * N #self._cache_lf_[(z, wave)] = x, phi diff --git a/ares/util/SetDefaultParameterValues.py b/ares/util/SetDefaultParameterValues.py index 7a191129f..9b53e2ede 100644 --- a/ares/util/SetDefaultParameterValues.py +++ b/ares/util/SetDefaultParameterValues.py @@ -513,6 +513,11 @@ def PopulationParameters(): "pop_type": 'galaxy', + "pop_target_volume": None, + "pop_target_redshift": None, + "pop_target_density": 0., + "pop_target_seed": None, + "pop_tunnel": None, "pop_sfr_model": 'fcoll', # or sfrd-func, sfrd-tab, sfe-func, sfh-tab, rates, @@ -651,8 +656,6 @@ def PopulationParameters(): "pop_fstar_max": 1.0, "pop_fstar_negligible": 1e-5, # relative to maximum - "pop_sfr": None, - "pop_facc": 0.0, "pop_fsmooth": 1.0, @@ -677,7 +680,7 @@ def PopulationParameters(): "pop_aging": False, "pop_enrichment": False, "pop_quench": None, - "pop_quench_by": 'mass', + "pop_quench_method": 'zreion', "pop_flag_sSFR": None, "pop_flag_tau": None, diff --git a/input/litdata/furlanetto2017.py b/input/litdata/furlanetto2017.py index 88461a1a6..507974dea 100644 --- a/input/litdata/furlanetto2017.py +++ b/input/litdata/furlanetto2017.py @@ -6,14 +6,14 @@ 'pop_fstar': None, 'pop_fstar_max': 0.1, # fstar <= this value - + # SFE (through mass loading factor) 'pop_sfr_model': 'mlf-func', 'pop_mlf': 'pq[0]', 'pq_func[0]': 'pl_evolN', 'pq_func_var[0]': 'Mh', 'pq_func_var2[0]': '1+z', - + ## # Steve's Equation 13. ## @@ -22,8 +22,8 @@ 'pq_func_par2[0]': -2./3., 'pq_func_par3[0]': 9., 'pq_func_par4[0]': -1., - - 'pop_L1600_per_sfr': 1e-28, + + 'pop_lum_per_sfr': 1e-28, } @@ -44,7 +44,7 @@ 'pq_func[1]': 'pl_evolN', 'pq_func_var[1]': 'Mh', 'pq_func_var2[1]': '1+z', - + 'pq_val_ceil[1]': 1.0, # fshock <= 1 # Steve's Equation 6 (from Faucher-Giguere+ 2011) @@ -54,6 +54,3 @@ 'pq_func_par3[1]': 4., 'pq_func_par4[1]': 0.38, } - - - diff --git a/input/litdata/harikane2022.py b/input/litdata/harikane2022.py new file mode 100755 index 000000000..7d3250276 --- /dev/null +++ b/input/litdata/harikane2022.py @@ -0,0 +1,90 @@ +""" +Finkelstein et al., 2015, ApJ, 810, 71 +""" + +import numpy as np + +info = \ +{ + 'reference': 'Harikane et al., 2022', + 'data': 'Table 7', + 'label': 'Harikane+', +} + +redshifts = [9, 12, 17] +wavelength = 1500. + +ULIM = -1e10 + +fits = {} + +fits['lf'] = {} + +fits['lf']['pars'] = \ +{ + 'Mstar': [], + 'pstar': [], + 'alpha': [], +} + +fits['lf']['err'] = \ +{ + 'Mstar': [], + 'pstar': [], + 'alpha': [], +} + +# Table 5: reminder, errors stored in (+/-) order, will be flipped in internal +# analysis routines to satisfy matplotlib's errorbar. +tmp_data = {} +tmp_data['lf'] = \ +{ + 9 : {'M': np.arange(-23.03, -17.03, 1), + 'phi': np.array([6.95e-5, 7.70e-5, 4.08e-5, 4.08e-5, 3.01e-4, 5.36e-4]), + 'err': [ULIM, ULIM, [9.59e-5, 3.92e-5], [9.59e-5, 3.92e-5], + [2.19e-4, 1.84e-4], [4.39e-4, 3.68e-4]], + }, + 12 : {'M': np.arange(-23.23, -17.23, 1), + 'phi': np.array([5.91e-6, 6.51e-6, 4.48e-6, 1.29e-5, 1.46e-5, 1.03e-4]), + 'err': [ULIM, ULIM, [10.34e-6, 3.82e-6], [1.72e-5, 0.87e-5], + [1.95e-5, 0.99e-5], [1.06e-4, 0.67e-4]], + }, + 17 : {'M': np.array([-23.59, -20.59]), + 'phi': np.array([2.15e-6, 5.39e-6]), + 'err': [ULIM, [7.16e-6, 3.57e-6]], + }, +} + +#for redshift in tmp_data['lf'].keys(): +# for i in range(len(tmp_data['lf'][redshift]['M'])): +# tmp_data['lf'][redshift]['phi'][i] *= 1e-3 +# +# if tmp_data['lf'][redshift]['err'][i] == ULIM: +# continue +# +# tmp_data['lf'][redshift]['err'][i][0] *= 1e-3 +# tmp_data['lf'][redshift]['err'][i][1] *= 1e-3 +# tmp_data['lf'][redshift]['err'][i] = \ +# tuple(tmp_data['lf'][redshift]['err'][i]) + +units = {'lf': 1.} + +data = {} +data['lf'] = {} +for key in tmp_data['lf']: + N = len(tmp_data['lf'][key]['M']) + mask = np.array([tmp_data['lf'][key]['err'][i] == ULIM for i in range(N)]) + + #mask = [] + #for element in tmp_data['lf'][key]['err']: + # if element == ULIM: + # mask.append(1) + # else: + # mask.append(0) + # + #mask = np.array(mask) + + data['lf'][key] = {} + data['lf'][key]['M'] = np.ma.array(tmp_data['lf'][key]['M'], mask=mask) + data['lf'][key]['phi'] = np.ma.array(tmp_data['lf'][key]['phi'], mask=mask) + data['lf'][key]['err'] = tmp_data['lf'][key]['err'] diff --git a/input/litdata/morishita2018.py b/input/litdata/morishita2018.py index 891973321..520ab9d18 100644 --- a/input/litdata/morishita2018.py +++ b/input/litdata/morishita2018.py @@ -1,7 +1,7 @@ """ Morishita et al., 2018, ApJ -Table 6. +Table 6. """ import numpy as np @@ -9,7 +9,7 @@ info = \ { 'reference': 'Morishita et al., 2018, MNRAS', - 'data': 'Table 3', + 'data': 'Table 3', 'label': 'Morishita+ (2018)', } @@ -29,7 +29,7 @@ 9.: {'M': [-23., -22., -21.], 'phi': [-5.9, -5.9, -5.4], 'err': [ULIM, (0.5, 0.8), (0.5, 0.8)], - }, + }, } units = {'lf': 'log10'} @@ -40,8 +40,8 @@ #mask = np.array(tmp_data['lf'][key]['err']) == ULIM N = len(tmp_data['lf'][key]['M']) mask = np.array([tmp_data['lf'][key]['err'][i] == ULIM for i in range(N)]) - + data['lf'][key] = {} - data['lf'][key]['M'] = np.ma.array(tmp_data['lf'][key]['M'], mask=mask) - data['lf'][key]['phi'] = np.ma.array(tmp_data['lf'][key]['phi'], mask=mask) + data['lf'][key]['M'] = np.ma.array(tmp_data['lf'][key]['M'], mask=mask) + data['lf'][key]['phi'] = np.ma.array(tmp_data['lf'][key]['phi'], mask=mask) data['lf'][key]['err'] = tmp_data['lf'][key]['err'] diff --git a/requirements.txt b/requirements.txt index c332b1cd1..a83943366 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy==1.22.0 +numpy>=1.22.2 matplotlib==2.2.4 scipy==1.2.1 h5py==2.9.0