<jupyter_start><jupyter_text>Assign halo mass to HSC galaxies using results of the A.S.A.P model<jupyter_code>test_dir = '../model/'
model_str = 'final_default'

# The configuration file
config_file = os.path.join(test_dir, 'asap_%s.yaml' % model_str)

# The results of the final sampling process
result_file = os.path.join(test_dir, 'asap_%s_sample.npz' % model_str)

# Initialize the model, load the data
cfg, params, obs_data, um_data = fitting.initial_model(config_file, verbose=True)

# Load in the final sampling results
(mod_result_samples, 
 mod_result_chains,
 mod_result_lnprob, 
 mod_result_best, _, _) = io.load_npz_results(result_file)

_, n_step, n_dim = mod_result_chains.shape

mod_result_best = np.nanmean(mod_result_chains[:, -int(n_step * 0.1):, :].reshape([-1, n_dim]), axis=0)

# Predict the stellar mass in inner and outer apertures
logms_inn, logms_tot, sig_logms, mask_use = predict_mstar_basic(
    um_data['um_mock'], mod_result_best, min_logms=10.5,
    logmh_col=cfg['um']['logmh_col'], min_scatter=cfg['um']['min_scatter'],
    pivot=cfg['um']['pivot_logmh'])

# Predict the SMFs and DeltaSigma profiles
um_smf_tot, um_smf_inn, um_dsigma = make_model_predictions(
    mod_result_best, cfg, obs_data, um_data)

# Check the likelihood for SMF and DeltaSigma profiles
lnlike_smf, lnlike_dsigma = ln_likelihood(
    mod_result_best, cfg, obs_data, um_data, sep_return=True)

# The useful part of the mock catalog
um_mock_use = um_data['um_mock'][mask_use]

print("\n# Best parameters: ", mod_result_best)
print("# ln(Likelihood) for SMFs   : %8.4f" % lnlike_smf)
print("# ln(Likelihood) for DSigma : %8.4f" % lnlike_dsigma)

print("\n", mod_result_best)
# Get the predicted aperture stellar mass and halo mass for central galaxies
um_mock_use = um_data['um_mock'][mask_use]

mask_cen = um_mock_use['upid'] == -1
um_mock_cen = um_mock_use[mask_cen]

logms_tot_cen = logms_tot[mask_cen]
logms_inn_cen = logms_inn[mask_cen]
logmh_vir_cen = um_mock_cen['logmh_vir']
logmh_peak_cen = um_mock_cen['logmh_peak']

# Aperture mass array for UM mock galaxies
um_mass = np.stack((logms_inn_cen, logms_tot_cen), axis=-1)<jupyter_output><empty_output><jupyter_text>Read in the HSC galaxy catalog<jupyter_code>hsc_mass = Table.read(
    os.path.join(cfg['obs']['dir'], 's16a_wide2_massive_fsps1_imgsub_use_short.fits'))

print("# There are %d HSC galaxies" % len(hsc_mass))
# Select relevant galaxies, exclude the ones with unreasonable mass and size
hsc_good = hsc_mass[(hsc_mass['r50_max'] >= 1.0) & (hsc_mass['r50_max'] <= 50.0) &
                    (hsc_mass['logm_max'] >= 11.2) & (hsc_mass['logm_max'] <= 12.5) &
                    (hsc_mass['logm_10'] >= 10.75) & (hsc_mass['logm_10'] <= 11.75)]

print("# There are %d good HSC massive galaxies" % len(hsc_good))

hsc_mhalo = copy.deepcopy(hsc_mass)

# Aperture mass array for HSC massive galaxies
hsc_aper_mass = np.stack((np.array(hsc_mass['logm_10']), np.array(hsc_mass['logm_max'])), axis=-1)<jupyter_output># There are 28860 good HSC massive galaxies<jupyter_text>Assign halo mass using random forest regressor<jupyter_code>from sklearn.ensemble import RandomForestRegressor

# Build the random forest regressor
mstar_forest = RandomForestRegressor(
    n_estimators=20, criterion='mse', max_depth=None, random_state=None, n_jobs=-1)

# Fit the data
mstar_forest.fit(um_mass, logmh_vir_cen)

# Predict halo mass
mvir_forest = mstar_forest.predict(hsc_aper_mass)

# Save the model
from sklearn.externals import joblib

joblib.dump(mstar_forest, '../data/um_asap_default_mvir_random_forest.pkl') 

# Add a column
hsc_mhalo.add_column(Column(data=mvir_forest, name='logmh_vir_forest'))
fig_1 = plt.figure(figsize=(6, 6))
fig_1.subplots_adjust(left=0.16, right=0.995, bottom=0.13, top=0.995, wspace=0.00, hspace=0.00)

ax1 = fig_1.add_subplot(111)

ax1.grid(linestyle='--', linewidth=2, alpha=0.3, zorder=0)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(20)
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(20)

# Scatter
um_mvir_forest = mstar_forest.predict(um_mass)

ax1.scatter(um_mock_cen['logmh_vir'], um_mvir_forest, s=2, alpha=0.2, c='darkgray')

# Mean and std
logmh_bin = np.linspace(12.5, 15.4, 12)

idx_logmh = np.digitize(um_mock_cen['logmh_vir'], logmh_bin)

logmh_mean_um = [np.nanmedian(um_mock_cen['logmh_vir'][idx_logmh == k]) 
                 for k in range(len(logmh_bin))]

logmh_mean_forest = [np.nanmedian(um_mvir_forest[idx_logmh == k]) 
                     for k in range(len(logmh_bin))]

sigmh_forest = [np.nanstd(um_mvir_forest[idx_logmh == k]) 
                for k in range(len(logmh_bin))]

print(sigmh_forest)

ax1.scatter(logmh_mean_um, logmh_mean_forest, s=180, alpha=0.8, edgecolor='k',
            facecolor='orangered')

ax1.errorbar(logmh_mean_um, logmh_mean_forest, yerr=np.asarray(sigmh_forest)/2.0, 
             fmt='o', ecolor='maroon', elinewidth=2.5, capsize=5, capthick=2.0,
             color='orangered', markersize=0, alpha=0.8,
             markeredgecolor='none', zorder=1)
    
ax1.plot(np.linspace(12.5, 15.9, 10), np.linspace(12.5, 15.9, 10), 
         linestyle='--', linewidth=2.5, c='k', alpha=0.5)

ax1.set_xlim(12.4, 15.4)
ax1.set_ylim(12.4, 15.4) 

_ = ax1.set_xlabel(r'$\log\ (M_{\rm Vir,\ True}/M_{\odot})$', fontsize=25)
_ = ax1.set_ylabel(r'$\log\ (M_{\rm Vir,\ Forest/M_{\odot}})$', fontsize=25)<jupyter_output>[0.2628392001884028, 0.10823030235140932, 0.12129966140660885, 0.13170737248216138, 0.1337910398815224, 0.13514004620298672, 0.1361953190274504, 0.13651239775185608, 0.11022524011191917, 0.10096141294811198, 0.09985206971063003, 0.0]<jupyter_text>Assign halo mass using best-fit Mhalo-Mmax-M10 relation<jupyter_code>from lts_planefit import lts_planefit

# Only use massive central galaxies
mask_fit = (logms_tot_cen >= 11.60)

x_mmax = logms_tot_cen[mask_fit]
y_minn = logms_inn_cen[mask_fit]
z_mhalo = um_mock_cen[mask_fit]['logmh_vir']

x_sig = sig_logms[mask_cen][mask_fit] * 0.0 + 0.03
y_sig = sig_logms[mask_cen][mask_fit] * 0.0 + 0.02
z_sig = sig_logms[mask_cen][mask_fit]

p = lts_planefit(x_mmax, y_minn, z_mhalo, x_sig, y_sig, z_sig, 
                 frac=0.95, pivotx=np.mean(x_mmax), pivoty=np.mean(y_minn))

c_plane, a_plane, b_plane = p.abc
a_pivot, b_pivot = np.mean(x_mmax), np.mean(y_minn)
fig_2 = plt.figure(figsize=(6, 6))
fig_2.subplots_adjust(left=0.16, right=0.995, bottom=0.13, top=0.995, wspace=0.00, hspace=0.00)

ax1 = fig_2.add_subplot(111)

ax1.grid(linestyle='--', linewidth=2, alpha=0.3, zorder=0)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(20)
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(20)
    
# Scatter plot
logmh_plane_mod = (a_plane * (logms_tot_cen - a_pivot) + 
                   b_plane * (logms_inn_cen - b_pivot) + c_plane)

ax1.scatter(um_mock_cen['logmh_vir'], logmh_plane_mod, s=2, alpha=0.2, color='darkgray')

# Mean and std
logmh_bin = np.linspace(12.5, 15.4, 12)

idx_logmh = np.digitize(logmh_plane_mod, logmh_bin)

logmh_mean_um = [np.nanmedian(um_mock_cen['logmh_vir'][idx_logmh == k]) 
                 for k in range(len(logmh_bin))]

logmh_mean_plane = [np.nanmedian(logmh_plane_mod[idx_logmh == k]) 
                    for k in range(len(logmh_bin))]

sigmh_plane = [np.nanstd(um_mock_cen['logmh_vir'][idx_logmh == k]) 
               for k in range(len(logmh_bin))]

ax1.scatter(logmh_mean_um, logmh_mean_plane, s=180, alpha=0.8, edgecolor='k',
            facecolor='orangered')

ax1.errorbar(logmh_mean_um, logmh_mean_plane, yerr=np.asarray(sigmh_plane)/2.0, 
             fmt='o', ecolor='maroon', elinewidth=2.5, capsize=5, capthick=2.0,
             color='orangered', markersize=0, alpha=0.8,
             markeredgecolor='none', zorder=1)

# One-to-one relation
ax1.plot(np.linspace(12.5, 15.9, 10), np.linspace(12.5, 15.9, 10), 
         linestyle='--', linewidth=2.5, c='k', alpha=0.5)

ax1.set_xlim(12.4, 15.4)
ax1.set_ylim(12.4, 15.4) 

_ = ax1.set_ylabel(r'$\log\ (M_{\rm Vir,\ Plane/M_{\odot}})$', fontsize=25)
_ = ax1.set_xlabel(r'$\log\ (M_{\rm Vir,\ True}/M_{\odot})$', fontsize=25)<jupyter_output><empty_output><jupyter_text>The best-fit aperture mass plane * 13.6826 + 3.541 * (logm_max - 11.71) - 2.542 * (logm_10 - 11.33)<jupyter_code># Best-fit plane
mvir_plane = a_plane * (hsc_mhalo['logm_max'] - a_pivot) + b_plane * (hsc_mhalo['logm_10'] - b_pivot) + c_plane

hsc_mhalo.add_column(Column(data=mvir_plane, name='logmh_vir_plane'))
def assign_mvir_plane(logm_max, logm_10, a_plane=3.541, a_pivot=11.71, b_plane=-2.542, b_pivot=11.33,
                      c_plane=13.6826):
    """Assign Mvir to galaxies with observed aperture masses."""
    return a_plane * (logm_max - a_pivot) + b_plane * (logm_10 - b_pivot) + c_plane<jupyter_output><empty_output><jupyter_text>Assign halo mass using more complex relation built by symbolic regression<jupyter_code>from gplearn.genetic import SymbolicRegressor

um_mvir = um_mock_cen['logmh_vir']

mask_fit = logms_tot_cen >= 11.55

est_gp = SymbolicRegressor(population_size=5000, 
                           const_range=(-50.0, 50.0),
                           metric='mse',
                           generations=2000, stopping_criteria=0.002,
                           p_crossover=0.9, p_subtree_mutation=0.01,
                           p_hoist_mutation=0.02, p_point_mutation=0.01,
                           max_samples=0.95, verbose=0, n_jobs=-1,
                           parsimony_coefficient=0.001, 
                           random_state=0)

est_gp.fit(um_mass[mask_fit], um_mvir[mask_fit])

print(est_gp._program)
with open('../data/um_asap_default_mvir_symbolic_regression.pkl', 'wb') as f:
    pickle.dump(est_gp, f)
    
f.close()
mvir_symbol = est_gp.predict(hsc_aper_mass)

hsc_mhalo['logmh_vir_symbol'] = mvir_symbol
#hsc_mhalo.add_column(Column(data=mvir_symbol, name='logmh_vir_symbol'))
fig_3 = plt.figure(figsize=(6, 6))
fig_3.subplots_adjust(left=0.16, right=0.995, bottom=0.13, top=0.995, wspace=0.00, hspace=0.00)

ax1 = fig_3.add_subplot(111)

ax1.grid(linestyle='--', linewidth=2, alpha=0.3, zorder=0)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(20)
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(20)
    
# Scatter plot

logmh_symbol_mod = est_gp.predict(um_mass) 

ax1.scatter(um_mock_cen['logmh_vir'], logmh_symbol_mod, s=2, alpha=0.1, c='skyblue')

# Mean and std
logmh_bin = np.linspace(12.5, 15.4, 12)

idx_logmh = np.digitize(logmh_symbol_mod, logmh_bin)

logmh_mean_um = [np.nanmedian(um_mock_cen['logmh_vir'][idx_logmh == k]) 
                 for k in range(len(logmh_bin))]

logmh_mean_symbol = [np.nanmedian(logmh_symbol_mod[idx_logmh == k]) 
                    for k in range(len(logmh_bin))]

sigmh_symbol = [np.nanstd(um_mock_cen['logmh_vir'][idx_logmh == k]) 
               for k in range(len(logmh_bin))]

ax1.scatter(logmh_mean_um, logmh_mean_symbol, s=180, alpha=0.8, edgecolor='k',
            facecolor='orangered')

ax1.errorbar(logmh_mean_um, logmh_mean_symbol, yerr=np.asarray(sigmh_symbol)/2.0, 
             fmt='o', ecolor='maroon', elinewidth=2.5, capsize=5, capthick=2.0,
             color='orangered', markersize=0, alpha=0.8,
             markeredgecolor='none', zorder=1)

# One-to-one relation
ax1.plot(np.linspace(12.5, 15.9, 10), np.linspace(12.5, 15.9, 10), 
         linestyle='--', linewidth=2.5, c='k', alpha=0.5)

ax1.set_xlim(12.4, 15.4)
ax1.set_ylim(12.4, 15.4) 

_ = ax1.set_xlabel(r'$\log\ (M_{\rm Vir,\ True}/M_{\odot})$', fontsize=25)
_ = ax1.set_ylabel(r'$\log\ (M_{\rm Vir,\ Symbol/M_{\odot}})$', fontsize=25)<jupyter_output><empty_output><jupyter_text>Predict virial radius<jupyter_code># Read the UM mock catalog with additional information
um_vagc = np.load(
    os.path.join(cfg['um']['dir'], 'um_smdpl_insitu_exsitu_0.7124_basic_logmp_11.5.npy'))

mask_mass = um_vagc[cfg['um']['halo_col']] >= cfg['um']['min_mvir']

um_vagc_use = um_vagc[mask_mass]

# Make sure the galaxies are lined up
assert np.all(um_vagc_use['halo_id'] == um_data['um_mock']['halo_id'])
from sklearn.ensemble import RandomForestRegressor

# Build the random forest regressor
rvir_forest = RandomForestRegressor(
    n_estimators=15, criterion='mse', max_depth=None, random_state=None, n_jobs=-1)

logr_vir_um_cen = np.log10(um_vagc_use[mask_use][mask_cen]['rvir_hlist'])

# Fit the data
rvir_forest.fit(um_mass, logr_vir_um_cen)

# Predict halo mass
logr_forest = rvir_forest.predict(hsc_aper_mass)

# Add a column
hsc_mhalo.add_column(Column(data=logr_forest, name='logr_vir_forest'))

# Save the model
#joblib.dump(mstar_forest, 'data/um_asap_default_rvir_random_forest.pkl') 
logmvir_hsc = (
    3.541 * (hsc_mhalo['logm_max'] - 11.71) - 
    2.542 * (hsc_mhalo['logm_10'] - 11.33) + 13.6826)
from colossus.halo import mass_so

rvir = mass_so.M_to_R(10.0 ** logmvir_hsc / 0.7, 0.4, 'vir')
plt.scatter(logr_forest, np.log10(rvir), alpha=0.3)

plt.plot([2.6, 3.2], [2.6, 3.2], c='r')
fig_4 = plt.figure(figsize=(6, 6))
fig_4.subplots_adjust(left=0.16, right=0.995, bottom=0.13, top=0.995, wspace=0.00, hspace=0.00)

ax1 = fig_4.add_subplot(111)

ax1.grid(linestyle='--', linewidth=2, alpha=0.3, zorder=0)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(20)
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(20)

# Scatter
um_rvir_forest = rvir_forest.predict(um_mass)

ax1.scatter(logr_vir_um_cen, um_rvir_forest, s=2, alpha=0.2, c='darkgray')

# Mean and std
logrvir_bin = np.linspace(2.1, 3.3, 8)

idx_logmh = np.digitize(logr_vir_um_cen, logrvir_bin)

logrvir_mean_um = [np.nanmedian(logr_vir_um_cen[idx_logmh == k]) 
                   for k in range(len(logrvir_bin))]

logrvir_mean_forest = [np.nanmedian(um_rvir_forest[idx_logmh == k]) 
                       for k in range(len(logrvir_bin))]

sigrvir_forest = [np.nanstd(um_rvir_forest[idx_logmh == k]) 
                  for k in range(len(logrvir_bin))]

ax1.scatter(logrvir_mean_um, logrvir_mean_forest, s=180, alpha=0.8, edgecolor='k',
            facecolor='orangered')

ax1.errorbar(logrvir_mean_um, logrvir_mean_forest, yerr=np.asarray(sigrvir_forest)/2.0, 
             fmt='o', ecolor='maroon', elinewidth=2.5, capsize=5, capthick=2.0,
             color='orangered', markersize=0, alpha=0.8,
             markeredgecolor='none', zorder=1)
    
ax1.plot(np.linspace(0.5, 3.9, 10), np.linspace(0.5, 3.9, 10), 
         linestyle='--', linewidth=2.5, c='k', alpha=0.5)

ax1.set_xlim(2.1, 3.3)
ax1.set_ylim(2.1, 3.3) 

_ = ax1.set_xlabel(r'$\log\ (R_{\rm Vir,\ True}/{\rm kpc})$', fontsize=25)
_ = ax1.set_ylabel(r'$\log\ (R_{\rm Vir,\ Forest/{\rm kpc}})$', fontsize=25)
fig_5 = plt.figure(figsize=(6, 6))
fig_5.subplots_adjust(left=0.16, right=0.995, bottom=0.13, top=0.995, wspace=0.00, hspace=0.00)

ax1 = fig_5.add_subplot(111)

ax1.grid(linestyle='--', linewidth=2, alpha=0.3, zorder=0)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(20)
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(20)
    
ms_mask = hsc_mhalo['logm_max'] >= 11.5
    
SC = ax1.scatter(logr_forest[ms_mask], np.log10(hsc_mhalo[ms_mask]['r50_100']), alpha=0.6, 
                 c=(hsc_mhalo[ms_mask]['logm_10']), s=30,
                 label='__no_label__', cmap='coolwarm', vmin=10.9, vmax=11.5)

# r_n = 0.015 R200
rvir_arr = np.linspace(100, 2000, 500)
r50_arr = 0.015 * rvir_arr

ax1.plot(np.log10(rvir_arr), np.log10(r50_arr), linestyle='--', linewidth=3.0,
         c='maroon', label=r'$r_{50} = 0.015 \times R_{\rm Vir}$')

#_ = ax1.text(3.00, 0.8, r'$\mathrm{Size:\ }\log M_{\mathrm{Vir,UM}}$', fontsize=16)

# Colorbar
cax = fig_5.add_axes([0.45, 0.18, 0.50, 0.03])
cbar = plt.colorbar(SC, cax=cax, orientation='horizontal')
cbar.solids.set_edgecolor("face")

_ = ax1.text(0.53, 0.11, r'$\log M_{\star,\ 10\ \rm kpc}$', 
             transform=ax1.transAxes, size=20)

for tick in cax.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in cax.yaxis.get_major_ticks():
    tick.label.set_fontsize(0)

ax1.legend(loc='upper left', fontsize=20)

ax1.set_xlim(2.61, 3.29)
ax1.set_ylim(0.62, 1.59)

_ = ax1.set_xlabel(r'$\log\ (R_{\rm Vir, \ UM}/\mathrm{kpc})$', fontsize=20)
_ = ax1.set_ylabel(r'$\log\ (R_{50}/\mathrm{kpc})$', fontsize=20)
from colossus.cosmology import cosmology
from colossus.halo import mass_so
from colossus.halo import mass_adv

params = {'flat': True, 'H0': 70.0, 'Om0': 0.30, 'Ob0': 0.049, 
          'sigma8': 0.81, 'ns': 0.95}
cosmo = cosmology.setCosmology('myCosmo', params)

m200c, r200c, c200c = mass_adv.changeMassDefinitionCModel(
    10.0 ** logmvir_hsc / 0.7, 0.4, 'vir', '200c')<jupyter_output><empty_output><jupyter_text>Save the HSC catalogs with Mvir<jupyter_code>hsc_mhalo.write('../data/s16a_wide2_massive_fsps1_imgsub_use_mvir_added.fits', format='fits',
                overwrite=True)
import emcee
import corner

def log_prob_D(params):
    m = params[:2]
    b, log_lambda = params[2:]
    v = np.append(-m, 1.0)
    
    # Enforce the log-uniform prior on lambda. 
    if not -5.0 < log_lambda < 5.0:
        return -np.inf
        
    # Compute \Sigma^2 and \Delta.
    Sigma2 = np.dot(np.dot(S, v), v) + np.exp(2*log_lambda)
    Delta = np.dot(X, v) - b
    
    # Compute the log likelihood up to a constant.
    ll = -0.5 * np.sum(Delta**2 / Sigma2 + np.log(Sigma2))
    return ll

# Run the MCMC.
sampler_D = emcee.EnsembleSampler(nwalkers, 4, log_prob_D)
p0 = np.append(m_true, [b_true, np.log(lambda_true)])
p0 = p0 + 1e-4 * np.random.randn(nwalkers, len(p0))
pos, _, _ = sampler_D.run_mcmc(p0, 500)
sampler_D.reset()
sampler_D.run_mcmc(pos, 5000)
samples_D = sampler_D.flatchain

tau = sampler_D.get_autocorr_time(c=4)
nsamples = len(samples_D)
print("{0:.0f} independent samples of m1".format(nsamples / tau[0]))
print("{0:.0f} independent samples of m2".format(nsamples / tau[1]))
print("{0:.0f} independent samples of b".format(nsamples / tau[2]))
print("{0:.0f} independent samples of ln(lambda)".format(nsamples / tau[3]))

corner.corner(samples_D, labels=["$m_1$", "$m_2$", "b", "$\ln\lambda$"],
              truths=np.append(m_true, [b_true, np.log(lambda_true)]));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15510697981409163, 'fraction_numerical': 0.06094274242340905, 'mean_word_length': 4.352369738004121, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 38, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '304360', 'license': 'Unlicense', 'n_tokens_mistral': 8411, 'n_tokens_neox': 7653, 'n_words': 1259, 'path': 'note/assign_halo_mass.ipynb', 'repo_name': 'dr-guangtou/asap', 'size': '18181'}