<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>
{'fraction_non_alphanumeric': 0.16246583662314, 'fraction_numerical': 0.0391740054661403, 'mean_word_length': 4.4088669950738915, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '636922', 'license': 'MIT', 'n_tokens_mistral': 1531, 'n_tokens_neox': 1467, 'n_words': 271, 'path': 'nb/tng/sfhs.ipynb', 'repo_name': 'changhoonhahn/GQP_mock_challenge', 'size': '3293'}