<jupyter_start><jupyter_text>Show the STARDUST2 resultsreferees want to be convinced that this is really 94% likely to be a SNIa.  Try the following: * show a sampling of type Ia light curves (should look like good fits)* show a sampling of 'best-fit' Type Ibc and II light curves (most are pretty bad)Next thing to try?* show some plots of chi2 or log likelihood spread over parameter space?<jupyter_code>classify.plot_maxlike_fit(results['salt2-extended'], templateset='snana')
plt.savefig("snRequiem_stardust_classify_Ia.png", overwrite=True)

sndat
from matplotlib import rcParams

rcParams['legend.shadow']=False
rcParams['legend.fontsize']=30
#rcParams['text.usetex']=True
classify.plot_fits(results, nshow=3, templateset='snana', zp=23.9, dpi=300)
fig = plt.gcf()
ax0, ax1 = fig.axes
ax0.set_xlim(57440, 57640)
ax1.set_xlim(57440, 57640)
                       
for ax, band in zip(fig.axes, ['f105w','f160w']):
    iband = np.where(sndat['band']==band)[0]
    ax.errorbar(sndat['time'][iband], sndat['flux'][iband], sndat['fluxerr'][iband],
                marker='D', color='k', ls=' ', ms=5)

ax0.set_ylabel(r'magnification-corrected flux density ($\mu$Jy)', usetex=False)
ax0.set_xlabel(r'observed time - model time delay', usetex=False)
ax1.set_xlabel(r'observed time - model time delay', usetex=False)


plt.savefig("snRequiem_stardust_classify_alltypes.pdf", 
            dpi=300, overwrite=True)
results['pIa'], results['pIbc'], results['pII']           
templatenamesdict = {'Ibc':[k for k in results.keys() 
                            if k in classify.SubClassDict_SNANA['ibc']],
                     'II':[k for k in results.keys() 
                           if k in classify.SubClassDict_SNANA['ii']]
                    }

# Set up the axes grid
fig = plt.figure(1, (6., 6.))
axgrid = ImageGrid(fig, 111, nrows_ncols=(3,2), 
                   axes_pad=0.1, aspect=False, share_all=True)


# get the list of names for all CC SN templates used
templatenamesdict = {'Ibc':[k for k in results.keys() 
                            if k in classify.SubClassDict_SNANA['ibc']],
                     'II':[k for k in results.keys() 
                           if k in classify.SubClassDict_SNANA['ii']]
                    }

# The actual SN data:
sn = results['salt2-extended']['sn']

# make a time array for plotting models
t0 =57500  # MJD near first data point
time = np.arange(t0-30, t0+120, 2)


# ====================
# TOP ROW OF AXES: 
axlist = [axgrid[0], axgrid[1]]

# Plot N=30 random samples from the Type Ia nested sampling sequence

# Grab 30 random samples from the Type Ia nested sampling sequence
# We draw from the second half of the sampling sequence to make sure we've settled in to near the likelihood peak
# vparamnames = results['salt2-extended']['res']['vparam_names']
samples = results['salt2-extended']['res']['samples']  # Each sample step gives ['t0', 'x0', 'x1', 'c']
isamples = randint(1200, 2012, size=30)

for band,ax in zip(['f105w','f160w'], axlist):

    # grab x1, c from the nested sampling
    modelfit = results['salt2-extended']['fit']
    for i in isamples:
        t0,x0,x1,c = samples[i]
        modelfit.set(t0=t0,x0=x0,x1=x1,c=c)
        
        # get the flux vs time for this model
        f = modelfit.bandflux(band, time, zp=23.9, zpsys='ab')

        # plot it
        ax.plot(time, f, color='k', alpha=0.1, lw=3)

    # plot the observations
    iband = np.where(sn['band']==band)
    ax.errorbar(sn['time'][iband], sn['flux'][iband], 
                sn['fluxerr'][iband], color='r', marker='d', ls=' ')

    if ax==axlist[-1]:
        ax.text(0.95,0.95, 'Ia', transform=ax.transAxes, ha='right', va='top')

# ====================
# Middle and bottom rows of axes
axesdict = {'Ibc':[axgrid[2], axgrid[3]],
            'II': [axgrid[4], axgrid[5]]}

# Plot max-likelihood light curves for all the Type Ibc and II templates
for sntype in ['Ibc', 'II']:
    tempnamelist = templatenamesdict[sntype]
    axlist = axesdict[sntype]
    for band,ax in zip(['f105w','f160w'], axlist):
        for tempname in tempnamelist:
            modelfit = results[tempname]['fit']

            # get the flux vs time for this model
            f = modelfit.bandflux(band, time, zp=23.9, zpsys='ab')
            # plot it
            ax.plot(time, f, color='k', alpha=0.1, lw=3)
            continue
        # plot the observations
        iband = np.where(sn['band']==band)
        ax.errorbar(sn['time'][iband], sn['flux'][iband], 
                    sn['fluxerr'][iband], color='r', marker='d', ls=' ')
            
        if ax==axlist[-1]:
            ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top')<jupyter_output><empty_output><jupyter_text>Light curves (Random samples from SALT2 or MLE across all CC templates)<jupyter_code>#####LIGHT CURVE WITH MLE FOR CC######
import plot_util
fig,axes,gs = plot_util.grid_plot(3,2,figsize=(14,10),sharex=True,sharey=True)
# The actual SN data:
sn = results['salt2-extended']['sn']

# make a time array for plotting models
t0_init =57500  # MJD near first data point
#time = np.arange(t0-30, t0+120, 2)
time = np.arange(-30, 120, 2)

# ====================
# TOP ROW OF AXES: 
#axlist = [axgrid[0], axgrid[1]]

# Plot N=30 random samples from the Type Ia nested sampling sequence

# Grab 30 random samples from the Type Ia nested sampling sequence
# We draw from the second half of the sampling sequence to make sure we've settled in to near the likelihood peak
# vparamnames = results['salt2-extended']['res']['vparam_names']
samples = results['salt2-extended']['res']['samples']  # Each sample step gives ['t0', 'x0', 'x1', 'c']
isamples = randint(0, samples.shape[0], size=30)


for band,ax in zip(['f105w','f160w'], axes[0][:]):

    # grab x1, c from the nested sampling
    modelfit = results['salt2-extended']['fit']
    for i in isamples:
        t0,x0,x1,c = samples[i]
        modelfit.set(t0=t0,x0=x0,x1=x1,c=c)
        
        # get the flux vs time for this model
        f = modelfit.bandflux(band, time+t0_init, zp=23.9, zpsys='ab')

        # plot it
        ax = plot_util.plot('plot',time, f, color='k', alpha=0.1, lw=3,ticklabelsize=20,
                           fontsize=20,ax=ax)
    if ax==axes[0][-1]:
        ax.text(0.95,0.95, 'Ia', transform=ax.transAxes, ha='right', va='top',fontsize=20)
    # plot the observations
    iband = np.where(sn['band']==band)
    ax.errorbar(sn['time'][iband]-t0_init, sn['flux'][iband], 
                sn['fluxerr'][iband], color='r', marker='d', ls=' ')
    ax.text(0.05,0.95, band.upper(), transform=ax.transAxes, ha='left', va='top',fontsize=20)

# ====================
# Middle and bottom rows of axes

i=0
# Plot max-likelihood light curves for all the Type Ibc and II templates
for sntype in ['Ibc', 'II']:
    tempnamelist = templatenamesdict[sntype]
    i+=1
    for band,ax in zip(['f105w','f160w'], axes[i][:]):
        for tempname in tempnamelist:
            modelfit = results[tempname]['fit']
            # get the flux vs time for this model
            f = modelfit.bandflux(band, time+t0_init, zp=23.9, zpsys='ab')
            # plot it
            ax = plot_util.plot('plot',time, f, color='k', alpha=0.1, lw=3,ticklabelsize=20,
                           fontsize=16,ax=ax)
            if ax==axes[i][-1]:
                ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top',fontsize=20)
            continue
            
        # plot the observations
        iband = np.where(sn['band']==band)
        ax.errorbar(sn['time'][iband]-t0_init, sn['flux'][iband], 
                   sn['fluxerr'][iband], color='r', marker='d', ls=' ')
        ax.text(0.05,0.95, band.upper(), transform=ax.transAxes, ha='left', va='top',fontsize=20)
    
    if ax==axes[-1]:
        #ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top')
        ax.set_xlabel('Observer Time (MJD-%i)'%t0_init,fontsize=20)
ticks=[]
for i in range(3):
    for tick in axes[i][-1].yaxis.get_major_ticks():
        ticks.append('')
    axes[i][-1].set_yticklabels(ticks)
    axes[i][-1].set_ylabel('')
    axes[i][0].set_ylim((0,1))
ticks=[]
for i in range(2):
    for j in range(2):
        for tick in axes[i][j].xaxis.get_major_ticks():
            ticks.append('')
        axes[i][j].set_xticklabels(ticks)
plt.xlabel('Observer Time (MJD-%i)'%t0_init,fontsize=20,labelpad=20)
plt.ylabel('Magnifcation-Corrected Flux Density ($\mu$Jy)',fontsize=20,labelpad=25)
plt.tight_layout()
plt.savefig('lightcurves.pdf',format='pdf')<jupyter_output>['time', 'band', 'flux', 'fluxerr', 'zp', 'zpsys', 'snid', 'img', 'mu', 'muerr', 'dt', 'index']<jupyter_text>Light curves with N random samples (from SALT2 or across all CC templates)<jupyter_code>#####LIGHT CURVE WITH N RANDOM SAMPLES######
import plot_util
fig,axes,gs = plot_util.grid_plot(3,2,figsize=(14,10),sharex=True,sharey=True)
# The actual SN data:
sn = results['salt2-extended']['sn']

# make a time array for plotting models
t0_init =57500  # MJD near first data point
#time = np.arange(t0-30, t0+120, 2)
time = np.arange(-30, 120, 2)
N_Samples = 100
# ====================
# TOP ROW OF AXES: 
#axlist = [axgrid[0], axgrid[1]]

# Plot N=30 random samples from the Type Ia nested sampling sequence

# Grab 30 random samples from the Type Ia nested sampling sequence
# We draw from the second half of the sampling sequence to make sure we've settled in to near the likelihood peak
# vparamnames = results['salt2-extended']['res']['vparam_names']
samples = results['salt2-extended']['res']['samples']  # Each sample step gives ['t0', 'x0', 'x1', 'c']
modelfit = results['salt2-extended']['fit']
weights = results['salt2-extended']['res']['weights']
finite_inds = np.where(np.isfinite(weights))[0]
samples = samples[finite_inds,:]
isamples = randint(0, samples.shape[0], size=N_Samples)

for band,ax in zip(['f105w','f160w'], axes[0][:]):

    # grab x1, c from the nested sampling
    
    for i in isamples:
        t0,x0,x1,c = samples[i]
        modelfit.set(t0=t0,x0=x0,x1=x1,c=c)

        # get the flux vs time for this model
        f = modelfit.bandflux(band, time+t0_init, zp=23.9, zpsys='ab')
        ax = plot_util.plot('plot',time, f, color='k', alpha=0.1, lw=3,ticklabelsize=20,
                           fontsize=20,ax=ax)
    if ax==axes[0][-1]:
        ax.text(0.95,0.95, 'Ia', transform=ax.transAxes, ha='right', va='top',fontsize=20)
    # plot the observations
    iband = np.where(sn['band']==band)
    ax.errorbar(sn['time'][iband]-t0_init, sn['flux'][iband], 
                sn['fluxerr'][iband], color='r', marker='d', ls=' ')
    ax.text(0.05,0.95, band.upper(), transform=ax.transAxes, ha='left', va='top',fontsize=20)

# ====================
# Middle and bottom rows of axes

i=0
# Plot max-likelihood light curves for all the Type Ibc and II templates
for sntype in ['Ibc', 'II']:
    tempnamelist = templatenamesdict[sntype]
    N_Per_Template = int(N_Samples/len(tempnamelist))
    i+=1
    for band,ax in zip(['f105w','f160w'], axes[i][:]):
        for tempname in tempnamelist:
            modelfit = results[tempname]['fit']
            samples = results[tempname]['res'].samples
            weights = results[tempname]['res'].weights
            finite_inds = np.where(np.isfinite(weights))[0]
            samples = samples[finite_inds,:]
            randInds = np.random.randint(0,samples.shape[0],N_Per_Template)
            for ind in randInds:
                modelfit.parameters[1:] = samples[ind,:]
                # get the flux vs time for this model
                f = modelfit.bandflux(band, time+t0_init, zp=23.9, zpsys='ab')
                # plot it
                ax = plot_util.plot('plot',time, f, color='k', alpha=0.05, lw=3,ticklabelsize=20,
                                   fontsize=20,ax=ax)
            if ax==axes[i][-1]:
                ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top',fontsize=20)
            continue
            
        # plot the observations
        iband = np.where(sn['band']==band)
        ax.errorbar(sn['time'][iband]-t0_init, sn['flux'][iband], 
                   sn['fluxerr'][iband], color='r', marker='d', ls=' ')
        ax.text(0.05,0.95, band.upper(), transform=ax.transAxes, ha='left', va='top',fontsize=20)
    
    if ax==axes[-1]:
        #ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top')
        ax.set_xlabel('Observer Time (MJD-%i)'%t0_init,fontsize=20)
ticks=[]
for i in range(3):
    for tick in axes[i][-1].yaxis.get_major_ticks():
        ticks.append('')
    axes[i][-1].set_yticklabels(ticks)
    axes[i][-1].set_ylabel('')
    axes[i][0].set_ylim((0,1))
ticks=[]
for i in range(2):
    for j in range(2):
        for tick in axes[i][j].xaxis.get_major_ticks():
            ticks.append('')
        axes[i][j].set_xticklabels(ticks)
plt.xlabel('Observer Time (MJD-%i)'%t0_init,fontsize=20,labelpad=20)
plt.ylabel('Magnifcation-Corrected Flux Density ($\mu$Jy)',fontsize=20,labelpad=25)
plt.tight_layout()
plt.savefig('lightcurves_NSAMPLES.pdf',format='pdf')<jupyter_output><empty_output><jupyter_text>Color curves (N from SALT2 or MLE across all CC templates)<jupyter_code>#####COLOR CURVE with MLE for CC#######
import plot_util
from copy import deepcopy
fig,axes,gs = plot_util.grid_plot(3,1,figsize=(10,14),sharex=True,sharey=True)
# The actual SN data:
sn=results['salt2-extended']['sn']
# temp=results['salt2-extended']['sn']
# temp.sort('index')
# iband=np.where(temp['band']=='f105w')[0]
# jband=np.where(temp['band']=='f160w')[0]
# sn=Table.read('mrg0138_phot.txt',format='ascii')
# sn['flux']=sn['psf_flux']
# sn['time']=temp['time']

# make a time array for plotting models
t0_init =57500  # MJD near first data point
time = np.arange(-30, 120, 2)


# ====================
# TOP ROW OF AXES: 
axlist = [axgrid[0], axgrid[1]]

# Plot N=30 random samples from the Type Ia nested sampling sequence

# Grab 30 random samples from the Type Ia nested sampling sequence
# We draw from the second half of the sampling sequence to make sure we've settled in to near the likelihood peak
# vparamnames = results['salt2-extended']['res']['vparam_names']
samples = results['salt2-extended']['res']['samples']  # Each sample step gives ['t0', 'x0', 'x1', 'c']
isamples = randint(0, samples.shape[0], size=30)


# grab x1, c from the nested sampling
modelfit = results['salt2-extended']['fit']
for i in isamples:
    t0,x0,x1,c = samples[i]
    modelfit.set(t0=t0,x0=x0,x1=x1,c=c)

    # get the flux vs time for this model
    f = modelfit.color('F105W','F160W', 'ab',time+t0_init)
    axes[0]=plot_util.plot('plot',time, f, color='k', alpha=0.1, lw=3,ticklabelsize=20,
                           fontsize=20,ax=axes[0])
    # plot it
    #axes[0].plot(time, f, color='k', alpha=0.1, lw=3)

# plot the observations
iband = np.where(sn['band']=='f105w')
jband = np.where(sn['band']=='f160w')
mags = -2.5*np.log10(sn['flux'])+23.9
magerrs=1.0857*sn['fluxerr']/sn['flux']

finalmagerrs = np.sqrt(magerrs[iband]**2+magerrs[jband]**2)

axes[0].errorbar(sn['time'][iband]-t0_init,mags[iband]-mags[jband], 
    finalmagerrs, color='r', marker='d', ls=' ')
axes[0].invert_yaxis()
axes[0].set_ylim((4,-1))
ticks=[]
for tick in axes[0].xaxis.get_major_ticks():
    ticks.append('')
axes[0].set_xticklabels(ticks)
axes[1].set_xticklabels(ticks)
axes[0].text(0.95,0.95, 'Ia', transform=axes[0].transAxes, ha='right', va='top',fontsize=20)
# ====================
# Middle and bottom rows of axes

i=1
# Plot max-likelihood light curves for all the Type Ibc and II templates
for sntype in ['Ibc', 'II']:
    tempnamelist = templatenamesdict[sntype]
    ax=axes[i]
    i+=1
    for tempname in tempnamelist:
        modelfit = results[tempname]['fit']

        # get the flux vs time for this model
        f = modelfit.color('F105W','F160W', 'ab',time+t0_init)
        # plot it
        ax = plot_util.plot('plot',time, f, color='k', alpha=0.1, lw=3,ticklabelsize=20,
                           fontsize=20,ax=ax)
        continue
    # plot the observations
    #iband = np.where(sn['band']==band)
    #ax.errorbar(sn['time'][iband], sn['flux'][iband], 
    #            sn['fluxerr'][iband], color='r', marker='d', ls=' ')
    
    ax.errorbar(sn['time'][iband]-t0_init,mags[iband]-mags[jband], 
        finalmagerrs, color='r', marker='d', ls=' ')
    ax.invert_yaxis()
    ax.set_ylim((4,-1))
    if ax==axes[-1]:
        #ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top')
        ax.set_xlabel('Observer Time (MJD-%i)'%t0_init,fontsize=20)
    ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top',fontsize=20)
plt.ylabel('F105W-F160W Color (Magnitudes)',fontsize=20,labelpad=25)
plt.tight_layout()
plt.savefig('colorcurves.pdf',format='pdf')<jupyter_output><ipython-input-33-70a667c8fe57>:61: UserWarning: FixedFormatter should only be used together with FixedLocator
  axes[0].set_xticklabels(ticks)
<ipython-input-33-70a667c8fe57>:62: UserWarning: FixedFormatter should only be used together with FixedLocator
  axes[1].set_xticklabels(ticks)
/Users/jpierel/miniconda3/envs/basic/lib/python3.8/site-packages/sncosmo/models.py:188: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)<jupyter_text>Color curves with N random samples (from SALT2 or across all CC templates)<jupyter_code>#####COLOR CURVE with N random samples#######
import plot_util
np.random.seed(3)
fig,axes,gs = plot_util.grid_plot(3,1,figsize=(10,14),sharex=True,sharey=True)
# The actual SN data:
sn = results['salt2-extended']['sn']
# temp=results['salt2-extended']['sn']
# temp.sort('index')
# iband=np.where(temp['band']=='f105w')[0]
# jband=np.where(temp['band']=='f160w')[0]
# sn=Table.read('mrg0138_phot.txt',format='ascii')
# sn['flux']=sn['psf_flux']
# sn['time']=temp['time']
# make a time array for plotting models
t0_init =57500  # MJD near first data point
time = np.arange(-30, 120, 2)
N_Samples = 250

# ====================
# TOP ROW OF AXES: 
axlist = [axgrid[0], axgrid[1]]

# Plot N=30 random samples from the Type Ia nested sampling sequence

# Grab 30 random samples from the Type Ia nested sampling sequence
# We draw from the second half of the sampling sequence to make sure we've settled in to near the likelihood peak
# vparamnames = results['salt2-extended']['res']['vparam_names']
samples = results['salt2-extended']['res']['samples']  # Each sample step gives ['t0', 'x0', 'x1', 'c']
weights = results['salt2-extended']['res']['weights']
finite_inds = np.where(np.isfinite(weights))[0]
samples = samples[finite_inds,:]
isamples = randint(0, samples.shape[0], size=N_Samples)


# grab x1, c from the nested sampling
modelfit = results['salt2-extended']['fit']
for i in isamples:
    t0,x0,x1,c = samples[i]
    modelfit.set(t0=t0,x0=x0,x1=x1,c=c)

    # get the flux vs time for this model
    f = modelfit.color('F105W','F160W', 'ab',time+t0_init)
    axes[0]=plot_util.plot('plot',time, f, color='k', alpha=0.05, lw=3,ticklabelsize=20,
                           fontsize=20,ax=axes[0])
    # plot it
    #axes[0].plot(time, f, color='k', alpha=0.1, lw=3)

# plot the observations
iband = np.where(sn['band']=='f105w')
jband = np.where(sn['band']=='f160w')
mags = -2.5*np.log10(sn['flux'])+23.9
magerrs=1.0857*sn['fluxerr']/sn['flux']

finalmagerrs = np.sqrt(magerrs[iband]**2+magerrs[jband]**2)

axes[0].errorbar(sn['time'][iband]-t0_init,mags[iband]-mags[jband], 
    finalmagerrs, color='r', marker='d', ls=' ')
axes[0].invert_yaxis()
axes[0].set_ylim((4,-1))
ticks=[]
for tick in axes[0].xaxis.get_major_ticks():
    ticks.append('')
axes[0].set_xticklabels(ticks)
axes[1].set_xticklabels(ticks)
axes[0].text(0.95,0.95, 'Ia', transform=axes[0].transAxes, ha='right', va='top',fontsize=20)
axes[0].set_xlim((time[0],time[-1]))
# ====================
# Middle and bottom rows of axes

i=1
# Plot max-likelihood light curves for all the Type Ibc and II templates
for sntype in ['Ibc', 'II']:
    tempnamelist = templatenamesdict[sntype]
    N_Per_Template = int(N_Samples/len(tempnamelist))
    ax=axes[i]
    i+=1
    for tempname in tempnamelist:
        modelfit = results[tempname]['fit']
        samples = results[tempname]['res'].samples
        weights = results[tempname]['res'].weights
        finite_inds = np.where(np.isfinite(weights))[0]
        samples = samples[finite_inds,:]
        randInds = np.random.randint(0,samples.shape[0],N_Per_Template)
        for ind in randInds:
            modelfit.parameters[1:] = samples[ind,:]
            # get the flux vs time for this model
            f = modelfit.color('F105W','F160W', 'ab',time+t0_init)
            # plot it
            ax = plot_util.plot('plot',time, f, color='k', alpha=0.05, lw=3,ticklabelsize=20,
                               fontsize=20,ax=ax)
    # plot the observations
    #iband = np.where(sn['band']==band)
    #ax.errorbar(sn['time'][iband], sn['flux'][iband], 
    #            sn['fluxerr'][iband], color='r', marker='d', ls=' ')
    
    ax.errorbar(sn['time'][iband]-t0_init,mags[iband]-mags[jband], 
        finalmagerrs, color='r', marker='d', ls=' ')
    ax.invert_yaxis()
    ax.set_ylim((4,-1))
    if ax==axes[-1]:
        #ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top')
        ax.set_xlabel('Observer Time (MJD-%i)'%t0_init,fontsize=20)
    ax.text(0.95,0.95, sntype, transform=ax.transAxes, ha='right', va='top',fontsize=20)
    ax.set_xlim((time[0],time[-1]))
plt.ylabel('F105W-F160W Color (Magnitudes)',fontsize=20,labelpad=25)
plt.tight_layout()
plt.savefig('colorcurves_NSAMPLES.pdf',format='pdf')<jupyter_output><ipython-input-32-effafa510726>:63: UserWarning: FixedFormatter should only be used together with FixedLocator
  axes[0].set_xticklabels(ticks)
<ipython-input-32-effafa510726>:64: UserWarning: FixedFormatter should only be used together with FixedLocator
  axes[1].set_xticklabels(ticks)
/Users/jpierel/miniconda3/envs/basic/lib/python3.8/site-packages/sncosmo/models.py:188: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
/Users/jpierel/miniconda3/envs/basic/lib/python3.8/site-packages/sncosmo/models.py:188: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
/Users/jpierel/miniconda3/envs/basic/lib/python3.8/site-packages/sncosmo/models.py:188: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
/Users/jpierel/miniconda3/envs/basic/lib/python3.8/site-packages/sncosmo/models.py:188: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)<jupyter_text>STOP<jupyter_code>def rescale_likelihood( a ):
	"""
	Rescale the likelihood array using 
	a Sorted Cumulative Sum function (e.g., for making contours)
	Construct an array "sumabove" such that the cell 
	at index i in sumabove is equal to the sum of all 
	cells from the input array "a" that have a
	cell value higher than a[i]
	"""
	# Collapse the array into 1 dimension
	sumabove = deepcopy(a).ravel()

	# Sort the raveled array by descending cell value
	iravelsorted = sumabove.argsort( axis=0 )[::-1]

	# Reassign each cell to be the cumulative sum of all
	# input array cells with a higher value :
	sumabove[iravelsorted] = sumabove[iravelsorted].cumsum()

	# Normalize by the max value in the array
	sumabove /= sumabove.max()

	# Now unravel back into shape of original array and return
	return( sumabove.reshape( a.shape ) )
import scipy,plot_util
plt.close()
fig,axes,gs = plot_util.grid_plot(3,1,figsize=(10,14),sharex=True,sharey=True)
t0_init =57500  # MJD near first data point
time = np.arange(-30, 120, 2)
i=1
grid_time=np.linspace(-30,120,1000)
grid_fluxes=np.linspace(0,1.5,1000)
band='f160w'
for sntype in ['Ibc', 'II']:
    ax=axes[i]
    i+=1
    tempnamelist = templatenamesdict[sntype]
    all_times=[]
    all_fluxes=[]
    all_probs=[]
    for tempname in tempnamelist:
        modelfit = results[tempname]['fit']
        # get the flux vs time for this model
        
        samples = results[tempname]['res'].samples
        weights = results[tempname]['res'].weights
        
        for k in np.random.randint(0,len(samples),size=100):
            modelfit.parameters[1:]=samples[k,:]
            f = modelfit.bandflux(band, grid_time+t0_init, zp=23.9, zpsys='ab')
            all_times=np.append(all_times,grid_time)
            all_fluxes=np.append(all_fluxes,f)
            all_probs=np.append(all_probs,[weights[i]]*len(grid_time))
    grid_probs=[]
   
    for j in range(len(grid_time)):
        
        temp=scipy.interpolate.griddata(np.array([all_times,all_fluxes]).T,
                                                             all_probs,
                    np.array([[grid_time[i]]*len(grid_fluxes),grid_fluxes]).T,
                                       fill_value=0)
        #temp = rescale_likelihood(temp)
        grid_probs.append(temp)
    
    #grid_probs=np.array(grid_probs).reshape([len(grid_probs),len(grid_probs)])
    grid_probs=rescale_likelihood(np.array(grid_probs))
    CS=ax.contourf(grid_time, 
                   grid_fluxes, 
                   np.array(grid_probs), colors='k',
                levels=[0,.68], alpha=.8, zorder=10)
    
    #ax.errorbar(sn['time'][iband]-t0_init,mags[iband]-mags[jband], 
    #    finalmagerrs, color='r', marker='d', ls=' ')
    plt.show()
    sys.exit()
print(grid_probs)
len(all_times)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17066656192309204, 'fraction_numerical': 0.03861110020032209, 'mean_word_length': 3.6707026233718585, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 26, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '29742', 'license': 'MIT', 'n_tokens_mistral': 10005, 'n_tokens_neox': 9194, 'n_words': 1979, 'path': 'Analysis/Classification/stardust2_lightcurve_classification.ipynb', 'repo_name': 'srodney/mrg0138_supernova-1', 'size': '25459'}