<jupyter_start><jupyter_text>Measure contribution of Major merger, minor merger, and smooth accretion for only the 'safe' samples...?. Because tree bad link more likely occur at major merger events, I guess the 'safe' samples have less major mergers than the total sample.   import analysis.Major_Minor_accretion as mma Largest halo in the refinement region is not the main cluster.Is the largest NP halo the main cluster? To check it, color halos in NP.No, max_np galaxy/halo is the main galaxy/halo.But 'rvir' value is wrong.and fixed.<jupyter_code>import utils.sampling as smp
import matplotlib.pyplot as plt
import tree
import pickle
import tree.halomodule as hmo
import numpy as np
from analysis.misc import load_cat
import scipy.stats
import tree.ctutils as ctu
from analysis.evol_lambda import MainPrg
import draw
import load
import analysis.evol_lambda as evl
import analysis.Major_Minor_accretion as mma
import analysis.misc as amsc
import tree.ctutils as ctu
import utils.match as mtc
# Read a single galaxy evolution catalog.

from analysis.MajorMinorAccretion_module import *
from analysis.all_plot_modules import *

verbose=True
# In[4]:
base = './'
cdir = ['catalog/', 'HM/', 'catalog_GM/', "easy_final/"][3]

clusters = ['01605', '07206', \
            '35663', '24954', '49096', \
            '05427', '05420', '29172', \
            '29176', '10002', '36415', \
            '06098', '39990', '36413', \
            '17891', '04466']
# parameters used for lambda_arr clipping.
#ind_upper = 20
#ind_lower = 20
#sig_upper = 2.0
#sig_lower = 2.0

 # 62: z = 1.666
nout_fi = 187

minimum_good_snap = 87

def body(clusters,
         dist_gal_scale_in=5,
         dist_gal_scale_out=10,
         dt_before=0.5,
         dt_after=0.5,
         dt_settle=0.5,
         load=False,
         nout_ini = 37,
         filter_small=True,
         min_mass_ratio = 0.05,
         measure_delta_savefig=False,
         find_merger_epoch_plot=False,
         cdir=""):

    suffix = "_{}_{}_{}_{}_{}_{}_{}".format(dist_gal_scale_in,
    dist_gal_scale_out,dt_before,dt_after,dt_settle,nout_ini,
                                           min_mass_ratio)
    if filter_small: 
        suffix = suffix + "_filtered_"

    if load:
        return pickle.load(open("main_prgs_final_augmented" + suffix + ".pickle", 'rb'))
        
    else:
        mpgs = []
        for cluster in clusters:
            print(cluster)
            wdir = base + cluster + '/'

            # Serialize catalogs. -> Only main galaxies
            # main galaxy list
            alltrees = ctu.load_tree(wdir, is_gal=True)
            ad = alltrees.data
            tn = ad[ad['nout'] == nout_fi]

            cat = load_cat(wdir + cdir + 'catalog' + str(nout_fi) + '.pickle')
            #idx_all = [tn['id'][tn['Orig_halo_id'] == id_final][0] for id_final in cat['id']]
            idx_all = cat['idx'][cat["idx"] > 0].astype(int) # why idx are float???

            mpg_tmp = []
            for i, idx in enumerate(idx_all):
                #print(i, idx)

                mpg_tmp.append(MainPrg(ad, idx))
            #    mpg_tmp =[MainPrg(ad, idx) for idx in idx_all]
            for nout in range(nout_ini, nout_fi + 1):
                cat = load_cat(wdir + cdir + 'catalog' + str(nout) + '.pickle')
                for gal in mpg_tmp:
                    gal.set_data(cat, nout)
                    gal.cluster = int(cluster)
        #        print(nout)
            # get rid of galaxies with too short tree.
            mpg_tmp = [gg for gg in mpg_tmp if sum(gg.data["reff"] > 0) > minimum_good_snap]
            for gal in mpg_tmp:
                gal.fill_missing_data()
                gal.clip_non_detection()
                gal.smoothed_lambda_org = mma.smooth(gal.data["lambda_r"], window_len=15)[:-1]
                gal.smoothed_r = mma.smooth(gal.data["reff"], window_len=15)[:-1]
                gal.smoothed_lambda = mma.smooth(l_at_smoothed_r(gal, npix_per_reff=5), window_len=15)[:-1]

            # save for each cluser
            with open(wdir + "main_prgs" + suffix + ".pickle", "wb") as f:
                pickle.dump(mpg_tmp, f)    
                
            # Find_merger_epochs needs smoothed_r
            find_merger_epochs(alltrees,
                               idx_all,
                               mpg_tmp,
                               nout_ini=nout_ini,
                               dist_gal_scale_in=dist_gal_scale_in,
                               dist_gal_scale_out=dist_gal_scale_out,
                               min_mass_ratio = min_mass_ratio,
                               mass_ratio='early',
                               verbose=False,
                               do_plot=find_merger_epoch_plot,
                               max_rgal=40,
                               pdf_fname=str(cluster) + "merger_ratio_epoch" + suffix + ".pdf")
            

            while len(mpg_tmp) > 0:
                mpgs.append(mpg_tmp.pop())

        if filter_small:
            for gal in mpgs:
                # Keep only the largest merger among multiple mergers
                filter_small_mergers(gal.merger)

        with open("main_prgs" + suffix + ".pickle", 'wb') as f:
            mpgs.pop(2)
            pickle.dump(mpgs, f)

        measure_delta(mpgs,
                      dt_before=dt_before,
                      dt_after=dt_after,
                      dt_settle=dt_settle,
                      nout_ini=nout_ini,
                      savefig=measure_delta_savefig,
                      figname="figs/measure_delta" + suffix)

        Maj_min_acc_ratio(mpgs, major_ratio=4)

        with open("main_prgs_final_augmented" + suffix + ".pickle", 'wb') as f:
            pickle.dump(mpgs, f)
    return mpgs
mpgs = body(clusters,
        dist_gal_scale_in=5,
        dist_gal_scale_out=10,
        dt_before=0.5,
        dt_after=0.5,
        dt_settle=0.5,
        load=False,
        nout_ini=37,
        filter_small=True,
        min_mass_ratio = 0.01, # 0.1로 테스트
        measure_delta_savefig=False,
        find_merger_epoch_plot=False,
        cdir=cdir)# z = 3
wdir = '/home/hoseung/Work/data/'

mpgs = pickle.load(open(wdir + "main_prgs_final_augmented_5_10_0.5_0.5_0.5_37_0.01_filtered_.pickle", "rb"))
#def kde_sci(mpgs
mstar_cut_hard = 5e9
mcut = 3.3e10
hist=False
shade=False
kde=True
norm_hist=False
detected=True
maj_ratio=4
excess=False
fname="figs/MajMinNon_contribution_per_gal_"
img_scale=1.5


fontsize_ticks = 6 * img_scale
fontsize_tick_label = 8 * img_scale
fontsize_legend = 5 * img_scale

from matplotlib.ticker import NullFormatter


l_dl_e, l_mr_e, l_mass_e = [], [], []
s_dl_e, s_mr_e, s_mass_e = [], [], []
l_dlt_g, l_dlo_g, l_dlM_g, l_dlm_g, l_mass_g = [], [], [], [], []
s_dlt_g, s_dlo_g, s_dlM_g, s_dlm_g, s_mass_g = [], [], [], [], []

M_changed = 0
m_changed = 0
no_merger_count = 0
count = 0
Maj_small = 0
for i, gal in enumerate(mpgs):
    mgal = gal.data["mstar"][0]
    if mgal > mstar_cut_hard:
        delta_lambda_tot = np.average(gal.data['lambda_r'][:5]) - np.average(gal.data['lambda_r'][-5:])
        delta_lambda_major = 0
        delta_lambda_minor = 0

        # Large
        if mgal > mcut:
            if hasattr(gal, "merger"):
                if gal.merger is not None:
                    l_dl_e.extend(gal.merger.delta_l)
                    l_mr_e.extend(gal.merger.mr)
                    for dl, mr in zip(gal.merger.delta_l, gal.merger.mr):
                        if (mr < maj_ratio) and (dl > -1):
                            delta_lambda_major = delta_lambda_major + dl
                        if (mr > maj_ratio) and (dl > -1):
                            delta_lambda_minor = delta_lambda_minor + dl

                delta_lambda_other = delta_lambda_tot - delta_lambda_major - delta_lambda_minor
                l_dlt_g.append(delta_lambda_tot)
                l_dlo_g.append(delta_lambda_other)
                l_dlM_g.append(delta_lambda_major)
                l_dlm_g.append(delta_lambda_minor)
        # small
        else:
            #s_mass_g.append(mgal)
            if hasattr(gal, "merger"):
                if gal.merger is not None:
                    s_dl_e.extend(gal.merger.delta_l)
                    s_mr_e.extend(gal.merger.mr)
                    for dl, mr in zip(gal.merger.delta_l, gal.merger.mr):
                        if (mr < maj_ratio) and (dl > -1):
                            delta_lambda_major = delta_lambda_major + dl
                        if (mr > maj_ratio) and (dl > -1):
                            delta_lambda_minor = delta_lambda_minor + dl


                delta_lambda_other = delta_lambda_tot - delta_lambda_major - delta_lambda_minor
                s_dlt_g.append(delta_lambda_tot)
                s_dlo_g.append(delta_lambda_other)
                s_dlM_g.append(delta_lambda_major)
                s_dlm_g.append(delta_lambda_minor)
# filter out a few galaxies with m < 5e9.
mpgs = [gal for gal in mpgs if gal.data["mstar"][0] > mstar_cut_hard]
dlt_all = np.array([gal.dlt for gal in mpgs])
dlM_all = np.array([gal.dlM for gal in mpgs])
dlm_all = np.array([gal.dlm for gal in mpgs])
dlo_all = np.array([gal.dlo for gal in mpgs])
mgals   = np.array([gal.data["mstar"][0] for gal in mpgs])
# All
fig, ax = plt.subplots()
ax.hist(dlt_all, histtype="step", label="total")
ax.hist(dlM_all, histtype="step", label="Major")
ax.hist(dlm_all, histtype="step", label="minor")
ax.hist(dlo_all, histtype="step", label="other")
plt.legend()
plt.show()
from matplotlib.ticker import NullFormatter

fig, axs = plt.subplots(2)
draw_kdes(dlM_all,
              dlm_all,
              dlo_all,
              dlt_all,
              axs[0],
              [1,2,3])

draw_kdes(dlM_all[dlM_all !=0],
              dlm_all[dlm_all !=0],
              dlo_all,
              dlt_all,
              axs[1],
              [1,2,3])

axs[0].yaxis.set_major_formatter(NullFormatter())
axs[0].set_ylabel("relative probability", fontsize=fontsize_tick_label)
axs[1].yaxis.set_major_formatter(NullFormatter())
axs[1].set_ylabel("relative probability", fontsize=fontsize_tick_label)
axs[0].set_xlabel(r"$\Delta\lambda")
axs[1].set_xlabel(r"$\Delta\lambda")
axs[0].legend()
plt.show()
def draw_kdes(dlM, dlm, dlo, dlt, ax
                  , nevents = None
                  , lw=1.5
                  , cov=0.2):
    
    dM = kde_den(dlM, cov=cov)
    dm = kde_den(dlm, cov=cov)
    do = kde_den(dlo, cov=cov)
    dtot = kde_den(dlt, cov=cov)

    xs=np.linspace(-0.7,0.7,51)
    dM_curve = dM(xs)
    dm_curve = dm(xs)
    do_curve = do(xs)
    dtot_curve = dtot(xs)

    nM = len(dlM)
    nm = len(dlm)
    no = len(dlo)
    ntot = len(dlt)


    Mlabel="Major"
    mlabel="Minor"
    olabel="non-merger"
    totlabel="Total"
    ax.plot(xs, dM_curve*nM/ntot, label=Mlabel, lw=lw, color="r")
    ax.plot(xs, dm_curve*nm/ntot, label=mlabel, lw=lw, color="g")        
    ax.plot(xs, do_curve*no/ntot, label=olabel, lw=lw, color="b")
    ax.plot(xs, dtot_curve, label=totlabel, lw=lw, color="black")

    ax.set_ylim([0, 1.15*ax.get_ylim()[1]])
np.array([dlM_all ==0]).sum()
fig, ax = plt.subplots()
#draw_kdes(dlM_all[dlM_all !=0], dlm_all[dlm_all !=0], dlo_all, dlt_all, ax2, lw=1, cov=0.2)
yp = 3.0
#ax.plot([0,0],[0,yp-0.4], color="black", lw=2)
ax.scatter([0,0],[3.5,4], marker="*", s=50)
#ax.plot([0, 0.2], [4,4], color="black", lw=1)
#ax.plot([0, 0.2], [5,5], color="black", lw=1)
ax.text(0.03, 3.9, "{} galaxies \nwithout Major merger".format(np.array([dlM_all ==0]).sum()))
ax.text(0.03, 3.4, "{} galaxies \nwithout minor merger".format(np.array([dlm_all ==0]).sum()))

#ax.annotate('Without Major mergers', xy=(0,yp -0.4 ), xytext=(0.1,0.0),
#            arrowprops={'arrowstyle': '-'}, va='center')
#ax.annotate('Without Major mergers', xy=(0,yp + 1), xytext=(0.1,0.0),
#            arrowprops={'arrowstyle': '-'}, va='center')
#ax.arrow(0, 0, 0, yp + 0.4, head_width=0.05, lw=1.5, head_length=0.1, fc='k', ec='k', arrowstyle='-')

#ax.arrow(0, yp + 0.4, 0, 1, head_width=0.05, lw=1.5, head_length=0.1, fc='k', ec='k')


draw_kdes(dlM_all[dlM_all !=0], dlm_all[dlm_all !=0], dlo_all, dlt_all, ax, lw=2, cov=0.2)

ax.set_ylim(0,4.5)
ax.set_yticks([0,1,2,3])
plt.show()
f, (ax, ax2) = plt.subplots(2, 1, sharex=True)

plot_type="hist"

# plot the same data on both axes
if plot_type=="kde":
    draw_kdes(dlM_all[dlM_all !=0], dlm_all[dlm_all !=0], dlo_all, dlt_all, ax, lw=1, cov=0.2)
    draw_kdes(dlM_all[dlM_all !=0], dlm_all[dlm_all !=0], dlo_all, dlt_all, ax2, lw=1, cov=0.2)
    
    
    
    
elif plot_type=="hist":
    nbins = 14
    hrange = [-0.75, 0.65]
    n, bins, patches = ax2.hist(dlt_all, histtype="step", label="total", bins=nbins, range=hrange)
    ax2.hist(dlM_all[dlM_all !=0], histtype="step", label="Major", bins=nbins, range=hrange)
    ax2.hist(dlm_all[dlm_all !=0], histtype="step", label="minor", bins=nbins, range=hrange)
    ax2.hist(dlo_all, histtype="step", label="other", bins=nbins, range=hrange)

    #ax2.hist(dlt_all, histtype="step", label="total")
    n2, bins2, patches2 = ax.hist(dlM_all[dlM_all ==0], histtype="step", label="Major"
                                  , bins=nbins, range=hrange)
    ax.hist(dlm_all[dlm_all ==0], histtype="step", label="minor", bins=nbins, range=hrange)
    #ax2.hist(dlo_all, histtype="step", label="other")
    
    ax.set_ylim(max(n) * 1.1, max(n2) * 1.1)  # outliers only
    ax2.set_ylim(0, max(n) * 1.1)  # most of the data
    
    

# zoom-in / limit the view to different portions of the data
#ax.set_ylim(3., 30)  # outliers only
#ax2.set_ylim(0, 3.)  # most of the data

# hide the spines between ax and ax2
ax.spines['bottom'].set_visible(False)
ax2.spines['top'].set_visible(False)
ax.xaxis.tick_top()
ax.tick_params(labeltop='off')  # don't put tick labels at the top
ax2.xaxis.tick_bottom()

plt.show()
plt.savefig(fname + "{:.1f}_11.png".format(np.log10(mcut)), dpi=200, bbox_inches="tight")
plt.savefig(fname + "{:.1f}_11.pdf".format(np.log10(mcut)), bbox_inches='tight') # eps does NOT support transparency!<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13591474245115454, 'fraction_numerical': 0.02763765541740675, 'mean_word_length': 2.379591836734694, 'pattern_counts': {'":': 2, '<': 8, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '350527', 'license': 'MIT', 'n_tokens_mistral': 5510, 'n_tokens_neox': 5184, 'n_words': 1023, 'path': 'scripts/scripts_ipynb/MajorMinorAccretion_contribution_per_gal_rev.ipynb', 'repo_name': 'Hoseung/pyRamAn', 'size': '14075'}