<jupyter_start><jupyter_text>Compare the outer envelope stellar mass and richness-based cluster finder DSigma profiles of HSC massive galaxies<jupyter_code># DeltaSigma profiles of HSC massive galaxies
topn_massive = pickle.load(open(os.path.join(res_dir, 'topn_galaxies_sum.pkl'), 'rb'))

# DeltaSigma profiles of redMaPPer and CAMIRA clusters
topn_cluster = pickle.load(open(os.path.join(res_dir, 'topn_clusters_cen_sum.pkl'), 'rb'))

# For clusters, but using both central and satellite galaxies
topn_cluster_all = pickle.load(open(os.path.join(res_dir, 'topn_clusters_sum.pkl'), 'rb'))<jupyter_output><empty_output><jupyter_text>DSigma profiles of mock galaxies<jupyter_code>sim_dsig = Table.read(os.path.join(sim_dir, 'sim_merge_all_dsig.fits'))<jupyter_output><empty_output><jupyter_text>Halo mass distributions<jupyter_code>sim_mhalo = Table.read(os.path.join(sim_dir, 'sim_merge_mhalo_hist.fits'))<jupyter_output><empty_output><jupyter_text>Main Figure<jupyter_code>dsig_cmp_2 = topn_cluster['cam_s16a_n_mem'][:3]
label_cmp_2 = r'$\rm CAMIRA$'

dsig_cmp = topn_cluster['redm_hsc_lambda']
label_cmp = r'$\rm redMaPPer$'
sub_cmp = r'{\rm redM}'

dsig_ref = topn_massive['logm_50_100']
label_ref = r'$M_{\star,\ [50, 100]}$'
sub_ref = r'{[50, 100]}'

fig = compare_dsigma_profiles(
    dsig_ref, dsig_cmp, sim_dsig, sim_mhalo, sig_type='jk', compare_to_model=True,
    label_ref=label_ref, label_cmp=label_cmp, sub_ref=sub_ref, sub_cmp=sub_cmp,
    marker_ref='o', marker_cmp='D', msize_ref=200, msize_cmp=160,
    dsig_cmp_2=dsig_cmp_2, label_cmp_2=label_cmp_2)
fig.savefig(os.path.join(fig_dir, 'fig_9.png'), dpi=120)
fig.savefig(os.path.join(fig_dir, 'fig_9.pdf'), dpi=120)<jupyter_output><empty_output><jupyter_text>Using the total stellar mass of CAMIRA<jupyter_code>dsig_cmp = topn_cluster['cam_s16a_logms'][:3]
label_cmp = r'$\rm CAMIRA$'
sub_cmp = r'{\rm CAM}'

#dsig_cmp = topn_cluster['redm_hsc_lambda']
#label_cmp = r'$\rm redMaPPer$'
#sub_cmp = r'{\rm redM}'

dsig_ref = topn_massive['logm_50_100'][:3]
label_ref = r'$M_{\star,\ [50, 100]}$'
sub_ref = r'{[50, 100]}'

fig = compare_dsigma_profiles(
    dsig_ref, dsig_cmp, sim_dsig, sim_mhalo, sig_type='jk', compare_to_model=True,
    label_ref=label_ref, label_cmp=label_cmp, sub_ref=sub_ref, sub_cmp=sub_cmp,
    marker_ref='o', marker_cmp='D', msize_ref=200, msize_cmp=160,
    dsig_cmp_2=None, label_cmp_2=None)<jupyter_output><ipython-input-41-5c5c3a4a4373>:22: RuntimeWarning: invalid value encountered in log10
  mod['r_mpc'], np.log10(ref_rand), fill_value='extrapolate')(r_mpc_obs)<jupyter_text>HSC redMaPPer: with spec-z v.s. all photo-z<jupyter_code>dsig_ref = topn_cluster['redm_hsc_lambda']
label_ref = r'${\rm redM:\ spec-}z$'
sub_ref = r'{\rm spec-}z'

dsig_cmp = topn_cluster['redm_hsc_photoz_lambda']
label_cmp = r'${\rm redM:\ phot-}z$'
sub_cmp = r'{\rm phot-}z'

fig = compare_dsigma_profiles(
    dsig_ref, dsig_cmp, sim_dsig, sim_mhalo, sig_type='jk', compare_to_model=False,
    label_ref=label_ref, label_cmp=label_cmp, sub_ref=sub_ref, sub_cmp=sub_cmp,
    marker_ref='o', marker_cmp='D', msize_ref=200, msize_cmp=160)
def compare_dsigma_profiles(dsig_ref, dsig_cmp, sim_dsig, sim_mhalo, sig_type='bt', compare_to_model=True,
                            label_ref=r'$\rm Ref$', label_cmp=r'$\rm Test$',
                            sub_ref=r'{\rm Ref}', sub_cmp=r'{\rm Test}',
                            cmap_list=None, color_bins=None, 
                            marker_ref='o', msize_ref=150,
                            marker_cmp='P', msize_cmp=180, 
                            show_best_cmp=False, middle_title=None, mvir_min=None,
                            dsig_cmp_2=None, label_cmp_2=r'$=rm Test 2$'):
    """Compare the Dsigma profiles."""

    def get_dsig_ratio(obs, ref, mod=None):
        """"""
        obs_rand = np.random.normal(
            loc=obs['dsigma'][0], scale=obs['dsig_err_{:s}'.format(sig_type)][0])

        if mod is not None:
            ref_rand = np.random.normal(
                loc=mod['dsig'], scale=(mod['dsig_err'] * err_factor))

            ref_inter = 10.0 ** (
                interpolate.interp1d(
                    mod['r_mpc'], np.log10(ref_rand), fill_value='extrapolate')(r_mpc_obs)
            )
            return obs_rand / ref_inter
        else:
            ref_rand = np.random.normal(
                loc=ref['dsigma'][0], scale=obs['dsig_err_{:s}'.format(sig_type)][0])
            return obs_rand / ref_rand


    # Color maps and bins
    if cmap_list is None:
        cmap_list = [
            palettable.colorbrewer.sequential.Blues_7_r,
            palettable.colorbrewer.sequential.OrRd_7_r,
            palettable.colorbrewer.sequential.YlGn_7_r,
            palettable.colorbrewer.sequential.Purples_7_r]

    if color_bins is None:
        color_bins = ["#377eb8", "#e41a1c", "#1b9e77", "#984ea3"]

    # Radius bin of the observed DSigma profiles
    r_mpc_obs = dsig_ref.meta['r_mpc']

    # ---- Start the figure ---- #
    fig = plt.figure(figsize=(16, 18))

    # Setup the figure
    n_col, n_bins = 3, len(dsig_ref)
    fig_y = int(4 * n_bins + 2)
    left, right = 0.08, 0.98
    if n_bins == 4:
        bottom, top = 0.055, 0.96
    elif n_bins == 3:
        bottom, top = 0.08, 0.96
    elif n_bins == 2:
        bottom, top = 0.11, 0.92
        
    x_space = 0.08
    x_size = (right - left - x_space * 1.05) / n_col
    y_size = (top - bottom) / n_bins
    
    fig = plt.figure(figsize=(16, fig_y))

    for bin_id in np.arange(len(dsig_ref)) + 1:
        # Setup the three columns
        ax1 = fig.add_axes([left, top - y_size * bin_id, x_size, y_size])
        ax2 = fig.add_axes([left + x_space + x_size, top - y_size * bin_id, x_size, y_size])
        ax3 = fig.add_axes([left + x_space + x_size * 2.04, top - y_size * bin_id, x_size, y_size])

        # Subplot title
        if bin_id == 1:
            ax1.set_title(r'$R \times \Delta\Sigma\ \rm Profile$', fontsize=38, pad=18)
            if middle_title is None:
                ax2.set_title(r'${\rm Richness\ v.s.}\ M_{\star,\ \rm Outer}$', fontsize=38, pad=18)
            else:
                ax2.set_title(middle_title, fontsize=38, pad=18)
            ax3.set_title(r'$M_{\rm vir}\ \rm Distribution$', fontsize=38, pad=18)

        # Color map
        cmap, color = cmap_list[bin_id - 1], color_bins[bin_id - 1]

        # MDPL halo mass information for this bin
        sim_dsig_bin = sim_dsig[sim_dsig['bin'] == bin_id - 1]
        sim_mhalo_bin = sim_mhalo[sim_mhalo['number_density_bin'] == bin_id - 1]

        # DSigma result for this bin
        dsig_ref_bin = dsig_ref[dsig_ref['bin_id'] == bin_id]
        dsig_cmp_bin = dsig_cmp[dsig_cmp['bin_id'] == bin_id]

        # Best fit DSigma profiles
        dsig_ref_best = sim_dsig_bin[
            np.argmin(
                np.abs(sim_dsig_bin['scatter'] - dsig_ref_bin['sig_med_{:s}'.format(sig_type)]))]
        dsig_cmp_best = sim_dsig_bin[
            np.argmin(
                np.abs(sim_dsig_bin['scatter'] - dsig_cmp_bin['sig_med_{:s}'.format(sig_type)]))]

        if dsig_ref_bin['sig_med_{:s}'.format(sig_type)] < 0.6:
            err_factor = 5.
        else:
            err_factor = 4.

        # Interpolated the reference model profile
        ref_model_inter = 10.0 ** (
            interpolate.interp1d(
                dsig_ref_best['r_mpc'], np.log10(dsig_ref_best['dsig']),
                fill_value='extrapolate')(r_mpc_obs)
        )

        if compare_to_model:
            ratio_sample = [
                get_dsig_ratio(
                    dsig_cmp_bin, dsig_ref_bin, mod=dsig_ref_best) for i in np.arange(2000)]
            ratio_cmp = dsig_cmp_bin['dsigma'][0] / ref_model_inter
        else:
            ratio_sample = [
                get_dsig_ratio(dsig_cmp_bin, dsig_ref_bin, mod=None) for i in np.arange(2000)]
            ratio_cmp = dsig_cmp_bin['dsigma'][0] / dsig_ref_bin['dsigma'][0]

        ratio_cmp_err_low = ratio_cmp - np.nanpercentile(ratio_sample, 16, axis=0)
        ratio_cmp_err_upp = np.nanpercentile(ratio_sample, 84, axis=0) - ratio_cmp
        
        if dsig_cmp_2 is not None:
            try:
                dsig_cmp_2_bin = dsig_cmp_2[dsig_cmp_2['bin_id'] == bin_id]
                dsig_cmp_2_best = sim_dsig_bin[
                    np.argmin(
                        np.abs(sim_dsig_bin['scatter'] - dsig_cmp_2_bin['sig_med_{:s}'.format(sig_type)]))]

                if compare_to_model:
                    ratio_sample = [
                        get_dsig_ratio(
                            dsig_cmp_2_bin, dsig_ref_bin, mod=dsig_ref_best) for i in np.arange(2000)]
                    ratio_cmp_2 = dsig_cmp_2_bin['dsigma'][0] / ref_model_inter
                else:
                    ratio_sample = [
                        get_dsig_ratio(dsig_cmp_2_bin, dsig_ref_bin, mod=None) for i in np.arange(2000)]
                    ratio_cmp_2 = dsig_cmp_2_bin['dsigma'][0] / dsig_ref_bin['dsigma'][0]

                ratio_cmp_2_err_low = ratio_cmp_2 - np.nanpercentile(ratio_sample, 16, axis=0)
                ratio_cmp_2_err_upp = np.nanpercentile(ratio_sample, 84, axis=0) - ratio_cmp_2
                show_cmp_2 = True
            except Exception:
                show_cmp_2 = False
        else:
            show_cmp_2 = False

        # Halo mass distribution from model
        if bin_id == 1:
            n_bins = 8
        elif bin_id == 2:
            n_bins = 4
        else:
            n_bins = 3

        mvir_true, hist_true, mvir_avg_true = catalog.rebin_mhalo_hist(
            sim_mhalo_bin, bin_id - 1, 0.0, n_bin=n_bins)

        mvir_ref, hist_ref, mvir_avg_ref = catalog.rebin_mhalo_hist(
            sim_mhalo_bin, bin_id - 1, dsig_ref_bin['sig_med_{:s}'.format(sig_type)], n_bin=20)
        mvir_cmp, hist_cmp, mvir_avg_cmp = catalog.rebin_mhalo_hist(
            sim_mhalo_bin, bin_id - 1, dsig_cmp_bin['sig_med_{:s}'.format(sig_type)], n_bin=20)

        # ----- Plot 1: R x DSigma plot ----- #
        ax1.set_xscale("log", nonpositive='clip')

        # MDPL: Best-fit
        ax1.fill_between(
            dsig_ref_best['r_mpc'],
            dsig_ref_best['r_mpc'] * (
                dsig_ref_best['dsig'] - dsig_ref_best['dsig_err'] * err_factor),
            dsig_ref_best['r_mpc'] * (
                dsig_ref_best['dsig'] + dsig_ref_best['dsig_err'] * err_factor),
            alpha=0.2, edgecolor='grey', linewidth=2.0,
            label=r'__no_label__', facecolor='grey', linestyle='-', rasterized=True)

        if show_best_cmp:
            ax1.fill_between(
                dsig_cmp_best['r_mpc'],
                dsig_cmp_best['r_mpc'] * (
                    dsig_cmp_best['dsig'] - dsig_cmp_best['dsig_err'] * err_factor),
                dsig_cmp_best['r_mpc'] * (
                    dsig_cmp_best['dsig'] + dsig_cmp_best['dsig_err'] * err_factor),
                alpha=0.15, edgecolor='grey', linewidth=2.0,
                label=r'__no_label__', facecolor='grey', linestyle='--', rasterized=True)

        # Reference DSigma profile
        ax1.errorbar(
            r_mpc_obs,
            r_mpc_obs * dsig_ref_bin['dsigma'][0],
            yerr=(r_mpc_obs * dsig_ref_bin['dsig_err_{:s}'.format(sig_type)][0]),
            ecolor=color, color=color, alpha=0.9, capsize=4, capthick=2.5, elinewidth=2.5,
            label='__no_label__', fmt='o', zorder=0)
        ax1.scatter(
            r_mpc_obs,
            r_mpc_obs * dsig_ref_bin['dsigma'][0],
            s=msize_ref, alpha=0.9, facecolor=color, edgecolor='w', marker=marker_ref, 
            linewidth=2.5, label=label_ref)

        # DSigma profiles to compare with
        ax1.errorbar(
            r_mpc_obs * 1.01,
            r_mpc_obs * dsig_cmp_bin['dsigma'][0],
            yerr=(r_mpc_obs * dsig_cmp_bin['dsig_err_{:s}'.format(sig_type)][0]),
            ecolor=cmap.mpl_colormap(0.5), color='w', alpha=0.9, capsize=4, capthick=2.5,
            elinewidth=2.5, label='__no_label__', fmt='o', zorder=0)
        ax1.scatter(
            r_mpc_obs * 1.01,
            r_mpc_obs * dsig_cmp_bin['dsigma'][0],
            s=msize_cmp, alpha=0.9, facecolor='w', edgecolor=cmap.mpl_colormap(0.5), 
            marker=marker_cmp, linewidth=3.0, label=label_cmp)

        y_max = np.max(
            [np.max(dsig_ref_best['r_mpc'] * dsig_ref_best['dsig']), 
             np.max(dsig_cmp_best['r_mpc'] * dsig_cmp_best['dsig'])]) * 1.5
        ax1.set_ylim(0.1, y_max)

        # Bin ID
        _ = ax1.text(
            0.08, 0.83, r'$\rm Bin\ {:1d}$'.format(bin_id), fontsize=35, transform=ax1.transAxes)
        
        if bin_id == 1:
            ax1.legend(loc='best', fontsize=20)

        if bin_id == len(dsig_ref):
            _ = ax1.set_xlabel(r'$R\ [\mathrm{Mpc}]$', fontsize=30)
        else:
            ax1.set_xticklabels([])
        _ = ax1.set_ylabel(r'$R \times \Delta\Sigma\ [10^{6}\ M_{\odot}/\mathrm{pc}]$', fontsize=30)


        # ----- Plot 2: Ratio of DSigma plot ----- #
        ax2.set_xscale("log", nonpositive='clip')

        ax2.axhline(
            1.0, linewidth=3.0, alpha=0.5, color='k', linestyle='--', label='__no_label__', )

        # Uncertainty of the model
        ax2.fill_between(
            dsig_ref_best['r_mpc'],
            1.0 - (dsig_ref_best['dsig_err'] * err_factor / dsig_ref_best['dsig']),
            1.0 + (dsig_ref_best['dsig_err'] * err_factor / dsig_ref_best['dsig']),
            alpha=0.2, edgecolor='none', linewidth=1.0, label='__no_label__',
            facecolor='grey', rasterized=True)

        ax2.errorbar(
            r_mpc_obs, ratio_cmp, yerr=[ratio_cmp_err_low, ratio_cmp_err_upp],
            ecolor=cmap.mpl_colormap(0.5), color='w', alpha=0.8, capsize=4, capthick=2.5,
            elinewidth=3.0, label='__no_label__', fmt='o', zorder=0)
        ax2.scatter(
            r_mpc_obs, ratio_cmp,
            s=msize_cmp, alpha=0.9, facecolor='w', edgecolor=cmap.mpl_colormap(0.5),
            marker=marker_cmp, linewidth=3.0, label=label_cmp)
        
        if show_cmp_2:
            ax2.errorbar(
                r_mpc_obs * 1.09, ratio_cmp_2, yerr=[ratio_cmp_2_err_low, ratio_cmp_2_err_upp],
                ecolor=cmap.mpl_colormap(0.5), color='w', alpha=0.6, capsize=4, capthick=2.5,
                elinewidth=3.0, label='__no_label__', fmt='o', zorder=0)
            ax2.scatter(
                r_mpc_obs * 1.09, ratio_cmp_2,
                s=165, alpha=0.7, facecolor=cmap.mpl_colormap(0.3), edgecolor='w',
                marker='s', linewidth=2.5, label=label_cmp_2)

        ax2.set_ylim(0.1, 1.89)

        if np.max(ratio_cmp) < 1.2:
            y_pos = 0.85
        else:
            y_pos = 0.15

        if bin_id == 1:
            ax2.legend(loc='best', fontsize=22)

        if bin_id == len(dsig_ref):
            _ = ax2.set_xlabel(r'$R\ [\mathrm{Mpc}]$', fontsize=30)
        else:
            ax2.set_xticklabels([])
            
        _ = ax2.set_ylabel(r'$\Delta\Sigma_{\rm Richness}/\Delta\Sigma_{[50, 100]}$', fontsize=30)


        # ----- Plot 3: Halo mass distribution plot ----- #

        # Histogram for sigma = 0.0
        ax3.fill_between(mvir_true, hist_true / hist_true.sum() / 1.7, color='grey',
                         step="pre", alpha=0.3, label=r'$\sigma_{\mathcal{M}|\mathcal{O}}=0.0$', zorder=0)
        ax3.axvline(mvir_avg_true, color='k', alpha=0.7, linewidth=4.0, linestyle='--')

        # Halo mass distribution for the reference sample
        ax3.fill_between(mvir_ref, hist_ref / hist_ref.sum(), color=color,
                         step="pre", alpha=0.5, label=label_ref, zorder=1)
        ax3.axvline(mvir_avg_ref, color=color, alpha=0.8, linewidth=4.0, linestyle='-.')

        # Halo mass distribution for the comparison sample
        ax3.fill_between(mvir_cmp, hist_cmp / hist_cmp.sum(), edgecolor=color, facecolor='none',
                         step="pre", alpha=0.8, linewidth=5, label=label_cmp, zorder=2)
        ax3.axvline(mvir_avg_cmp, color=color, alpha=0.8, linewidth=4.0, linestyle=':')

        if mvir_min is None:
            mvir_min = 11.85
        ax3.set_xlim(mvir_min, 15.35)

        if bin_id == 1:
            ax3.legend(loc='best', fontsize=25)

        ax3.axhline(0.0, linewidth=3.0, c='grey', alpha=0.7)
        ax3.set_yticklabels([])
        
        if bin_id == len(dsig_ref):
            _ = ax3.set_xlabel(r'$\log(M_{\rm vir}/M_{\odot})\ [\rm dex]$', fontsize=30)
        else:
            ax3.set_xticklabels([])
            
        for tick in ax1.xaxis.get_major_ticks():
            tick.label.set_fontsize(30)
        for tick in ax1.yaxis.get_major_ticks():
            tick.label.set_fontsize(30)

        for tick in ax2.xaxis.get_major_ticks():
            tick.label.set_fontsize(30)
        for tick in ax2.yaxis.get_major_ticks():
            tick.label.set_fontsize(30)

        for tick in ax3.xaxis.get_major_ticks():
            tick.label.set_fontsize(30)
        for tick in ax3.yaxis.get_major_ticks():
            tick.label.set_fontsize(30)

    return fig<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13733075435203096, 'fraction_numerical': 0.03170974737705879, 'mean_word_length': 2.456645056726094, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 25, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '78865', 'license': 'MIT', 'n_tokens_mistral': 6787, 'n_tokens_neox': 6417, 'n_words': 1019, 'path': 'notebooks/figure/fig9.ipynb', 'repo_name': 'mattkwiecien/jianbing', 'size': '17061'}