<jupyter_start><jupyter_text>Code used in the interpretation of AT2019dge.MCMC approach learnt from [this notebook](http://localhost:8889/notebooks/fit_early_power_law.ipynb).<jupyter_code>import numpy as np
import pandas as pd
import emcee
import corner
import time
import scipy.optimize as op
from allsn_info import get_at2019dge
from helper.models import model_arnett_Ltph
from multiprocessing import Pool
from helper import phys
from helper.mcmcfit import mylinear_fit<jupyter_output><empty_output><jupyter_text>Estimate Ejecta MassEstimatte $t_s$, assuming $t_{\rm rise}\approx 2\,$d.<jupyter_code>tmin_d = 2/(0.11*np.log(1 + 9*10) + 0.36)
tmin_d
tmax_d = 2/(0.11*np.log(1 + 9*0.01) + 0.36)
tmax_d
(tmin_d*86400)**2 * 1e+9 * phys.c / 0.1 / phys.sm
(tmax_d*86400)**2 * 1e+9 * phys.c / 0.1/ phys.sm
    Lpeak = 5 * 10**(42)
    tpeak = 2 * 86400
    ts = 8.8*86400
    beta = 4/3
    L0 = Lpeak * beta**2 * tpeak**2 / 2 / ts**2 / (1 - (1 + beta*tpeak/ts)*np.exp(-beta * tpeak / ts))
    epsilon_Ni = 3.9e+10 # erg / g / s
    M_Ni = L0 / (epsilon_Ni) / phys.sm
M_Ni<jupyter_output><empty_output><jupyter_text>Test radioactivity<jupyter_code>import matplotlib
import matplotlib.pyplot as plt
fs = 14
matplotlib.rcParams['font.size']=fs
ms = 8
matplotlib.rcParams['lines.markersize']=ms

%matplotlib notebook
def get_refineLbbaxis(ax):
    ax.set_ylabel('$L$ log'+r'$_{10}\rm(erg\,s^{-1})$')
    ax.set_xlabel('Time since explosion')
    ax.xaxis.set_major_locator(plt.MultipleLocator(10))
    ax.xaxis.set_minor_locator(plt.MultipleLocator(2))
    ax.yaxis.set_major_locator(plt.MultipleLocator(0.5))
    ax.yaxis.set_minor_locator(plt.MultipleLocator(0.1))
    ax.tick_params(direction='in', axis='both', which = 'both', top=True, right=True)
    ax.tick_params(which = 'major', length = 4)
    ax.tick_params(which = 'minor', length = 2)
def wygoda19_model(Mni, t, t0=30):
    """
    t in day
    Mni in the unit of Msun
    """
    Qgamma = 1e+43 * (6.45 * np.exp(-t/8.76) + 1.38 * np.exp(-t / 111.4)) * Mni
    Qpos = 4.64e+41 * (-np.exp(-t/8.76) + np.exp(-t / 111.4)) * Mni
    Q = Qgamma*(1 - np.exp(-(t0/t)**2)) + Qpos
    return Q
result = get_at2019dge()
tb0 = result['tb']
z = result['z']
data = pd.read_csv('../data/otherSN/Yao2020/bbdata.csv')
t_data = data['t'].values 
L_data = data['Lbb'].values
L_unc_data = data['Lbb_unc'].values
lgL_data = data['lgLbb'].values
lgL_unc_data = data['lgLbb_unc'].values
ind_rough = lgL_unc_data==0
# t_max = result["t_max"]
tb0 = tb0[tb0['filter'].values=='r']
tb0 = tb0[tb0.instrument!="P60+SEDM"]
tb0 = tb0[tb0.tmax_of+3.3 > max(t_data[~ind_rough])]
t_quasi = tb0["tmax_of"].values + 3.3 # time since explosion
Lquasi = tb0["Llambda"].values * tb0['wave'].values
Lquasi_unc = tb0["Llambda_unc"].values * tb0['wave'].values
lgLquasi = np.log10(Lquasi)
lgLquasi_unc = Lquasi_unc / Lquasi / np.log(10)
tnebula = np.linspace(25, 80, 200)
L30_002 = wygoda19_model(0.02, tnebula, t0=40)
lgL30_002 = np.log10(L30_002)
L30_001 = wygoda19_model(0.008, tnebula, t0=40)
lgL30_001 = np.log10(L30_001)
L30_0015 = wygoda19_model(0.004, tnebula, t0=40)
lgL30_0015 = np.log10(L30_0015)
tphot = np.linspace(0.1, 25, 100)

kappa = 0.06
Lambda = 6/5
Mej = 0.05 * phys.sm
beta = 13.8
Ek = 1e+50 # this is an order-or-magnitude estimate
print (np.sqrt(kappa / beta / phys.c) * (Lambda * Mej**3 / Ek)**(1/4) / 24 / 3600)

texp = -2
Lphot = model_arnett_Ltph(tphot-texp, taum_ = 5, Mni_ = 0.014)
lgLphot = np.log10(Lphot)<jupyter_output><empty_output><jupyter_text>Model the radiactive decay in the photospheric phase<jupyter_code>plt.figure(figsize=(6, 5))
ax = plt.subplot(111)
ax.errorbar(t_quasi, lgLquasi, lgLquasi_unc, fmt='--o', color = "grey", markerfacecolor='none', zorder = 1, markersize=7)
ax.errorbar(t_data[~ind_rough], lgL_data[~ind_rough], lgL_unc_data[~ind_rough], fmt='ok', zorder = 1, markersize=5)
ax.plot(t_data[ind_rough], lgL_data[ind_rough], 'ok', mfc='white', zorder = 1, markersize=5)
ax.plot(t_data, lgL_data, ':k', zorder = 1)

ax.plot(tphot-texp, lgLphot, 'c:', label = "Arnett model")
ax.plot(tnebula, lgL30_002, 'b--', label=r"$M_{\rm Ni}=0.02 M_\odot$")
ax.plot(tnebula, lgL30_001, 'b-', label=r"$M_{\rm Ni}=0.008 M_\odot$")
ax.plot(tnebula, lgL30_0015, 'b:', label=r"$M_{\rm Ni}=0.004 M_\odot$")
get_refineLbbaxis(ax)
plt.tight_layout(h_pad=0)
plt.legend(loc = "upper right", fontsize= fs)
ax.set_xlim(-5, 70)
ax.set_ylim(40., 43.)
ax.text(35, 40.1, 'nebular phase', fontsize = fs-2, color="k")
ax.text(-4, 40.1, 'photospheric phase', fontsize = fs-2, color="c")
plt.savefig('../paper/figures/Lbb.pdf')
plt.show()<jupyter_output><empty_output><jupyter_text>Nebula phaseUse MCMC fun to constrain $M_{\rm Ni}$ and $t_0$.<jupyter_code>x = np.hstack([t_data[[~ind_rough]], t_quasi])
y = np.hstack([lgL_data[~ind_rough], lgLquasi])
ey = np.hstack([lgL_unc_data[~ind_rough], lgLquasi_unc])
ix = x > 20
x = x[ix]
y = y[ix]
ey = ey[ix]
def w19_lnlike(theta, t, lgL, lgL_unc):
    lgMNi, t0 = theta
    model = wygoda19_model(10**lgMNi, t, t0)
    lgmodel = np.log10(model)
    chi2_term = -1/2*np.sum((lgL - lgmodel)**2/lgL_unc**2)
    error_term = np.sum(np.log(1/np.sqrt(2*np.pi*lgL_unc**2)))
    ln_l = chi2_term + error_term
    return ln_l

w19_nll = lambda *args: -w19_lnlike(*args)

#Define priors on parameters 
def w19_lnprior(theta):
    lgMNi, t0 = theta
    if (-4 < lgMNi < -1 and 5 < t0 < 100):
        return 0.0
    return -np.inf


def w19_lnprob(theta, t, lgL, lgL_unc):
    lp = w19_lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + w19_lnlike(theta, t, lgL, lgL_unc)
ml_guess = np.array([-2, 30])
ndim = len(ml_guess)
nwalkers = 100
nfac = [1e-3]*ndim
max_samples = 300
check_tau = 20

#initial position of walkers
pos = [ml_guess + ml_guess * nfac * np.random.randn(ndim) for i in range(nwalkers)]
with Pool(8) as pool:
    sampler = emcee.EnsembleSampler(nwalkers, ndim, w19_lnprob, 
                                    args=(x, y, ey),
                                    pool=pool)
    index = 0
    autocorr = np.empty(max_samples)
    old_tau = np.inf
    for sample in sampler.sample(pos, iterations=max_samples, thin_by=250, progress=True):
        if sampler.iteration % check_tau:
            continue
        
        tstart = time.time()
        tau = sampler.get_autocorr_time(tol=0)
        tend = time.time()
        autocorr[index] = np.mean(tau)
        index += 1
        steps_so_far = index*check_tau
        print('''After {:d} steps, 
autocorrelation takes {:.3f} s ({} total FFTs)                
acceptance fraction = {:.4f}, and
tau = {}'''.format(steps_so_far, 
                   tend-tstart, nwalkers*ndim,
                   np.mean(sampler.acceptance_fraction), 
                   tau))

        # Check convergence
        converged = np.all(tau * 100 < sampler.iteration)
        converged &= np.all(np.abs(old_tau - tau) / tau < 0.01)
        if converged:
            break
        old_tau = tau
converged
#define funtion to make corner plot
def makeCorner(sampler, nburn, paramsNames, quantiles=[0.16, 0.5, 0.84], truths=[]):
    samples = sampler.get_chain(discard=nburn, flat=True)
    if len(truths) > 0:
        f = corner.corner(samples, labels = paramsNames, quantiles = quantiles, 
                          truths=truths, plot_datapoints=False, title_kwargs = {"fontsize": fs})
    else:
        f = corner.corner(samples, labels = paramsNames, quantiles = quantiles, 
                          show_titles=True, plot_datapoints=False, title_kwargs = {"fontsize": fs})
    
#define function to plot walker chains  
def plotChains(sampler, nburn, paramsNames, nplot=nwalkers):
    Nparams = len(paramsNames)
    fig, ax = plt.subplots(Nparams+1,1, figsize = (8,2*(Nparams+1)), sharex = True)
    fig.subplots_adjust(hspace = 0)
    ax[0].set_title('Chains', fontsize=fs)
    xplot = np.arange(sampler.get_chain().shape[0])

    selected_walkers = np.random.choice(range(sampler.get_chain().shape[1]), nplot, replace=False)
    for i,p in enumerate(paramsNames):
        for w in selected_walkers:
            burn = ax[i].plot(xplot[:nburn], sampler.get_chain()[:nburn,w,i], 
                              alpha = 0.4, lw = 0.7, zorder = 1)
            ax[i].plot(xplot[nburn:], sampler.get_chain(discard=nburn)[:,w,i], 
                       color=burn[0].get_color(), alpha = 0.8, lw = 0.7, zorder = 1)
            
            ax[i].set_ylabel(p)
            if i==Nparams-1:
                ax[i+1].plot(xplot[:nburn], sampler.get_log_prob()[:nburn,w], 
                             color=burn[0].get_color(), alpha = 0.4, lw = 0.7, zorder = 1)
                ax[i+1].plot(xplot[nburn:], sampler.get_log_prob(discard=nburn)[:,w], 
                             color=burn[0].get_color(), alpha = 0.8, lw = 0.7, zorder = 1)
                ax[i+1].set_ylabel('ln P')
            
    return ax
gr_paramsNames=['lg' +r'$M_\mathrm{Ni}$', 
                r"$t_0$"]

plotChains(sampler, 25, gr_paramsNames, nplot=35)
plt.tight_layout()
makeCorner(sampler, 25, gr_paramsNames)
print (10**(-2.09))
print (10**(-2.09+0.06)-10**(-2.09))
print (10**(-2.09)-10**(-2.09-0.04))
tneb = np.linspace(20, 60, 100)
Lneb = wygoda19_model(0.008, tneb, t0=37)<jupyter_output><empty_output><jupyter_text>Photospheric phase<jupyter_code>from helper.models import arnett_nll, arnett_lnprob
x = np.hstack([t_data[[~ind_rough]], t_quasi])
y = np.hstack([lgL_data[~ind_rough], lgLquasi])
ey = np.hstack([lgL_unc_data[~ind_rough], lgLquasi_unc])
ix = (x < 20)&(x > 10)
x = x[ix]
y = y[ix]
ey = ey[ix]
ml_guess = np.array([5, -1.77, 0])
ndim = len(ml_guess)
nwalkers = 100
nfac = [1e-3]*ndim
max_samples = 300
check_tau = 20

#initial position of walkers
pos = [ml_guess + ml_guess * nfac * np.random.randn(ndim) for i in range(nwalkers)]
plt.figure()
plt.errorbar(x, y, ey)
Lmodel = model_arnett_Ltph(x, taum_ = 5, Mni_ = 0.017)
lgLmodel = np.log10(Lmodel)
plt.plot(x, lgLmodel)
with Pool(8) as pool:
    sampler = emcee.EnsembleSampler(nwalkers, ndim, arnett_lnprob, 
                                    args=(x, y, ey),
                                    pool=pool)
    index = 0
    autocorr = np.empty(max_samples)
    old_tau = np.inf
    for sample in sampler.sample(pos, iterations=max_samples, thin_by=250, progress=True):
        if sampler.iteration % check_tau:
            continue
        
        tstart = time.time()
        tau = sampler.get_autocorr_time(tol=0)
        tend = time.time()
        autocorr[index] = np.mean(tau)
        index += 1
        steps_so_far = index*check_tau
        print('''After {:d} steps, 
autocorrelation takes {:.3f} s ({} total FFTs)                
acceptance fraction = {:.4f}, and
tau = {}'''.format(steps_so_far, 
                   tend-tstart, nwalkers*ndim,
                   np.mean(sampler.acceptance_fraction), 
                   tau))

        # Check convergence
        converged = np.all(tau * 100 < sampler.iteration)
        converged &= np.all(np.abs(old_tau - tau) / tau < 0.01)
        if converged:
            break
        old_tau = tau<jupyter_output>7%|▋         | 4999/75000 [04:54<1:11:38, 16.28it/s]/anaconda3/lib/python3.6/site-packages/emcee/autocorr.py:38: RuntimeWarning: invalid value encountered in true_divide
  acf /= acf[0]
/anaconda3/lib/python3.6/site-packages/emcee/autocorr.py:43: RuntimeWarning: invalid value encountered in less
  m = np.arange(len(taus)) < c * taus
/anaconda3/lib/python3.6/site-packages/emcee/autocorr.py:101: RuntimeWarning: invalid value encountered in greater
  flag = tol * tau_est > n_t
/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:27: RuntimeWarning: invalid value encountered in less
/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in less
  7%|▋         | 5003/75000 [04:55<1:14:11, 15.73it/s]<jupyter_text>Let us subtract a radioactivity model from the bolometric light curve.<jupyter_code>plt.figure(figsize=(6,5))
plt.errorbar(tphase, L_data, L_unc_data, fmt = ".r")
plt.errorbar(t_quasi, Lquasi, Lquasi_unc, fmt='r:o', markerfacecolor='none', zorder = 1, markersize=7)
plt.plot(tneb, Lneb, 'k-')
plt.plot(tphot, Lphot)
plt.ylim(1e+40, 1e+43)
plt.semilogy()
L_resi = L_data - model_arnett_Ltph(tphase, taum_ = 1.5, Mni_ = 0.008)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1559117910326308, 'fraction_numerical': 0.05362519326226707, 'mean_word_length': 3.3335684062059237, 'pattern_counts': {'":': 2, '<': 40, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '590782', 'license': 'MIT', 'n_tokens_mistral': 5440, 'n_tokens_neox': 4943, 'n_words': 951, 'path': 'playground/.ipynb_checkpoints/radioactivity-Copy1-checkpoint.ipynb', 'repo_name': 'yaoyuhan/AT2019dge', 'size': '12289'}