From 9f9e9836a22a07b9b3b5dbaffc7dbd8b01cc983e Mon Sep 17 00:00:00 2001 From: Ulrich Feindt Date: Fri, 15 Feb 2019 21:55:36 +0100 Subject: [PATCH 1/4] fixed nugent template late-time issues; added option to not use MW dust --- setup.py | 6 ++--- simsurvey/__init__.py | 2 +- simsurvey/simultarget.py | 54 +++++++++++++++++++++++++++++++--------- 3 files changed, 46 insertions(+), 16 deletions(-) diff --git a/setup.py b/setup.py index 8427041..cf990d6 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # -# Copyright (C) 2015-18 simsurvey Developers +# Copyright (C) 2015-19 simsurvey Developers DESCRIPTION = "simsurvey: Basic tools for simulating astronomical surveys" LONG_DESCRIPTION = """\ @@ -13,8 +13,8 @@ MAINTAINER_EMAIL = 'ulrich.feindt@fysik.su.se' URL = 'https://github.com/ufeindt/simsurvey/' LICENSE = 'BSD (3-clause)' -DOWNLOAD_URL = 'https://github.com/ufeindt/simsurvey/tarball/0.4.4' -VERSION = '0.4.4' +DOWNLOAD_URL = 'https://github.com/ufeindt/simsurvey/tarball/0.4.5' +VERSION = '0.4.5' try: from setuptools import setup, find_packages diff --git a/simsurvey/__init__.py b/simsurvey/__init__.py index 1792bcb..d49cab7 100644 --- a/simsurvey/__init__.py +++ b/simsurvey/__init__.py @@ -1,7 +1,7 @@ """ Module base on astrobject, sncosmo and astropy to prepare a future survey. """ -__version__ = "0.4.4" +__version__ = "0.4.5" from .simultarget import * from .simulsurvey import * diff --git a/simsurvey/simultarget.py b/simsurvey/simultarget.py index 9ad831b..0c05ad4 100644 --- a/simsurvey/simultarget.py +++ b/simsurvey/simultarget.py @@ -58,7 +58,7 @@ def generate_transients(zrange,**kwargs): """ return get_transient_generator(zrange,**kwargs).transients -def generate_lightcurves(zrange, obs, **kwargs): +def generate_lightcurves(zrange, obs, trim_observations=True, **kwargs): """ This module calls get_transient_generator to create the TransientGenerator object and then generates lightcurves based @@ -68,7 +68,7 @@ def generate_lightcurves(zrange, obs, **kwargs): """ tr = get_transient_generator(zrange,**kwargs) - return tr.get_lightcurves(obs) + return tr.get_lightcurves(obs, trim_observations=trim_observations) ####################################### @@ -83,7 +83,7 @@ class TransientGenerator( BaseObject ): PROPERTIES = ["transient_coverage", "event_coverage"] - SIDE_PROPERTIES = ["sfd98_dir", "ratefunc", "model", "err_mwebv"] + SIDE_PROPERTIES = ["sfd98_dir", "ratefunc", "model", "err_mwebv", "apply_mwebv"] DERIVED_PROPERTIES = ["simul_parameters", "mwebv", "mwebv_sfd98", "has_mwebv_sfd98", "lightcurve_parameters"] @@ -99,7 +99,7 @@ def __init__(self, empty=False, **kwargs): def create(self, zrange=(0.0, 0.2), ratekind="basic", ratefunc=None, ntransient=None, transient=None, template=None, load=False, mjd_range=(57754.0,58849.0), - ra_range=(0,360), dec_range=(-90,90), + ra_range=(0,360), dec_range=(-90,90), apply_mwebv=True, mw_exclusion=0, sfd98_dir=None, transientprop=None, err_mwebv=0.01): """ """ @@ -114,7 +114,8 @@ def create(self, zrange=(0.0, 0.2), ratekind="basic", ratefunc=None, transientprop = {} self.set_sfd98_dir(sfd98_dir) - + self.set_apply_mwebv(apply_mwebv) + if not load: self.set_event_parameters(update=False, **{"ra_range":ra_range, "dec_range":dec_range, @@ -292,10 +293,14 @@ def get_lightcurve_full_param(self, *args, **kwargs): full_out = kwargs.get("full_out", True) for i in range(*range_args(self.ntransient, *args)): out = dict(z=self.zcmb[i], t0=self.mjd[i], - mwebv=(self.mwebv[i] - if self.has_mwebv_sfd98 else 0), **{p: v[i] for p, v in self.lightcurve.items()}) + if self.apply_mwebv: + if self.has_mwebv_sfd98: + out["mwebv"] = self.mwebv[i] + else: + out["mwebv"] = 0 + if full_out: out["ra"] = self.ra[i] out["dec"] = self.dec[i] @@ -306,11 +311,12 @@ def get_lightcurve_full_param(self, *args, **kwargs): yield out - def get_lightcurves(self, obs, **kwargs): + def get_lightcurves(self, obs, trim_observations=True, **kwargs): """Realize lightcurves based on the randomized lightcurve parameters and a single set of observations""" params = self.get_lightcurve_full_param(full_out=False) - return sncosmo.realize_lcs(obs, self.model, params, **kwargs) + return sncosmo.realize_lcs(obs, self.model, params, + trim_observations=trim_observations, **kwargs) # --------------------------- # # - Plots Methods - # @@ -718,11 +724,21 @@ def set_model(self, model): # if model.__class__ is not sncosmo.models.Model: # raise TypeError("model must be sncosmo.model.Model") - if "mwebv" not in model.param_names: + if "mwebv" not in model.param_names and self.apply_mwebv: model.add_effect(sncosmo.CCM89Dust(), 'mw', 'obs') self._side_properties["model"] = model + @property + def apply_mwebv(self): + """Option to switch of MW dust completely (need at wavelengths > 3.3 micron)""" + return self._side_properties["apply_mwebv"] + + def set_apply_mwebv(self, apply_mwebv): + """ + """ + self._side_properties["apply_mwebv"] = apply_mwebv + @property def err_mwebv(self): """Assumed error of dustmap; will be applied to lightcurve creation""" @@ -1015,7 +1031,14 @@ def model_Ibc_snana(self): def model_IIn_nugent(self): """ """ - return sncosmo.Model(source='nugent-sn2n', + sncosmo.get_source('nugent-sn2n', version='2.1') + p, w, f = sncosmo.read_griddata_ascii( + os.path.join(sncosmo.builtins.get_cache_dir(), + 'sncosmo/models/nugent/sn2n_flux.v2.1.dat') + ) + mask = (p < 150) + + return sncosmo.Model(source=sncosmo.TimeSeriesSource(p[mask], w, f[mask]), effects=[sncosmo.CCM89Dust()], effect_names=['host'], effect_frames=['rest']) @@ -1030,7 +1053,14 @@ def model_IIn_snana(self): def model_IIP_nugent(self): """ """ - return sncosmo.Model(source='nugent-sn2p', + sncosmo.get_source('nugent-sn2p', version='1.2') + p, w, f = sncosmo.read_griddata_ascii( + os.path.join(sncosmo.builtins.get_cache_dir(), + 'sncosmo/models/nugent/sn2p_flux.v1.2.dat') + ) + mask = (p < 130) + + return sncosmo.Model(source=sncosmo.TimeSeriesSource(p[mask], w, f[mask]), effects=[sncosmo.CCM89Dust()], effect_names=['host'], effect_frames=['rest']) From 02f6b16310c187e22abc3eccb2fa6aed3e9251e3 Mon Sep 17 00:00:00 2001 From: Ulrich Feindt Date: Tue, 2 Apr 2019 15:33:32 +0200 Subject: [PATCH 2/4] untested v0.5.0 --- simsurvey/simulsurvey.py | 66 +++++++++++++++++++++++----------------- simsurvey/simultarget.py | 16 +++------- 2 files changed, 43 insertions(+), 39 deletions(-) diff --git a/simsurvey/simulsurvey.py b/simsurvey/simulsurvey.py index 2911ac1..908c66f 100644 --- a/simsurvey/simulsurvey.py +++ b/simsurvey/simulsurvey.py @@ -552,8 +552,9 @@ class SurveyPlan( BaseObject ): DERIVED_PROPERTIES = [] def __init__(self, time=None, ra=None, dec=None, band=None, skynoise=None, - obs_field=None, zp=None, comment=None, width=7.295, height=7.465, - fields=None, empty=False, load_opsim=None, **kwargs): + obs_field=None, obs_ccd,= None, zp=None, comment=None, + width=7.295, height=7.465, fields=None, empty=False, + load_opsim=None, **kwargs): """ Parameters: ---------- @@ -565,12 +566,12 @@ def __init__(self, time=None, ra=None, dec=None, band=None, skynoise=None, return self.create(time=time,ra=ra,dec=dec,band=band,skynoise=skynoise, - obs_field=obs_field, zp=zp, comment=comment, + obs_field=obs_field, obs_ccd=obs_ccd, zp=zp, comment=comment, width=width, height=height, fields=fields, load_opsim=load_opsim, **kwargs) def create(self, time=None, ra=None, dec=None, band=None, skynoise=None, - obs_field=None, zp=None, comment=None, + obs_field=None, obs_ccd=None, zp=None, comment=None, width=7.295, height=7.465, fields=None, load_opsim=None, **kwargs): """ @@ -584,7 +585,8 @@ def create(self, time=None, ra=None, dec=None, band=None, skynoise=None, if load_opsim is None: self.add_observation(time, band, skynoise, ra=ra, dec=dec, - zp=zp, comment=comment, field=obs_field) + zp=zp, comment=comment, field=obs_field, + ccd=obs_ccd) else: self.load_opsim(load_opsim, **kwargs) @@ -631,12 +633,14 @@ def add_observation(self, time, band, skynoise, ra=None, dec=None, field=None, if zp is None: zp = np.array([np.nan for r in ra]) + if ccd is None: + ccd = np.array([np.nan for r in ra]) if comment is None: comment = np.array(['' for r in ra]) - new_obs = Table(data=[time, band, zp, skynoise, ra, dec, field, comment], + new_obs = Table(data=[time, band, zp, skynoise, ra, dec, field, ccd, comment], names=['time', 'band', 'zp', 'skynoise', - 'RA', 'Dec', 'field', 'comment']) + 'RA', 'Dec', 'field', 'ccd', 'comment']) if self._properties['pointings'] is None: self._properties['pointings'] = new_obs @@ -769,13 +773,20 @@ def get_non_field_obs(self, ra, dec, progress_bar=False, notebook=False): ccd = [np.array([], dtype=int) for r in ra] if single_coord: - if tmp['field']: + add = False + if (tmp['field'] and + (np.isnan(obs['ccd']) or + (not np.isnan(obs['ccd']) and tmp['ccd'] == obs['ccd']))): observed = True out = np.append(out, [k]) if ccd is not None: ccd = np.append(ccd, [tmp['ccd']]) + else: - for l in np.where(tmp['field'])[0]: + mask = tmp['field'] + if not np.isnan(obs['ccd']): + mask = mask & (tmp['ccd'] == obs['ccd']) + for l in np.where(mask)[0]: observed = True out[l] = np.append(out[l], [k]) if ccd is not None: @@ -804,7 +815,10 @@ def observed_on(self, fields=None, ccds=None, if fields is not None: for k, l in enumerate(fields): - mask = (self.pointings['field'] == l) + if ccds is None: + mask = (self.pointings['field'] == l) + else: + mask = (self.pointings['field'] == l) & self.pointing['ccd'] == ccds[k]) out['time'].extend(self.pointings['time'][mask].quantity.value) out['band'].extend(self.pointings['band'][mask]) out['zp'].extend(self.pointings['zp'][mask]) @@ -889,10 +903,10 @@ class LightcurveCollection( BaseObject ): __nature__ = "LightcurveCollection" PROPERTIES = ['lcs', 'meta', 'meta_rejected'] - SIDE_PROPERTIES = ['threshold', 'n_samenight', 'p_bins'] + SIDE_PROPERTIES = ['threshold', 'n_det', 'p_bins'] DERIVED_PROPERTIES = ['stats'] - def __init__(self, threshold=5., n_samenight=2, + def __init__(self, threshold=5., n_det=2, p_bins=np.arange(-30, 71, 5), empty=False, **kwargs): """ @@ -903,7 +917,7 @@ def __init__(self, threshold=5., n_samenight=2, """ self.__build__() self.set_threshold(threshold) - self.set_n_samenight(n_samenight) + self.set_n_det(n_det) self.set_p_bins(p_bins) self._prep_stats_() @@ -973,7 +987,7 @@ def save_tar(self, filename, notebook=False): conffile = os.path.join(file_base, 'config') f = open(conffile, 'w') print('threshold: %.3f' % self.threshold, file=f) - print('n_samenight: %i' % self.n_samenight, file=f) + print('n_det: %i' % self.n_det, file=f) print('p_bins:', self.p_bins, file=f) f.close() tar.add(conffile) @@ -1142,7 +1156,7 @@ def _add_lc_stats_(self, lc): """ """ p0, p1, dt = get_p_det_last(lc, thr=self.threshold, - n_samenight=self.n_samenight) + n_det=self.n_det) if p0 < 1e11 and p1 > -1e11: self._derived_properties['stats']['p_det'] = np.append( self._derived_properties['stats']['p_det'], p0 @@ -1276,14 +1290,14 @@ def set_threshold(self, threshold): self._side_properties["threshold"] = threshold @property - def n_samenight(self): + def n_det(self): """""" - return self._side_properties["n_samenight"] + return self._side_properties["n_det"] - def set_n_samenight(self, n_samenight): + def set_n_det(self, n_det): """ """ - self._side_properties["n_samenight"] = n_samenight + self._side_properties["n_det"] = n_det @property def p_bins(self): @@ -1317,18 +1331,14 @@ def get_tab_p_binned(self, key): # ========================= # # = Lightcurve statistics = # # ========================= # -def get_p_det_last(lc, thr=5., n_samenight=2): +def get_p_det_last(lc, thr=5., n_det=2): """ """ mask_det = lc['flux']/lc['fluxerr'] > thr - idx_nights = identify_nights(lc['time']) - - if np.sum(mask_det) > 1: - mult_det = [k_ for k_, idx_ in enumerate(idx_nights) - if np.sum(mask_det[idx_]) >= n_samenight] - if len(mult_det) > 0: - k__ = idx_nights[mult_det[0]][mask_det[idx_nights[mult_det[0]]]] - p0 = lc['time'][k__][1] - lc.meta['t0'] + + if np.sum(mask_det) > n_det: + k__ = np.where(mask_det)[0] + p0 = lc['time'][k__][n_det-1] - lc.meta['t0'] if k__[0] > 0: dt = lc['time'][k__[0]] - lc['time'][k__[0] - 1] else: diff --git a/simsurvey/simultarget.py b/simsurvey/simultarget.py index 0c05ad4..2ec2a19 100644 --- a/simsurvey/simultarget.py +++ b/simsurvey/simultarget.py @@ -154,7 +154,7 @@ def save(self, filename): """ """ prop_save = ["transient_coverage", "event_coverage"] - side_save = ["err_mwebv"] + side_save = ["err_mwebv", "apply_mwebv"] deri_save = ["simul_parameters", "mwebv", "mwebv_sfd98", "lightcurve_parameters", "has_mwebv_sfd98"] @@ -875,7 +875,7 @@ def rate_Ia_basic(self,z): More realistic value for low-z than sncosmo default (comoving volumetric rate at each redshift in units of yr^-1 Mpc^-3.) """ - return 3.e-5 + return 3.e-5 * (1 + z) def rate_Ia_basiclow(self,z): @@ -891,21 +891,21 @@ def rate_Ibc_basic(self,z): [TODO: Add source] (comoving volumetric rate at each redshift in units of yr^-1 Mpc^-3.) """ - return 2.25e-5 + return 2.25e-5 * (1 + z) def rate_IIn_basic(self,z): """ [TODO: Add source] (comoving volumetric rate at each redshift in units of yr^-1 Mpc^-3.) """ - return 7.5e-6 + return 7.5e-6 * (1 + z) def rate_IIP_basic(self,z): """ [TODO: Add source] (comoving volumetric rate at each redshift in units of yr^-1 Mpc^-3.) """ - return 1.2e-4 + return 1.2e-4 * (1 + z) # ========================== # # = *** Rates = # @@ -984,12 +984,6 @@ def set_model(self, model): # ============================ # # = Model definitions = # # ============================ # - # Models to be implemented - # - Ibc snana - # - IIn snana - # - IIP snana - # - SpectralIndexSource - # - ExpandingBlackBodySource def model_Ia_salt2(self): """ From 1254237f9764d53b5a664731da6bb6f93604fd0b Mon Sep 17 00:00:00 2001 From: Ulrich Feindt Date: Tue, 2 Apr 2019 15:51:44 +0200 Subject: [PATCH 3/4] bumped version number --- setup.py | 4 ++-- simsurvey/__init__.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/setup.py b/setup.py index cf990d6..e6f3290 100644 --- a/setup.py +++ b/setup.py @@ -13,8 +13,8 @@ MAINTAINER_EMAIL = 'ulrich.feindt@fysik.su.se' URL = 'https://github.com/ufeindt/simsurvey/' LICENSE = 'BSD (3-clause)' -DOWNLOAD_URL = 'https://github.com/ufeindt/simsurvey/tarball/0.4.5' -VERSION = '0.4.5' +DOWNLOAD_URL = 'https://github.com/ufeindt/simsurvey/tarball/0.5.0' +VERSION = '0.5.0' try: from setuptools import setup, find_packages diff --git a/simsurvey/__init__.py b/simsurvey/__init__.py index d49cab7..40bd435 100644 --- a/simsurvey/__init__.py +++ b/simsurvey/__init__.py @@ -1,7 +1,7 @@ """ Module base on astrobject, sncosmo and astropy to prepare a future survey. """ -__version__ = "0.4.5" +__version__ = "0.5.0" from .simultarget import * from .simulsurvey import * From 9630d93e0c6257118b082b4040ec3708cb921b03 Mon Sep 17 00:00:00 2001 From: Ulrich Feindt Date: Wed, 3 Apr 2019 15:57:37 +0200 Subject: [PATCH 4/4] bug fixes and testing done --- simsurvey/simulsurvey.py | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/simsurvey/simulsurvey.py b/simsurvey/simulsurvey.py index 908c66f..eb74282 100644 --- a/simsurvey/simulsurvey.py +++ b/simsurvey/simulsurvey.py @@ -552,7 +552,7 @@ class SurveyPlan( BaseObject ): DERIVED_PROPERTIES = [] def __init__(self, time=None, ra=None, dec=None, band=None, skynoise=None, - obs_field=None, obs_ccd,= None, zp=None, comment=None, + obs_field=None, obs_ccd=None, zp=None, comment=None, width=7.295, height=7.465, fields=None, empty=False, load_opsim=None, **kwargs): """ @@ -616,7 +616,7 @@ def set_fields(self, ra=None, dec=None, ccds=None, **kwargs): # self._update_field_radec() def add_observation(self, time, band, skynoise, ra=None, dec=None, field=None, - zp=None, comment=None): + ccd=None, zp=None, comment=None): """ """ if ra is None and dec is None and field is None: @@ -815,10 +815,12 @@ def observed_on(self, fields=None, ccds=None, if fields is not None: for k, l in enumerate(fields): - if ccds is None: - mask = (self.pointings['field'] == l) - else: - mask = (self.pointings['field'] == l) & self.pointing['ccd'] == ccds[k]) + mask = (self.pointings['field'] == l) + mask2 = np.isnan(self.pointings['ccd']) + if ccds is not None and np.any(~mask2): + mask3 = (self.pointings['ccd'] == ccds[k]) + mask = mask & (mask2 | mask3) + out['time'].extend(self.pointings['time'][mask].quantity.value) out['band'].extend(self.pointings['band'][mask]) out['zp'].extend(self.pointings['zp'][mask]) @@ -1338,15 +1340,11 @@ def get_p_det_last(lc, thr=5., n_det=2): if np.sum(mask_det) > n_det: k__ = np.where(mask_det)[0] - p0 = lc['time'][k__][n_det-1] - lc.meta['t0'] - if k__[0] > 0: - dt = lc['time'][k__[0]] - lc['time'][k__[0] - 1] - else: - dt = 1e12 + p0 = lc['time'][k__][n_det-1] - lc.meta['t0'] + if k__[0] > 0: + dt = lc['time'][k__[0]] - lc['time'][k__[0] - 1] else: - p0 = 1e12 - dt = 1e12 - + dt = 1e12 p1 = lc['time'].max() - lc.meta['t0'] else: p0 = 1e12 @@ -1379,6 +1377,7 @@ def get_lc_max(lc, band): if len(lc_b) > 0: max_flux = np.max(lc_b['flux']) zp = lc_b['zp'][lc_b['flux'] == max_flux] - return -2.5 * np.log10(max_flux) + zp - else: - return 99. + if max_flux > 0: + return -2.5 * np.log10(max_flux) + zp + + return 99.