<jupyter_start><jupyter_text>Read in TNG SFHs<jupyter_code>f_sfh = os.path.join(os.environ['HOME'], 'data', 'gqp_mc', 'tng', 'IQ_TNG_SFRHs-corrected.hdf5')
fsfh = h5py.File(f_sfh, 'r')

sfhs = fsfh['SFRH'][...] # stellar mass formed in bins of time and metallicity
# we don't care much about metallicity bins so lets add them up.
# also the last min in SFH doesn't make sense
sfhs_t = np.sum(sfhs, axis=2)[:,:-1] # SFH(t)
# read in corresponding bin edges
_t_low, _t_high = np.loadtxt(os.path.join(os.environ['HOME'], 'data', 'gqp_mc', 'tng', 't_sfh.txt'), delimiter=',', unpack=True, usecols=[0,1])
t_low = _t_low[:-1]
t_high = _t_high[:-1]
t_mid = 0.5 * (t_low + t_high)
dt = t_high - t_low
dt<jupyter_output><empty_output><jupyter_text>Recent SFHs are binned way too finely. Lets rebin in cosmic time<jupyter_code>def rebin(edges0, sfh0, edges1): 
    '''given edges and SFH rebin to edges1
    '''
    sfh0 = np.atleast_2d(sfh0)
    
    nbin = len(edges1) - 1
    sfh1 = np.empty((sfh0.shape[0], nbin))
    for i in range(nbin): 
        xp0, xp1 = edges1[i], edges1[i+1]
        
        # add up old bins that entirely fall inside
        i0 = np.where(edges0 >= xp0)[0][0]
        i1 = np.where(edges0 < xp1)[0][-1]
        
        sfh1[:,i] = np.sum(sfh0[:,i0:i1], axis=1)
        # add partial contributiones 
        if edges0[i0] > xp0: 
            sfh1[:,i] += sfh0[:,i0-1] * (edges0[i0] - xp0)/(edges0[i0] - edges0[i0-1])
        
        if edges0[i1] < xp1:
            sfh1[:,i] += sfh0[:,i1] * (xp1 - edges0[i1])/(edges0[i1+1] - edges0[i1])
    return sfh1
t_low_coarse = 0.25 * np.arange(55)
t_high_coarse = t_low_coarse + 0.25
t_mid_coarse = 0.5 * (t_low_coarse + t_high_coarse)
dt_coarse = t_high_coarse - t_low_coarse
sfhs_t_coarse = rebin(np.concatenate([t_low, [t_high[-1]]]), sfhs_t, 0.25*np.arange(56))
print((np.sum(sfhs_t, axis=1) - np.sum(sfhs_t_coarse, axis=1)).max())
sfrhs_t = sfhs_t_coarse / dt_coarse # divide by bin width to get SFR
Msum = np.trapz(sfrhs_t, t_mid_coarse) # noramlization 

normed_sfrhs_t = sfrhs_t / Msum[:,None]<jupyter_output><empty_output><jupyter_text>Lets smooth out the SFHs to remove any binning artifacts<jupyter_code>def smooth_SFH(y, smooth_std):
    g = Box1DKernel(smooth_std, mode='center')
    y_smooth = convolve(y, g, boundary='extend')
    return y_smooth
smooth_normed_sfrhs_t = np.array([smooth_SFH(normed_sfrhs_t[i,:], 10) for i in range(normed_sfrhs_t.shape[0])])
fig = plt.figure(figsize=(10, 6))
sub = fig.add_subplot(111)
for i in np.random.choice(np.arange(sfhs.shape[0]), size=1000): 
    sub.plot(t_mid_coarse, smooth_normed_sfrhs_t[i])
sub.set_xlabel(r'$t_{\rm look~back}$', fontsize=25)
sub.set_xlim(0., 14)
sub.set_ylim(0., 1)
n_comp = 6
model = NMF(n_components=n_comp, init='random', random_state=0)
W = model.fit_transform(smooth_normed_sfrhs_t[Msum > 1e1])
H = model.components_

fig = plt.figure(figsize=(10,6))
sub = fig.add_subplot(111)
for _H in H: 
    sub.plot(t_mid_coarse, _H)
sub.set_xlabel(r'$t_{\rm look~back}$', fontsize=25)
sub.set_xlim(0., 13.7)
sub.set_ylim(0, None)<jupyter_output><empty_output><jupyter_text>Lets save the components<jupyter_code>np.savetxt('t_sfh.tng%icomp.txt' % n_comp, t_low_coarse)
np.savetxt('NMF_basis.sfh.tng%icomp.txt' % n_comp, H.T)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup