<jupyter_start><jupyter_text>delta M Vs delta Lambda<jupyter_code>fig, ax = plt.subplots()
dts = [1,3,5,10]
alldm =[[]]*len(dts)
alldl =[[]]*len(dts)
allmavg=[[]]*len(dts)
alllavg=[[]]*len(dts)
for i, dt in enumerate(dts):
    dm =[]
    dl =[]
    mavg =[]
    lavg =[]
    
    for gal in mpgs:
        ind = gal.data['lambda_r'] > 0
        mstar = gal.data['mstar'][ind]
        sbam = gal.data['lambda_r'][ind]
        #dm.extend((mstar[:-dt] - mstar[dt:])/mstar[dt:])
        #dl.extend(sbam[:-dt] - sbam[dt:])
        mavg = [sum(mstar[i*dt:(i+1)*dt])/dt for i in range(int(len(mstar)/dt))]
        lavg = [sum(sbam[i*dt:(i+1)*dt])/dt for i in range(int(len(sbam)/dt))]
        mavg = np.array(mavg)
        lavg = np.array(lavg)
        dm.extend((mavg[:-1] - mavg[1:])/mavg[1:])
        dl.extend(lavg[:-1] - lavg[1:])
    alldm[i] = dm
    alldl[i] = dl
    allmavg[i] = mavg
    alllavg[i] = lavg

    
    #dm = np.array((mavg[:1]-mavg[1:])/mavg[1:])
    #dl = np.array(lavg[:1]-lavg[1:])
do_contour=True
do_3d = False

fig, ax = plt.subplots(2,2)
ax = ax.ravel()


for i, dt in enumerate(dts):
    dm = np.array(alldm[i])
    dl = np.array(alldl[i])
    
    if do_contour:
        import scipy.stats as st
        xmin, xmax = -0.1, 0.1
        ymin, ymax = -0.1, 0.1

        xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
        positions = np.vstack([xx.ravel(), yy.ravel()])
        values = np.vstack([dm, dl])
        kernel = st.gaussian_kde(values)
        f = np.reshape(kernel(positions).T, xx.shape)
        f /= max(f.ravel())
        levels=[0.002, 0.02, 0.05, 0.2, 0.4, 0.95]
        if do_3d:
            from mpl_toolkits.mplot3d import Axes3D
            ax[i] = fig.gca(projection='3d')
            surf = ax[i].plot_surface(xx, yy, f, rstride=1, cstride=1, cmap='jet',
                           linewidth=0, antialiased=False)
        else:
            cfset = ax[i].contour(xx, yy, f, levels=levels, cmap='winter', lw=5)
    else:
        ax[i].scatter(dm,dl, edgecolors='none', alpha=0.1)
    ax[i].set_xlim([xmin,xmax])
    ax[i].set_ylim([ymin,ymax])
    ax[i].scatter([0],[0], s=[100], linewidth='3', marker='*', c='red')
    ax[i].set_title("dt = {}".format(dt))
    #ax[i].hist2d(dm,dl,bins = 20, range=([-0.2,0.4],[-0.2,0.2]), cmap='PuBu')
plt.tight_layout()
plt.show()
dl.mean()<jupyter_output><empty_output><jupyter_text>Mstar, Lambda all<jupyter_code>fig, ax = plt.subplots(2)
# * (gal.data['mstar'] > 1e11)
for gal in mpgs:
    ind = np.where(gal.data['lambda_r'] > 0)[0]
    mstar_ini = np.average(gal.data['mstar'][ind[-5:]])
    l_ini = np.average(gal.data['lambda_r'][ind[-5:]])
    ax[0].plot(gal.nouts[ind], np.log10(gal.data['mstar'][ind]/mstar_ini), c="grey", alpha=0.05)
    ax[0].set_ylim([-1,1])
    ax[1].plot(gal.nouts[ind], gal.data['lambda_r'][ind] - l_ini, c="grey", alpha=0.05)
    ax[1].set_ylim([-0.5,0.5])
plt.show()
for gal in mpgs:
    print(gal.data['idx'][0])<jupyter_output>241740
241535
241448
241363
241338
241306
241258
241228
241202
241176
241140
241565
241482
241403
241347
241312
241269
241240
241216
241183
241157
241581
241498
241411
241351
241318
241294
241243
241219
241191
241162
241685
241534
241446
241361
241337
241303
241255
241227
241200
241175
241646
241506
241418
241352
241327
241296
241248
241221
241193
241165
241564
241480
241400
241346
241311
241266
241238
241214
241181
241154
241649
241520
241426
241358
241328
241300
241252
241224
241196
241169
241665
241530
241437
241359
241330
241302
241253
241225
241197
241173
241566
241489
241410
241350
241314
241275
241242
241218
241189
241159
241545
241465
241369
241342
241309
241263
241235
241210
241179
241151
241563
241467
241376
241344
241310
241265
241237
241213
241180
241153
241540
241453
241368
241339
241308
241260
241233
241206
241178
241146
183678
183555
183482
183432
183403
183359
183580
183503
183442
183417
183386
183335
183607
183526
183450
183424
183393
183341
183588
183525
183448
183423[...]<jupyter_text>Merger epoch with Mstar vs lambda<jupyter_code>import tree.ctutils as ctu
from tree import treeutils
import numpy as np
import pickle

# Calculate merger event parameters
def find_merger(atree, idx=None, aexp_min=0.0):
    """
        find indices of merger event from a tree.
        (Full tree or main progenitor trunk)
    """
    if idx == None:
        idx = atree['id'][0]

    nprg = 1
    merger_list=[]

    i = 0
    while nprg > 0:
        idx = ctu.get_progenitors(atree, idx, main=True)[0]
        ind = np.where(atree['id'] == idx)[0]
        if atree['aexp'][ind] < aexp_min:
            break
        nprg = ctu.get_npr(atree, idx)

        if nprg > 1:
            merger_list.append(i)
        i +=1
    return merger_list



def merger_mass_ratio(atree, idx=None):
    """
    return mass ratio of the given merger event
    """
    if idx == None:
        idx = atree['id'][0]

    prgs = ctu.get_progenitors(atree, idx)

    # only for mergers
    if len(prgs) > 1:
        i_prgs = [np.where(atree['id'] == i)[0] for i in prgs]
        mass = []
        for iprg in i_prgs:
            mass.append(atree['m'])
    else:
        print("This is not a merger")
        return 0


def merger_properties_main_prg(atree, idx):
    """
        Calculate merger mass ratio for "one" merger event.

    if idx == None:
        if nout == None:
            print("Both idx and nout are missing")
            return
    else:
        if nout == None:
            nout = np.where(atree['id'] == idx)[0]

    idx = atree['id'][ind]
    """

    #prgs = get_progenitors(atree, idx)
    #if len(prgs) > 1:
    #    i_prgs = [np.where(atree['id'] == i)[0] for i in prgs]

    i_prgs = np.where(atree['desc_id'] == idx)[0]

    print(i_prgs)
    id_prgs = atree['id'][i_prgs]
    mass_prgs = atree['m'][i_prgs]

    #mass_prgs_norm = mass_prgs / sum(mass_prgs)

    return mass_prgs

    
def load_tree(wdir, is_gal=False, no_dump=False):
    import pickle
    from tree import treemodule
    import tree.ctutils as ctu

    alltrees = treemodule.CTree()


    if is_gal:
        # Galaxy tree
        tree_path = 'GalaxyMaker/Trees/'
    else:
        # halo tree
        tree_path = 'halo/Trees/'

    try:
        alltrees = pickle.load(open(wdir + tree_path + "extended_tree.pickle", "rb" ))
        print("Loaded an extended tree")
    except:
        alltrees = treemodule.CTree()
        alltrees.load(filename= wdir + tree_path + 'tree_0_0_0.dat')
        if not no_dump:
            # Fix nout -----------------------------------------------------
            nout_max = alltrees.data['nout'].max()
            alltrees.data['nout'] += 187 - nout_max
            print("------ NOUT fixed")
            alltrees.data = ctu.augment_tree(alltrees.data, wdir, is_gal=is_gal)
            print("------ tree data extended")

    return alltrees<jupyter_output><empty_output><jupyter_text>merger epochs<jupyter_code># multi page PDF
from matplotlib.backends.backend_pdf import PdfPages

fig, ax = plt.subplots(2, sharex=True)
plt.subplots_adjust(hspace=0.001)
with PdfPages('multipage_pdf.pdf') as pdf:
    for i in inds[0:3]:
        gal = mpgs[i]
        ax[0].scatter(gal.nouts, np.log10(gal.data['mstar']))
        ax[0].set_xlim([50,190])
        #ax[0].set_ylim([8,13])
        ax[1].plot(gal.nouts, gal.data['lambda_r'], 'r-')
        #ax[1].set_xlim([50,190])
        #ax[1].set_ylim([0,1])
        #print(gal.merger.nout, gal.merger.mr)
        for mr, xx in zip(gal.merger.mr, gal.merger.nout):
            ax[0].axvline(xx, linestyle=':')
            ax[0].annotate("{:.1f}".format(mr), xy=(xx,0.8))
            ax[1].axvline(xx, linestyle=':')
            ax[1].annotate("{:.1f}".format(mr), xy=(xx,0.8))
        pdf.savefig()
        ax[0].clear()
        ax[1].clear()
    plt.close()

#plt.show()<jupyter_output><empty_output><jupyter_text>1. 하나의 은하를 골라서 모든 머저를 표시하고 minor merger와 major merger중 어느게 더 전체 delta lambda에 기여를 많이 하는지 확인. 2. 모든 은하에 대해서, major merger가 많은 시기와 delta lambda가 큰 시기가 일치하는지 확인3. 은하단 계산 추가로 시작. (BCG를 10kpc, 15kpc 등으로 고정한 것도 추가)4. 나중에 트리가 바뀔 수도 있음.<jupyter_code>inds[0]
#%%
import matplotlib.pyplot as plt
# plot each galaxy.
# stellar mass growth and lambda_r as a function of time.

# The exponent (also called ass "offset") in the figure (1e11)
# overlaps with lookback time tick labels.
# And moving the offset around is not easy. 
# So, manually divide the values.

# compile mstar and lambda
mstar = []
lstar = []
for gal in mpgs:
    ind = np.where(gal.data['lambda_r'] > 0)[0]
    mstar = np.average(gal.data['mstar'][ind[-5:]])
    lstar = np.average(gal.data['lambda_r'][ind[-5:]])

mm = mstar/1e10

plt.close()
plt.ioff()

def make_patch_spines_invisible(ax):
    """
        Useful for plotting multiple variables (more than two twinx())
    """
    ax.set_frame_on(True)
    ax.patch.set_visible(False)
    #for sp in ax.spines.itervalues():  
    # Changed in Python3
    for sp in ax.spines.values():
        sp.set_visible(False)

for i, idgal in enumerate(cat['final_gal']):
#for i, idgal in enumerate([1618]):
    if mm[i][0] < 0.2 :
        print(idgal, mm[i][0], mm[i][-1])
        continue
    print(idgal, "!!!!")

    plt.rcParams["figure.figsize"] = [12,10]
    fig, axes = plt.subplots(3)
#    plt.figure(num=1, figsize=[10,20])
    fig.suptitle("ID: " + str(idgal).zfill(5), fontsize=18)#, y=1.01)
    lns1 = axes[0].plot(nouts[::-1], l_r[i], label=r"$\lambda_{R}$")
    axes[0].set_xticks(z_pos)
    axes[0].set_xticklabels(z_target_str)
    plt.subplots_adjust(left = 0.1, right = 0.9, \
                        wspace = 0.1, hspace = 0.0, \
                        bottom = 0.1, top = 0.85)
    
    axes[0].set_xlim([37,187])
    axes[0].set_ylim([0,1.0])
    axes[0].set_ylabel(r"$\lambda_{R}$")
    axes[0].set_xlabel("redshift")
    
#    ax2 = axes[0].twinx()
    lns2 = axes[1].plot(nouts[::-1], mm[i], 'r-', label="stellar mass")
    axes[1].set_ylim([0, 1.3*max(mm[i])])
    axes[1].set_xlim([37,187])
    axes[1].set_ylabel(r"Stellar mass $[10^{10}M_{\odot}]$")
    axes[1].get_yaxis().get_offset_text().set_y(1)
    
#    ax3 = ax1.twinx() # Reff
#    ax3.spines["right"].set_position(("axes", 1.2))
#    make_patch_spines_invisible(ax3)
    # Second, show the right spine.
#    ax3.spines["right"].set_visible(True)
    axes[2].set_ylabel("Reff [kpc]")
    axes[2].set_xlim([37,187])
    lns3 = axes[2].plot(nouts[::-1], reff[i], 'g-', label='Reff')

    # hide x axes so that subplots stick together.
    plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)
    
    ax4 = axes[0].twiny()
    ax4.set_xlabel("Lookback time", labelpad=10)
    ax4.set_xticks(lbt_pos)
    ax4.set_xticklabels(lbt_target_str)
    lns = lns1+lns2+lns3
    labs = [l.get_label() for l in lns]
    axes[0].legend(lns, labs, loc=0)
    # logend location codes:
    # 0 ~ 10 
    # best, ur, ul, lr, ll, r, cl, cr, lower c, upper c, center
    #
       
#    plt.show()
    plt.savefig(wdir + 'catalog/' + str(idgal).zfill(5) + '.png')
    plt.close()
cat = pickle.load(open('./10002/catalog_GM/catalog105.pickle', 'rb'))
gg = mpgs[0]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16656162780228684, 'fraction_numerical': 0.10533897542090573, 'mean_word_length': 2.831666091755778, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '667303', 'license': 'MIT', 'n_tokens_mistral': 5378, 'n_tokens_neox': 4778, 'n_words': 934, 'path': 'scripts/scripts_ipynb/0422_figs.ipynb', 'repo_name': 'Hoseung/pyRamAn', 'size': '11107'}