Source code for lsst.sims.maf.metrics.snNSNMetric

import numpy as np
from lsst.sims.maf.metrics import BaseMetric
from lsst.sims.maf.utils.snNSNUtils import Load_Reference, Telescope, LCfast
from lsst.sims.maf.utils.snNSNUtils import SN_Rate, CovColor
import pandas as pd
import numpy.lib.recfunctions as rf
import time
from scipy.interpolate import interp1d
import numpy.lib.recfunctions as nlr

__all__ = ['SNNSNMetric']

[docs]class SNNSNMetric(BaseMetric): """ Estimate (nSN,zlim) of type Ia supernovae. Parameters --------------- metricName : str, opt metric name (default : SNSNRMetric) mjdCol : str, opt mjd column name (default : observationStartMJD) RACol : str,opt Right Ascension column name (default : fieldRA) DecCol : str,opt Declinaison column name (default : fieldDec) filterCol : str,opt filter column name (default: filter) m5Col : str, opt five-sigma depth column name (default : fiveSigmaDepth) exptimeCol : str,opt exposure time column name (default : visitExposureTime) nightCol : str,opt night column name (default : night) obsidCol : str,opt observation id column name (default : observationId) nexpCol : str,opt number of exposure column name (default : numExposures) vistimeCol : str,opt visit time column name (default : visitTime) season : list,opt list of seasons to process (float)(default: -1 = all seasons) zmin : float,opt min redshift for the study (default: 0.0) zmax : float,opt max redshift for the study (default: 1.2) pixArea: float, opt pixel area (default: 9.6) verbose: bool,opt verbose mode (default: False) ploteffi: bool,opt to plot observing efficiencies vs z (default: False) n_bef: int, opt number of LC points LC before T0 (default:5) n_aft: int, opt number of LC points after T0 (default: 10) snr_min: float, opt minimal SNR of LC points (default: 5.0) n_phase_min: int, opt number of LC points with phase<= -5(default:1) n_phase_max: int, opt number of LC points with phase>= 20 (default: 1) """ def __init__(self, metricName='SNNSNMetric', mjdCol='observationStartMJD', RACol='fieldRA', DecCol='fieldDec', filterCol='filter', m5Col='fiveSigmaDepth', exptimeCol='visitExposureTime', nightCol='night', obsidCol='observationId', nexpCol='numExposures', vistimeCol='visitTime', season=[-1], zmin=0.0, zmax=1.2, pixArea=9.6, verbose=False, ploteffi=False, n_bef=4, n_aft=10, snr_min=5., n_phase_min=1, n_phase_max=1, templateDir=None, zlim_coeff=-1., **kwargs): self.mjdCol = mjdCol self.m5Col = m5Col self.filterCol = filterCol self.RACol = RACol self.DecCol = DecCol self.exptimeCol = exptimeCol self.seasonCol = 'season' self.nightCol = nightCol self.obsidCol = obsidCol self.nexpCol = nexpCol self.vistimeCol = vistimeCol self.pixArea = pixArea self.zlim_coeff = zlim_coeff cols = [self.nightCol, self.m5Col, self.filterCol, self.mjdCol, self.obsidCol, self.nexpCol, self.vistimeCol, self.exptimeCol] super(SNNSNMetric, self).__init__( col=cols, metricDtype='object', metricName=metricName, **kwargs) self.season = season # LC selection parameters self.n_bef = n_bef # nb points before peak self.n_aft = n_aft # nb points after peak self.snr_min = snr_min # SNR cut for points before/after peak self.n_phase_min = n_phase_min # nb of point with phase <=-5 self.n_phase_max = n_phase_max # nb of points with phase >=20 # loading reference LC files lc_reference = Load_Reference(templateDir=templateDir).ref self.lcFast = {} telescope = Telescope(airmass=1.2) for key, vals in lc_reference.items(): self.lcFast[key] = LCfast(vals, key[0], key[1], telescope, self.mjdCol, self.RACol, self.DecCol, self.filterCol, self.exptimeCol, self.m5Col, self.seasonCol, self.nexpCol, self.snr_min) # loading parameters self.zmin = zmin # zmin for the study self.zmax = zmax # zmax for the study self.zStep = 0.05 # zstep self.daymaxStep = 3. # daymax step self.min_rf_phase = -20. # min ref phase for LC points selection self.max_rf_phase = 40. # max ref phase for LC points selection self.min_rf_phase_qual = -15. # min ref phase for bounds effects self.max_rf_phase_qual = 25. # max ref phase for bounds effects # snrate self.rateSN = SN_Rate( min_rf_phase=self.min_rf_phase_qual, max_rf_phase=self.max_rf_phase_qual) # verbose mode - useful for debug and code performance estimation self.verbose = False self.ploteffi = False # supernovae parameters self.params = ['x0', 'x1', 'daymax', 'color'] # r = [(-1.0, -1.0)] self.bad = np.rec.fromrecords([(-1.0, -1.0)], names=['nSN', 'zlim']) # self.bad = {'nSN': -1.0, 'zlim': -1.0}
[docs] def run(self, dataSlice, slicePoint=None): """ run method of the metric Parameters --------------- dataSlice: array data to process """ idarray = None healpixID = -1 if slicePoint is not None: if 'nside' in slicePoint.keys(): import healpy as hp self.pixArea = hp.nside2pixarea( slicePoint['nside'], degrees=True) r = [] names = [] healpixID = hp.ang2pix( slicePoint['nside'], np.rad2deg(slicePoint['ra']), np.rad2deg(slicePoint['dec']), nest=True, lonlat=True) for kk, vv in slicePoint.items(): r.append(vv) names.append(kk) idarray = np.rec.fromrecords([r], names=names) else: idarray = np.rec.fromrecords([0., 0.], names=['RA', 'Dec']) # Two things to do: concatenate data (per band, night) and estimate seasons dataSlice = rf.drop_fields(dataSlice, ['season']) dataSlice = self.coadd(pd.DataFrame(dataSlice)) dataSlice = self.getseason(dataSlice) # get the seasons seasons = self.season # if seasons = -1: process the seasons seen in data if self.season == [-1]: seasons = np.unique(dataSlice[self.seasonCol]) # get redshift range for processing zRange = list(np.arange(self.zmin, self.zmax, self.zStep)) if zRange[0] < 1.e-6: zRange[0] = 0.01 self.zRange = np.unique(zRange) # season infos dfa = pd.DataFrame(np.copy(dataSlice)) season_info = dfa.groupby(['season']).apply( lambda x: self.seasonInfo(x)).reset_index() #print('season info', season_info) # select seasons of at least 30 days idx = season_info['season_length'] >= 60. season_info = season_info[idx] # check wether requested seasons can be processed test_season = season_info[season_info['season'].isin(seasons)] # if len(test_season) == 0: if test_season.empty: return nlr.merge_arrays([idarray, self.bad], flatten=True) else: seasons = test_season['season'] # print('test_seas', seasons) # print('hh', season_info) # get season length depending on the redshift dur_z = season_info.groupby(['season']).apply( lambda x: self.duration_z(x)).reset_index() # remove dur_z with negative season lengths idx = dur_z['season_length'] >= 10. dur_z = dur_z[idx] # generating simulation parameters gen_par = dur_z.groupby(['z', 'season']).apply( lambda x: self.calcDaymax(x)).reset_index() if gen_par.empty: return nlr.merge_arrays([idarray, self.bad], flatten=True) resdf = pd.DataFrame() for seas in seasons: vara_df = self.run_season( dataSlice, [seas], gen_par, dur_z) #print('res', seas, vara_df) if vara_df is not None: resdf = pd.concat((resdf, vara_df)) # final result: median zlim for a faint sn # and nsn_med for z<zlim if resdf.empty: return nlr.merge_arrays([idarray, self.bad], flatten=True) resdf = resdf.round({'zlim': 3, 'nsn_med': 3}) x1_ref = -2.0 color_ref = 0.2 idx = np.abs(resdf['x1']-x1_ref) < 1.e-5 idx &= np.abs(resdf['color']-color_ref) < 1.e-5 idx &= resdf['zlim'] > 0 if not resdf[idx].empty: zlim = resdf[idx]['zlim'].median() nSN = resdf[idx]['nsn_med'].sum() resd = np.rec.fromrecords([(nSN, zlim, healpixID)], names=[ 'nSN', 'zlim', 'healpixID']) res = nlr.merge_arrays([idarray, resd], flatten=True) else: res = nlr.merge_arrays([idarray, self.bad], flatten=True) return res
[docs] def reducenSN(self, metricVal): # At each slicepoint, return the sum nSN value. return np.sum(metricVal['nSN'])
[docs] def reducezlim(self, metricVal): # At each slicepoint, return the median zlim return np.median(metricVal['zlim'])
[docs] def coadd(self, data): """ Method to coadd data per band and per night Parameters --------------- data: pandas df of observations Returns ----------- coadded data (pandas df) """ keygroup = [self.filterCol, self.nightCol] data.sort_values(by=keygroup, ascending=[ True, True], inplace=True) coadd_df = data.groupby(keygroup).agg({self.nexpCol: ['sum'], self.vistimeCol: ['sum'], self.exptimeCol: ['sum'], self.mjdCol: ['mean'], self.RACol: ['mean'], self.DecCol: ['mean'], self.m5Col: ['mean']}).reset_index() coadd_df.columns = [self.filterCol, self.nightCol, self.nexpCol, self.vistimeCol, self.exptimeCol, self.mjdCol, self.RACol, self.DecCol, self.m5Col] coadd_df.loc[:, self.m5Col] += 1.25 * \ np.log10(coadd_df[self.vistimeCol]/30.) coadd_df.sort_values(by=[self.filterCol, self.nightCol], ascending=[ True, True], inplace=True) return coadd_df.to_records(index=False)
[docs] def getseason(self, obs, season_gap=80., mjdCol='observationStartMJD'): """ Method to estimate seasons Parameters -------------- obs: numpy array array of observations season_gap: float, opt minimal gap required to define a season (default: 80 days) mjdCol: str, opt col name for MJD infos (default: observationStartMJD) Returns ---------- original numpy array with seasonnumber appended """ # check wether season has already been estimated if 'season' in obs.dtype.names: return obs obs.sort(order=mjdCol) seasoncalc = np.ones(obs.size, dtype=int) if len(obs) > 1: diff = np.diff(obs[mjdCol]) flag = np.where(diff > season_gap)[0] if len(flag) > 0: for i, indx in enumerate(flag): seasoncalc[indx+1:] = i+2 obs = rf.append_fields(obs, 'season', seasoncalc) return obs
[docs] def seasonInfo(self, grp): """ Method to estimate seasonal info (cadence, season length, ...) Parameters -------------- grp: pandas df group Returns --------- pandas df with the cfollowing cols: """ df = pd.DataFrame([len(grp)], columns=['Nvisits']) df['MJD_min'] = grp[self.mjdCol].min() df['MJD_max'] = grp[self.mjdCol].max() df['season_length'] = df['MJD_max']-df['MJD_min'] df['cadence'] = 0. if len(grp) > 5: to = grp.groupby(['night'])[self.mjdCol].median().sort_values() df['cadence'] = np.mean(to.diff()) return df
[docs] def duration_z(self, grp): """ Method to estimate the season length vs redshift This is necessary to take into account boundary effects when estimating the number of SN that can be detected daymin, daymax = min and max MJD of a season T0_min(z) = daymin-(1+z)*min_rf_phase_qual T0_max(z) = daymax-(1+z)*max_rf_phase_qual season_length(z) = T0_max(z)-T0_min(z) Parameters -------------- grp: pandas df group data to process: season infos Returns ---------- pandas df with season_length, z, T0_min and T0_max cols """ daymin = grp['MJD_min'].values daymax = grp['MJD_max'].values dur_z = pd.DataFrame(self.zRange, columns=['z']) dur_z['T0_min'] = daymin-(1.+dur_z['z'])*self.min_rf_phase_qual dur_z['T0_max'] = daymax-(1.+dur_z['z'])*self.max_rf_phase_qual dur_z['season_length'] = dur_z['T0_max']-dur_z['T0_min'] return dur_z
[docs] def calcDaymax(self, grp): """ Method to estimate T0 (daymax) values for simulation. Parameters -------------- grp: group (pandas df sense) group of data to process with the following cols: T0_min: T0 min value (per season) T0_max: T0 max value (per season) Returns ---------- pandas df with daymax, min_rf_phase, max_rf_phase values """ T0_max = grp['T0_max'].values T0_min = grp['T0_min'].values num = (T0_max-T0_min)/self.daymaxStep if T0_max-T0_min > 10: df = pd.DataFrame(np.linspace( T0_min, T0_max, int(num)), columns=['daymax']) else: df = pd.DataFrame([-1], columns=['daymax']) df['min_rf_phase'] = self.min_rf_phase_qual df['max_rf_phase'] = self.max_rf_phase_qual return df
[docs] def run_season(self, dataSlice, season, gen_par, dura_z): """ Method to run on seasons Parameters -------------- dataSlice: numpy array, opt data to process (scheduler simulations) seasons: list(int) list of seasons to process Returns --------- effi_seasondf: pandas df efficiency curves zlimsdf: pandas df redshift limits and number of supernovae """ time_ref = time.time() if self.verbose: print('#### Processing season', season) groupnames = ['season', 'x1', 'color'] gen_p = gen_par[gen_par['season'].isin(season)] if gen_p.empty: if self.verbose: print('No generator parameter found') return None dur_z = dura_z[dura_z['season'].isin(season)] obs = pd.DataFrame(np.copy(dataSlice)) obs = obs[obs['season'].isin(season)] obs = obs.sort_values(by=['night']) #print('data here', obs.columns) #print(obs[['night', 'filter', 'observationStartMJD', 'fieldRA', 'fieldDec']]) """ import matplotlib.pyplot as plt plt.plot(dataSlice['fieldRA'], dataSlice['fieldDec'], 'ko') print('data', len(dataSlice)) plt.show() """ # simulate supernovae and lc if self.verbose: print("SN generation") print(season, obs) sn = self.genSN(obs.to_records( index=False), gen_p.to_records(index=False)) if self.verbose: idx = np.abs(sn['x1']+2) < 1.e-5 idx &= np.abs(sn['z']-0.2) < 1.e-5 sel = sn[idx] sel = sel.sort_values(by=['z', 'daymax']) print('sn and lc', len(sn), sel.columns, sel[['x1', 'color', 'z', 'daymax', 'Cov_colorcolor', 'n_bef', 'n_aft']]) # from these supernovae: estimate observation efficiency vs z effi_seasondf = self.effidf(sn) # zlims can only be estimated if efficiencies are ok idx = effi_seasondf['z'] <= 0.2 x1ref = -2.0 colorref = 0.2 idx &= np.abs(effi_seasondf['x1']-x1ref) < 1.e-5 idx &= np.abs(effi_seasondf['color']-colorref) < 1.e-5 sel = effi_seasondf[idx] if np.mean(sel['effi']) > 0.02: # estimate zlims zlimsdf = self.zlims(effi_seasondf, dur_z, groupnames) # estimate number of medium supernovae zlimsdf['nsn_med'], zlimsdf['var_nsn_med'] = zlimsdf.apply(lambda x: self.nsn_typedf( x, 0.0, 0.0, effi_seasondf, dur_z), axis=1, result_type='expand').T.values else: return None if self.verbose: print('#### SEASON processed', time.time()-time_ref, season) return zlimsdf
[docs] def genSN(self, obs, gen_par): """ Method to simulate LC and supernovae Parameters --------------- obs: numpy array array of observations(from scheduler) gen_par: numpy array array of parameters for simulation """ time_ref = time.time() # LC estimation sn_tot = pd.DataFrame() lc_tot = pd.DataFrame() for key, vals in self.lcFast.items(): time_refs = time.time() gen_par_cp = np.copy(gen_par) if key == (-2.0, 0.2): idx = gen_par_cp['z'] < 0.9 gen_par_cp = gen_par_cp[idx] lc = vals(obs, gen_par_cp, bands='grizy') if self.verbose: print('End of simulation', key, time.time()-time_refs) if self.verbose: print('End of simulation after concat', key, time.time()-time_refs) # estimate SN sn = pd.DataFrame() if len(lc) > 0: sn = self.process(pd.DataFrame(lc)) if self.verbose: print('End of supernova', time.time()-time_refs) if not sn.empty: sn_tot = pd.concat([sn_tot, pd.DataFrame(sn)], sort=False) if self.verbose: print('End of supernova - all', time.time()-time_ref) return sn_tot
[docs] def process(self, tab): """ Method to process LC: sigma_color estimation and LC selection Parameters -------------- tab: pandas df of LC points with the following cols: flux: flux fluxerr: flux error phase: phase snr_m5: Signal-to-Noise Ratio time: time(MJD) mag: magnitude m5: five-sigma depth magerr: magnitude error exposuretime: exposure time band: filter zp: zero-point season: season number healpixID: pixel ID pixRA: pixel RA pixDec: pixel Dec z: redshift daymax: T0 flux_e_sec: flux(in photoelec/sec) flux_5: 5-sigma flux(in photoelec/sec) F_x0x0, ...F_colorcolor: Fisher matrix elements x1: x1 SN color: color SN n_aft: number of LC points before daymax n_bef: number of LC points after daymax n_phmin: number of LC points with a phase < -5 n_phmax: number of LC points with a phase > 20 Returns ---------- """ # now groupby tab = tab.round({'daymax': 3, 'z': 3, 'x1': 2, 'color': 2}) groups = tab.groupby( ['daymax', 'season', 'z', 'x1', 'color']) tosum = [] for ia, vala in enumerate(self.params): for jb, valb in enumerate(self.params): if jb >= ia: tosum.append('F_'+vala+valb) tosum += ['n_aft', 'n_bef', 'n_phmin', 'n_phmax'] # apply the sum on the group sums = groups[tosum].sum().reset_index() # select LC according to the number of points bef/aft peak idx = sums['n_aft'] >= self.n_aft idx &= sums['n_bef'] >= self.n_bef idx &= sums['n_phmin'] >= self.n_phase_min idx &= sums['n_phmax'] >= self.n_phase_max if self.verbose: print('selection parameters', self.n_bef, self.n_aft, self.n_phase_min, self.n_phase_max) finalsn = pd.DataFrame() goodsn = pd.DataFrame(sums.loc[idx]) # estimate the color for SN that passed the selection cuts if len(goodsn) > 0: goodsn.loc[:, 'Cov_colorcolor'] = CovColor(goodsn).Cov_colorcolor finalsn = pd.concat([finalsn, goodsn], sort=False) badsn = pd.DataFrame(sums.loc[~idx]) # Supernovae that did not pass the cut have a sigma_color=10 if len(badsn) > 0: badsn.loc[:, 'Cov_colorcolor'] = 100. finalsn = pd.concat([finalsn, badsn], sort=False) return finalsn
[docs] def effidf(self, sn_tot, color_cut=0.04): """ Method estimating efficiency vs z for a sigma_color cut Parameters --------------- sn_tot: pandas df data used to estimate efficiencies color_cut: float, opt color selection cut(default: 0.04) Returns ---------- effi: pandas df with the following cols: season: season pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID x1: SN stretch color: SN color z: redshift effi: efficiency effi_err: efficiency error(binomial) """ sndf = pd.DataFrame(sn_tot) listNames = ['season', 'x1', 'color'] groups = sndf.groupby(listNames) # estimating efficiencies effi = groups[['Cov_colorcolor', 'z']].apply( lambda x: self.effiObsdf(x, color_cut)).reset_index(level=list(range(len(listNames)))) # this is to plot efficiencies and also sigma_color vs z if self.ploteffi: import matplotlib.pylab as plt fig, ax = plt.subplots() figb, axb = plt.subplots() self.plot(ax, effi, 'effi', 'effi_err', 'Observing Efficiencies', ls='-') sndf['sigma_color'] = np.sqrt(sndf['Cov_colorcolor']) self.plot(axb, sndf, 'sigma_color', None, '$\sigma_{color}$') # get efficiencies vs z plt.show() return effi
[docs] def plot(self, ax, effi, vary, erry=None, legy='', ls='None'): """ Simple method to plot vs z Parameters -------------- ax: axis where to plot effi: pandas df data to plot vary: str variable(column of effi) to plot erry: str, opt error on y-axis(default: None) legy: str, opt y-axis legend(default: '') """ grb = effi.groupby(['x1', 'color']) yerr = None for key, grp in grb: x1 = grp['x1'].unique()[0] color = grp['color'].unique()[0] if erry is not None: yerr = grp[erry] ax.errorbar(grp['z'], grp[vary], yerr=yerr, marker='o', label='(x1,color)=({},{})'.format(x1, color), lineStyle=ls) ftsize = 15 ax.set_xlabel('z', fontsize=ftsize) ax.set_ylabel(legy, fontsize=ftsize) ax.xaxis.set_tick_params(labelsize=ftsize) ax.yaxis.set_tick_params(labelsize=ftsize) ax.legend(fontsize=ftsize)
[docs] def effiObsdf(self, data, color_cut=0.04): """ Method to estimate observing efficiencies for supernovae Parameters -------------- data: pandas df - grp data to process Returns ---------- pandas df with the following cols: - cols used to make the group - effi, effi_err: observing efficiency and associated error """ # reference df to estimate efficiencies df = data.loc[lambda dfa: np.sqrt(dfa['Cov_colorcolor']) < 100000., :] # selection on sigma_c<= 0.04 df_sel = df.loc[lambda dfa: np.sqrt( dfa['Cov_colorcolor']) <= color_cut, :] # make groups (with z) group = df.groupby('z') group_sel = df_sel.groupby('z') # Take the ratio to get efficiencies rb = (group_sel.size()/group.size()) err = np.sqrt(rb*(1.-rb)/group.size()) var = rb*(1.-rb)*group.size() rb = rb.array err = err.array var = var.array rb[np.isnan(rb)] = 0. err[np.isnan(err)] = 0. var[np.isnan(var)] = 0. return pd.DataFrame({group.keys: list(group.groups.keys()), 'effi': rb, 'effi_err': err, 'effi_var': var})
[docs] def zlims(self, effi_seasondf, dur_z, groupnames): """ Method to estimate redshift limits Parameters -------------- effi_seasondf: pandas df season: season pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID x1: SN stretch color: SN color z: redshift effi: efficiency effi_err: efficiency error (binomial) dur_z: pandas df with the following cols: season: season z: redshift T0_min: min daymax T0_max: max daymax season_length: season length groupnames: list(str) list of columns to use to define the groups Returns ---------- pandas df with the following cols: pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID season: season number x1: SN stretch color: SN color zlim: redshift limit """ res = effi_seasondf.groupby(groupnames).apply( lambda x: self.zlimdf(x, dur_z)).reset_index(level=list(range(len(groupnames)))) return res
[docs] def zlimdf(self, grp, duration_z): """ Method to estimate redshift limits Parameters -------------- grp: pandas df group efficiencies to estimate redshift limits; columns: season: season pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID x1: SN stretch color: SN color z: redshift effi: efficiency effi_err: efficiency error (binomial) duration_z: pandas df with the following cols: season: season z: redshift T0_min: min daymax T0_max: max daymax season_length: season length Returns ---------- pandas df with the following cols: zlimit: redshift limit """ zlimit = 0.0 # z range for the study zplot = np.arange(self.zmin, self.zmax, 0.01) # print(grp['z'], grp['effi']) if len(grp['z']) <= 3: return pd.DataFrame({'zlim': [zlimit]}) # 'status': [int(status)]}) # interpolate efficiencies vs z effiInterp = interp1d( grp['z'], grp['effi'], kind='linear', bounds_error=False, fill_value=0.) if self.zlim_coeff < 0.: # in that case zlim is estimated from efficiencies # first step: identify redshift domain with efficiency decrease zlimit = self.zlim_from_effi(effiInterp, zplot) #status = self.status['ok'] else: zlimit = self.zlim_from_cumul( grp, duration_z, effiInterp, zplot) return pd.DataFrame({'zlim': [zlimit]})
# 'status': [int(status)]})
[docs] def zlim_from_cumul(self, grp, duration_z, effiInterp, zplot, rate='cte'): """ Method to estimate the redshift limit from the cumulative The redshift limit is estimated to be the z value corresponding to: frac(NSN(z<zlimit))=zlimi_coeff Parameters --------------- grp: pandas group data to process duration_z: array duration as a function of the redshift effiInterp: interp1d interpolator for efficiencies zplot: interp1d interpolator for redshift values rate: str, opt rate to estimate the number of SN to estimate zlimit rate = cte: rate independent of z rate = SN_rate: rate from SN_Rate class Returns ---------- zlimit: float the redshift limit """ if rate == 'SN_rate': # get rate season = np.median(grp['season']) idx = duration_z['season'] == season seas_duration_z = duration_z[idx] durinterp_z = interp1d( seas_duration_z['z'], seas_duration_z['season_length'], bounds_error=False, fill_value=0.) # estimate the rates and nsn vs z zz, rate, err_rate, nsn, err_nsn = self.rateSN(zmin=self.zmin, zmax=self.zmax, duration_z=durinterp_z, survey_area=self.pixArea) # rate interpolation rateInterp = interp1d(zz, nsn, kind='linear', bounds_error=False, fill_value=0) else: # this is for a rate z-independent nsn = np.ones(len(zplot)) rateInterp = interp1d(zplot, nsn, kind='linear', bounds_error=False, fill_value=0) nsn_cum = np.cumsum(effiInterp(zplot)*rateInterp(zplot)) if nsn_cum[-1] >= 1.e-5: nsn_cum_norm = nsn_cum/nsn_cum[-1] # normalize zlim = interp1d(nsn_cum_norm, zplot) zlimit = zlim(self.zlim_coeff).item() if self.ploteffi: self.plot_NSN_cumul(grp, nsn_cum_norm, zplot) else: zlimit = 0. return zlimit
[docs] def plot_NSN_cumul(self, grp, nsn_cum_norm, zplot): """ Method to plot the NSN cumulative vs redshift Parameters -------------- grp: pandas group data to process """ import matplotlib.pylab as plt fig, ax = plt.subplots() x1 = grp['x1'].unique()[0] color = grp['color'].unique()[0] ax.plot(zplot, nsn_cum_norm, label='(x1,color)=({},{})'.format(x1, color)) ftsize = 15 ax.set_ylabel('NSN ($z<$)', fontsize=ftsize) ax.set_xlabel('z', fontsize=ftsize) ax.xaxis.set_tick_params(labelsize=ftsize) ax.yaxis.set_tick_params(labelsize=ftsize) ax.set_xlim((0.0, 0.8)) ax.set_ylim((0.0, 1.05)) ax.plot([0., 1.2], [0.95, 0.95], ls='--', color='k') plt.legend(fontsize=ftsize) plt.show()
[docs] def zlim_from_effi(self, effiInterp, zplot): """ Method to estimate the redshift limit from efficiency curves The redshift limit is defined here as the redshift value beyond which efficiency decreases up to zero. Parameters --------------- effiInterp: interpolator use to get efficiencies zplot: numpy array redshift values Returns ----------- zlimit: float the redshift limit """ # get efficiencies effis = effiInterp(zplot) # select data with efficiency decrease idx = np.where(np.diff(effis) < -0.005)[0] # Bail out if there is no data if np.size(idx) == 0: return 0 z_effi = np.array(zplot[idx], dtype={ 'names': ['z'], 'formats': [float]}) # from this make some "z-periods" to avoid accidental zdecrease at low z z_gap = 0.05 seasoncalc = np.ones(z_effi.size, dtype=int) diffz = np.diff(z_effi['z']) flag = np.where(diffz > z_gap)[0] if len(flag) > 0: for i, indx in enumerate(flag): seasoncalc[indx+1:] = i+2 z_effi = rf.append_fields(z_effi, 'season', seasoncalc) # now take the highest season (end of the efficiency curve) idd = z_effi['season'] == np.max(z_effi['season']) zlimit = np.min(z_effi[idd]['z']) return zlimit
[docs] def zlimdf_deprecated(self, grp, duration_z): """ Method to estimate redshift limits Parameters -------------- grp: pandas df group efficiencies to estimate redshift limits; columns: season: season pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID x1: SN stretch color: SN color z: redshift effi: efficiency effi_err: efficiency error (binomial) duration_z: pandas df with the following cols: season: season z: redshift T0_min: min daymax T0_max: max daymax season_length: season length Returns ---------- pandas df with the following cols: zlimit: redshift limit """ zlimit = 0.0 # z range for the study zplot = list(np.arange(self.zmin, self.zmax, 0.001)) # print(grp['z'], grp['effi']) if len(grp['z']) <= 3: return pd.DataFrame({'zlim': [zlimit]}) # 'status': [int(status)]}) # interpolate efficiencies vs z effiInterp = interp1d( grp['z'], grp['effi'], kind='linear', bounds_error=False, fill_value=0.) # get rate season = np.median(grp['season']) idx = duration_z['season'] == season seas_duration_z = duration_z[idx] # print('hhh1', seas_duration_z) durinterp_z = interp1d( seas_duration_z['z'], seas_duration_z['season_length'], bounds_error=False, fill_value=0.) # estimate the rates and nsn vs z zz, rate, err_rate, nsn, err_nsn = self.rateSN(zmin=self.zmin, zmax=self.zmax, duration_z=durinterp_z, survey_area=self.pixArea) # rate interpolation rateInterp = interp1d(zz, nsn, kind='linear', bounds_error=False, fill_value=0) # estimate the cumulated number of SN vs z nsn_cum = np.cumsum(effiInterp(zplot)*rateInterp(zplot)) if nsn_cum[-1] >= 1.e-5: nsn_cum_norm = nsn_cum/nsn_cum[-1] # normalize zlim = interp1d(nsn_cum_norm, zplot) zlimit = zlim(0.95).item() # status = self.status['ok'] if self.ploteffi: import matplotlib.pylab as plt fig, ax = plt.subplots() x1 = grp['x1'].unique()[0] color = grp['color'].unique()[0] ax.plot(zplot, nsn_cum_norm, label='(x1,color)=({},{})'.format(x1, color)) ftsize = 15 ax.set_ylabel('NSN ($z<$)', fontsize=ftsize) ax.set_xlabel('z', fontsize=ftsize) ax.xaxis.set_tick_params(labelsize=ftsize) ax.yaxis.set_tick_params(labelsize=ftsize) ax.set_xlim((0.0, 1.2)) ax.set_ylim((0.0, 1.05)) ax.plot([0., 1.2], [0.95, 0.95], ls='--', color='k') plt.legend(fontsize=ftsize) plt.show() return pd.DataFrame({'zlim': [zlimit]})
# 'status': [int(status)]})
[docs] def nsn_typedf(self, grp, x1, color, effi_tot, duration_z, search=True): """ Method to estimate the number of supernovae for a given type of SN Parameters -------------- grp: pandas series with the following infos: pixRA: pixelRA pixDec: pixel Dec healpixID: pixel ID season: season x1: SN stretch color: SN color zlim: redshift limit x1, color: SN params to estimate the number effi_tot: pandas df with columns: season: season pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID x1: SN stretch color: SN color z: redshift effi: efficiency effi_err: efficiency error (binomial) duration_z: pandas df with the following cols: season: season z: redshift T0_min: min daymax T0_max: max daymax season_length: season length Returns ---------- nsn: float number of supernovae """ # get rate season = np.median(grp['season']) idx = duration_z['season'] == season seas_duration_z = duration_z[idx] # print('hhh2', seas_duration_z) durinterp_z = interp1d( seas_duration_z['z'], seas_duration_z['season_length'], bounds_error=False, fill_value=0.) if search: effisel = effi_tot.loc[lambda dfa: ( dfa['x1'] == x1) & (dfa['color'] == color), :] else: effisel = effi_tot nsn, var_nsn = self.nsn(effisel, grp['zlim'], durinterp_z) return (nsn, var_nsn)
[docs] def nsn(self, effi, zlim, duration_z): """ Method to estimate the number of supernovae Parameters -------------- effi: pandas df grp of efficiencies season: season pixRA: RA of the pixel pixDec: Dec of the pixel healpixID: pixel ID x1: SN stretch color: SN color z: redshift effi: efficiency effi_err: efficiency error (binomial) zlim: float redshift limit value duration_z: pandas df with the following cols: season: season z: redshift T0_min: min daymax T0_max: max daymax season_length: season length Returns ---------- nsn, var_nsn : float number of supernovae (and variance) with z<zlim """ if zlim < 1.e-3: return (-1.0, -1.0) dz = 0.001 zplot = list(np.arange(self.zmin, self.zmax, dz)) # interpolate efficiencies vs z effiInterp = interp1d( effi['z'], effi['effi'], kind='linear', bounds_error=False, fill_value=0.) # estimate the cumulated number of SN vs z zz, rate, err_rate, nsn, err_nsn = self.rateSN(zmin=self.zmin, zmax=self.zmax, dz=dz, duration_z=duration_z, survey_area=self.pixArea) nsn_cum = np.cumsum(effiInterp(zplot)*nsn) nsn_interp = interp1d(zplot, nsn_cum) nsn = nsn_interp(zlim).item() return [nsn, 0.0]