diff --git a/__init__.py b/__init__deprecated.py similarity index 100% rename from __init__.py rename to __init__deprecated.py diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..743a913a --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,30 @@ +[build-system] +requires = ["setuptools>=61.0.0"] +build-backend = "setuptools.build_meta" + +[project] +name = "sn_simulation" +version = "1.1.0" +requires-python = ">=3.5" +dependencies = ["dustmaps",] +authors = [ + {name = "Philippe Gris", email = "philippe.gris@clermont.in2p3.fr"} +] +maintainers = [ + {name = "Philippe Gris", email = "philippe.gris@clermont.in2p3.fr"} +] +description = "Simulations for SNe Ia" +readme = "README.rst" +license = {file = "LICENSE"} +classifiers = ["Programming Language :: Python",] + +[tool.setuptools.packages.find] +where = [".","sn_simu_input"] +include = ['sn_simulator', 'sn_simu_wrapper', 'sn_simu_input'] + +[tool.setuptools.package-data] +"*" = ["*.txt"] + +[project.urls] +Homepage = "http://github.com/lsstdesc/sn_simulation" +Repository = "https://github.com/lsstdesc/sn_simulation" \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 00000000..ec784615 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +###### Requirements without Version Specifiers ###### +#numpy diff --git a/setup.py b/setup.py deleted file mode 100644 index acd8f0f4..00000000 --- a/setup.py +++ /dev/null @@ -1,25 +0,0 @@ -from setuptools import setup - -# get the version here -pkg_vars = {} - -with open("version.py") as fp: - exec(fp.read(), pkg_vars) - -setup( - name='sn_simulation', - version= pkg_vars['__version__'], - description='Simulations for supernovae', - url='http://github.com/lsstdesc/sn_simulation', - author='Philippe Gris', - author_email='philippe.gris@clermont.in2p3.fr', - license='BSD', - packages=['sn_simulator', 'sn_simu_wrapper'], - python_requires='>=3.5', - zip_safe=False, - install_requires=[ - 'sn_tools>=0.1', - 'sn_stackers>=0.1', - 'dustmaps' - ], -) diff --git a/setup_deprecated.py b/setup_deprecated.py new file mode 100644 index 00000000..d3765ed9 --- /dev/null +++ b/setup_deprecated.py @@ -0,0 +1,34 @@ +from setuptools import setup + +# get the version here +pkg_vars = {} + +with open("version.py") as fp: + exec(fp.read(), pkg_vars) + +setup( + name='sn_simulation', + version=pkg_vars['__version__'], + description='Simulations for supernovae', + url='http://github.com/lsstdesc/sn_simulation', + author='Philippe Gris', + author_email='philippe.gris@clermont.in2p3.fr', + license='BSD', + packages=['sn_simulator', 'sn_simu_wrapper', 'sn_simu_input'], + # All files from folder sn_simu_input + package_data={'sn_simu_input': ['*.txt']}, + python_requires='>=3.5', + zip_safe=False, + install_requires=[ + 'sn_tools>=0.1', + 'dustmaps', + 'sn_telmodel>=0.1' + ], +) + +# pyproject.toml +dependencies = [ + "getObsAtmo@git+https://github.com/lsstdesc/getObsAtmo#egg=getObsAtmo", + "rubin_sim@git+https://github.com/lsst/rubin_sim.git#egg=rubin_sim", + "dustmaps", +] diff --git a/sn_simu_input/__init__.py b/sn_simu_input/__init__.py new file mode 100644 index 00000000..58f3ace6 --- /dev/null +++ b/sn_simu_input/__init__.py @@ -0,0 +1 @@ +from .version import __version__ diff --git a/sn_simu_input/config_simulation.txt b/sn_simu_input/config_simulation.txt new file mode 100644 index 00000000..b17245ba --- /dev/null +++ b/sn_simu_input/config_simulation.txt @@ -0,0 +1,171 @@ +## production Id +ProductionIDSimu prodid str #Production Id + +## SN parameters +SN Id 100 int #SN Id +SN type SN_Ia str #SN type +# parameter dist +SN modelPar nameFile x1_color_G10.csv str #SN dist filename +SN modelPar dirFile reference_files str #SN dist model ref dir +SN modelPar x1sigma 0 int #shift for x1 distrib parameter +SN modelPar colorsigma 0 int #shift for color distrib parameter +SN modelPar mbsigma 0 int #shift for mb SN parameter +SN modelPar mbsigmafile sigma_mb_from_simu_Ny_40.hdf5 str #file for mbsigma application +#x1 +SN x1 type unique str #SN x1 type +SN x1 min -2.0 float #SN x1 min +SN x1 max 0.2 float #SN x1 max +SN x1 step 0.1 float #SN x1 step +# color +SN color type unique str #SN color type +SN color min 0.2 float #SN color min +SN color max 1.0 float #SN color max +SN color step 0.1 float #SN color step +# redshift +SN z type unique str #SN z type +SN z min 0.01 float #SN z min +SN z max 1.0 float #SN z max +SN z step 0.01 float #SN z step +SN z rate combined str #SN z rate +SN z weight sn_rate str #SN z weight(sn_rate/flat) +SN z maxsimu 1.1 float # SN max z for rate +SN z minsimu 0.01 float # SN max z for rate +SN z sigmaz 1.e-5 float # redshift error to be propagated to zsim +# daymax +SN daymax type unique str #SN daymax type +SN daymax step 1. float #SN daymax step +SN daymax restricted 0 int #to restrict daymax values +#rf cuts +SN minRFphase -20. float #SN min rf phase +SN maxRFphase 60. float #SN max rf phase +SN minRFphaseQual -10. float #SN min rf phase qual +SN maxRFphaseQual 35. float #SN max rf phase qual +# other parameters +SN absmag -19.0906 float #SN absmag +SN band bessellB str #SN band +SN magsys vega str #SN magsys +SN differentialFlux 0 int #SN diff flux +SN salt2Dir SALT2_Files str #SN SALT2 dir +SN blueCutoffu 380. float #SN blue cutoff u-band +SN redCutoffu 800. float #SN red cutoff u-band +SN blueCutoffg 380. float #SN blue cutoff g-band +SN redCutoffg 800. float #SN red cutoff g-band +SN blueCutoffr 380. float #SN blue cutoff r-band +SN redCutoffr 800. float #SN red cutoff r-band +SN blueCutoffi 360. float #SN blue cutoff i-band +SN redCutoffi 800. float #SN red cutoff i-band +SN blueCutoffz 380. float #SN blue cutoff z-band +SN redCutoffz 800. float #SN red cutoff z-band +SN blueCutoffy 380. float #SN blue cutoff y-band +SN redCutoffy 800. float #SN red cutoff y-band +SN ebvofMW -1. float #SN E(B-V) +SN NSNfactor 1 int #NSN*factor for simulation +SN NSNabsolute -1 int #absolute number of SN to produce +SN sigmaInt 0.0 float #SN intrinsic dispersion +SN nspectra 0 int #Number of spectra to generate +SN smearFlux 1 int #LC flux smearing +SN alpha 0.13 float #nuisance parameter (alpha) +SN beta 3.1 float #nuisance parameter (beta) +SN x0 griddata 0 int #to extract x0 from griddata +##Simulation parameters(from file) +SN simuFile None str # simulation parameter file + +## Cosmology parameters +Cosmology Model w0waCDM str #cosmology model +Cosmology Om 0.30 float #cosmology Omegam +Cosmology Ol 0.70 float #cosmology Omegall +Cosmology H0 70.0 float #cosmology H0 +Cosmology w0 -1.0 float #cosmology w0 +Cosmology wa 0.0 float #cosmology wa + +#Light curve +LC coadd 1 int #to coadd LC fluxes + +## Instrument +InstrumentSimu name LSST str #instrument name +InstrumentSimu telescope_dir throughputs str #main telescope dir +InstrumentSimu telescope_tag 1.9 str #throughputs tag version +InstrumentSimu throughputDir baseline str #instrument throughput dir +InstrumentSimu atmosDir atmos str #instrument atmos dir +InstrumentSimu atmosType const str #instrument atmos/const/dep +#airmass +InstrumentSimu airmass 1.2 float #instrument airmass +InstrumentSimu sigma_airmass 0.0 float #instrument sigma_airmass +InstrumentSimu round_airmass 1 int #instrument airmass rounding +InstrumentSimu min_airmass 1.0 float #airmass min after smearing +InstrumentSimu max_airmass 2.9 float #airmass max after smearing +#aerosol +InstrumentSimu aerosol 0.05 float #instrument aerosol +InstrumentSimu sigma_aerosol 0.0 float #instrument sigma_aerosol +InstrumentSimu round_aerosol 4 int #instrument aerosol rounding +InstrumentSimu min_aerosol 0. float #aerosol min after smearing +InstrumentSimu max_aerosol 0.5 float #aerosol max after smearing +#pwv +InstrumentSimu pwv 5.0 float #instrument precipit. water vapor +InstrumentSimu sigma_pwv 0.0 float #instrument sigma precipit. water vapor +InstrumentSimu round_pwv 3 int #instrument precipit. water vapor rounding +InstrumentSimu min_pwv 0.0 float #pwv min after smearing +InstrumentSimu max_pwv 25. float #pwv max after smearing +#ozone +InstrumentSimu ozone 300. float #ozone +InstrumentSimu sigma_ozone 0. float #sigma_ozone +InstrumentSimu round_ozone 3 int #ozone rounding +InstrumentSimu min_ozone 0.0 float #ozone min after smearing +InstrumentSimu max_ozone 600. float #ozone max after smearing +#number of trial to estimate zp +InstrumentSimu ntrial_zp 1 int #ntrial to estimate zp vs airmass + + + + + +##Observations +Observations filename fullDbName str #observation file name +Observations fieldtype WFD str #observations field type +Observations fieldname all str #observations field name (DD only) +Observations coadd 0 int #observations coaddition per night +Observations season -1 str #observations seasons + +## Simulator +Simulator name sn_simulator.sn_cosmo str #simulator name +Simulator model salt3 str #simulator model +Simulator version 2.0 float #simulator version +Simulator errorModel 0 int #simulator error model + +## Reference_Files +ReferenceFiles TemplateDir Template_LC str #dir for templates ref files +ReferenceFiles GammaDir reference_files str #dir for gamma ref files +ReferenceFiles GammaFile gamma.hdf5 str #gamma ref file name +ReferenceFiles DustCorrDir Template_Dust str #dir for template dust files +ReferenceFiles fluxpixelDir reference_files str #dir for template dust files +# ReferenceFiles zpDir reference_files str #dir for zp_airmass ref file +# ReferenceFiles zpFile zp_airmass_v1.9.npy str #zp_airmass ref file + +## Host +Host 0 int #Host + +# Display LC +Display LC display 0 int #display LC +Display LC time 5. float #display LC persistency time + +##Output: +OutputSimu directory Output_Simu str #Output directory +OutputSimu save 1 int #output save file in sn_simulator +OutputSimu savefromwrapper 0 int #output save file in simuwrapper +OutputSimu throwempty 1 int #do not save empty LC +OutputSimu throwafterdump 1 int #remove LC after dumping on file +OutputSimu clean 1 int #to clean output dir before processing + +#Multiprocessing +MultiprocessingSimu nproc 1 int #multiprocessing number of procs + +# Pixelisation +Pixelisation nside 64 int #pixelisation nside Healpix + +#Web path +WebPathSimu https://me.lsst.eu/gris/DESC_SN_pipeline str #web path for reference files + +#include saturation effects +saturation effect 0 int #to include saturation effects +saturation psf single_gauss str #PSF for saturation effects +saturation ccdfullwell 90.e3 float #ccd full well value diff --git a/sn_simu_input/sncosmo_builtins.txt b/sn_simu_input/sncosmo_builtins.txt new file mode 100644 index 00000000..c6e22595 --- /dev/null +++ b/sn_simu_input/sncosmo_builtins.txt @@ -0,0 +1,121 @@ +name version type subtype subclass reference website notes +nugent-sn1a 1.2 SN Ia TimeSeriesSource [N02] a +nugent-sn91t 1.1 SN Ia TimeSeriesSource [S04] a +nugent-sn91bg 1.1 SN Ia TimeSeriesSource [N02] a +nugent-sn1bc 1.1 SN Ib/c TimeSeriesSource [L05] a +nugent-hyper 1.2 SN Ib/c TimeSeriesSource [L05] a +nugent-sn2p 1.2 SN IIP TimeSeriesSource [G99] a +nugent-sn2l 1.2 SN IIL TimeSeriesSource [G99] a +nugent-sn2n 2.1 SN IIn TimeSeriesSource [G99] a +s11-2004hx 1.0 SN IIL/P TimeSeriesSource [S11] b [1] +s11-2005lc 1.0 SN IIP TimeSeriesSource [S11] b [1] +s11-2005hl 1.0 SN Ib TimeSeriesSource [S11] b [1] +s11-2005hm 1.0 SN Ib TimeSeriesSource [S11] b [1] +s11-2005gi 1.0 SN IIP TimeSeriesSource [S11] b [1] +s11-2006fo 1.0 SN Ic TimeSeriesSource [S11] b [1] +s11-2006jo 1.0 SN Ib TimeSeriesSource [S11] b [1] +s11-2006jl 1.0 SN IIP TimeSeriesSource [S11] b [1] +hsiao 1.0 SN Ia TimeSeriesSource [H07] c [2] +hsiao 2.0 SN Ia TimeSeriesSource [H07] c [2] +hsiao 3.0 SN Ia TimeSeriesSource [H07] c [2] +hsiao-subsampled 3.0 SN Ia TimeSeriesSource [H07] c [2] +salt2 2.0 SN Ia SALT2Source [G10] d +salt2 2.4 SN Ia SALT2Source [B14b] d +salt2-extended 1.0 SN Ia SALT2Source b [3] +salt2-extended 2.0 SN Ia SALT2Source +salt2-extended-h17 1.0 SN Ia SALT2Source [H17] e [4] +snf-2011fe 1.0 SN Ia TimeSeriesSource [P13] f +snana-2004fe 1.0 SN Ic TimeSeriesSource g [5] +snana-2004gq 1.0 SN Ic TimeSeriesSource g [5] +snana-sdss004012 1.0 SN Ic TimeSeriesSource g [5] +snana-2006fo 1.0 SN Ic TimeSeriesSource g [5] +snana-sdss014475 1.0 SN Ic TimeSeriesSource g [5] +snana-2006lc 1.0 SN Ic TimeSeriesSource g [5] +snana-2007ms 1.0 SN II-pec TimeSeriesSource g [5] +snana-04d1la 1.0 SN Ic TimeSeriesSource g [5] +snana-04d4jv 1.0 SN Ic TimeSeriesSource g [5] +snana-2004gv 1.0 SN Ib TimeSeriesSource g [5] +snana-2006ep 1.0 SN Ib TimeSeriesSource g [5] +snana-2007y 1.0 SN Ib TimeSeriesSource g [5] +snana-2004ib 1.0 SN Ib TimeSeriesSource g [5] +snana-2005hm 1.0 SN Ib TimeSeriesSource g [5] +snana-2006jo 1.0 SN Ib TimeSeriesSource g [5] +snana-2007nc 1.0 SN Ib TimeSeriesSource g [5] +snana-2004hx 1.0 SN IIP TimeSeriesSource g [5] +snana-2005gi 1.0 SN IIP TimeSeriesSource g [5] +snana-2006gq 1.0 SN IIP TimeSeriesSource g [5] +snana-2006kn 1.0 SN IIP TimeSeriesSource g [5] +snana-2006jl 1.0 SN IIP TimeSeriesSource g [5] +snana-2006iw 1.0 SN IIP TimeSeriesSource g [5] +snana-2006kv 1.0 SN IIP TimeSeriesSource g [5] +snana-2006ns 1.0 SN IIP TimeSeriesSource g [5] +snana-2007iz 1.0 SN IIP TimeSeriesSource g [5] +snana-2007nr 1.0 SN IIP TimeSeriesSource g [5] +snana-2007kw 1.0 SN IIP TimeSeriesSource g [5] +snana-2007ky 1.0 SN IIP TimeSeriesSource g [5] +snana-2007lj 1.0 SN IIP TimeSeriesSource g [5] +snana-2007lb 1.0 SN IIP TimeSeriesSource g [5] +snana-2007ll 1.0 SN IIP TimeSeriesSource g [5] +snana-2007nw 1.0 SN IIP TimeSeriesSource g [5] +snana-2007ld 1.0 SN IIP TimeSeriesSource g [5] +snana-2007md 1.0 SN IIP TimeSeriesSource g [5] +snana-2007lz 1.0 SN IIP TimeSeriesSource g [5] +snana-2007lx 1.0 SN IIP TimeSeriesSource g [5] +snana-2007og 1.0 SN IIP TimeSeriesSource g [5] +snana-2007ny 1.0 SN IIP TimeSeriesSource g [5] +snana-2007nv 1.0 SN IIP TimeSeriesSource g [5] +snana-2007pg 1.0 SN IIP TimeSeriesSource g [5] +snana-2006ez 1.0 SN IIn TimeSeriesSource g [5] +snana-2006ix 1.0 SN IIn TimeSeriesSource g [5] +snana-2004fe 2.0 SN Ic TimeSeriesSource +snana-2004gq 2.0 SN Ic TimeSeriesSource +snana-sdss004012 2.0 SN Ic TimeSeriesSource +snana-2006fo 2.0 SN Ic TimeSeriesSource +snana-sdss014475 2.0 SN Ic TimeSeriesSource +snana-2006lc 2.0 SN Ic TimeSeriesSource +snana-2007ms 2.0 SN II-pec TimeSeriesSource +snana-04d1la 2.0 SN Ic TimeSeriesSource +snana-04d4jv 2.0 SN Ic TimeSeriesSource +snana-2004gv 2.0 SN Ib TimeSeriesSource +snana-2006ep 2.0 SN Ib TimeSeriesSource +snana-2007y 2.0 SN Ib TimeSeriesSource +snana-2004ib 2.0 SN Ib TimeSeriesSource +snana-2005hm 2.0 SN Ib TimeSeriesSource +snana-2006jo 2.0 SN Ib TimeSeriesSource +snana-2007nc 2.0 SN Ib TimeSeriesSource +snana-2004hx 2.0 SN IIP TimeSeriesSource +snana-2005gi 2.0 SN IIP TimeSeriesSource +snana-2006gq 2.0 SN IIP TimeSeriesSource +snana-2006kn 2.0 SN IIP TimeSeriesSource +snana-2006jl 2.0 SN IIP TimeSeriesSource +snana-2006iw 2.0 SN IIP TimeSeriesSource +snana-2006kv 2.0 SN IIP TimeSeriesSource +snana-2006ns 2.0 SN IIP TimeSeriesSource +snana-2007iz 2.0 SN IIP TimeSeriesSource +snana-2007nr 2.0 SN IIP TimeSeriesSource +snana-2007kw 2.0 SN IIP TimeSeriesSource +snana-2007ky 2.0 SN IIP TimeSeriesSource +snana-2007lj 2.0 SN IIP TimeSeriesSource +snana-2007lb 2.0 SN IIP TimeSeriesSource +snana-2007ll 2.0 SN IIP TimeSeriesSource +snana-2007nw 2.0 SN IIP TimeSeriesSource +snana-2007ld 2.0 SN IIP TimeSeriesSource +snana-2007md 2.0 SN IIP TimeSeriesSource +snana-2007lz 2.0 SN IIP TimeSeriesSource +snana-2007lx 2.0 SN IIP TimeSeriesSource +snana-2007og 2.0 SN IIP TimeSeriesSource +snana-2007ny 2.0 SN IIP TimeSeriesSource +snana-2007nv 2.0 SN IIP TimeSeriesSource +snana-2007pg 2.0 SN IIP TimeSeriesSource +whalen-z15b 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z15d 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z15g 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z25b 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z25d 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z25g 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z40b 1.0 Pop III TimeSeriesSource [Whalen13] [6] +whalen-z40g 1.0 Pop III TimeSeriesSource [Whalen13] [6] +mlcs2k2 1.0 SN Ia MLCS2k2Source [Jha07] [7] +snemo2 1.0 SN Ia SNEMOSource [Saunders18] h +snemo7 1.0 SN Ia SNEMOSource [Saunders18] h +snemo15 1.0 SN Ia SNEMOSource [Saunders18] h diff --git a/sn_simu_input/version.py b/sn_simu_input/version.py new file mode 100644 index 00000000..aed115e2 --- /dev/null +++ b/sn_simu_input/version.py @@ -0,0 +1 @@ +__version__='v1.0.0' diff --git a/sn_simu_wrapper/config_simulation.py b/sn_simu_wrapper/config_simulation.py new file mode 100644 index 00000000..4056c619 --- /dev/null +++ b/sn_simu_wrapper/config_simulation.py @@ -0,0 +1,93 @@ +from sn_tools.sn_io import recursive_merge + +class ConfigSimulation: + """ + class to load a set of parameters (txt file) + and make a dict out of them + + Parameters + -------------- + type_sn: str + SN type + model_sn: str + model for simulation + config_file: str + configuration file + + """ + + + def __init__(self,type_sn,model_sn,config_file): + + # first: load all possible parameters + cdict = self.configDict(config_file) + + # second: keep keys only valids for (type_sn, model_sn) + + self.conf_dict = self.select(cdict,type_sn,model_sn) + + def configDict(self,config_file): + """ + Method to load a txt configuration file + and makes a dict out of it + + Parameters + -------------- + config_file: str + config file (txt) to transform to a dict + + Returns + ---------- + the dict + + """ + + ffile = open(config_file, 'r') + line = ffile.read().splitlines() + ffile.close() + + params = {} + for i,ll in enumerate(line): + if ll!='' and ll[0]!='#': + lla = ll.split('#') + lspl = lla[0].split(' ') + lspl = ' '.join(lspl).split() #remove empty char here + n = len(lspl) + mystr = '' + myclose = '' + for keya in [lspl[i] for i in range(n-2)]: + mystr += '{\''+keya+ '\':' + myclose +=' }' + + if lspl[n-1]!= 'str': + dd = '{} {} {}'.format(mystr,eval('{}({})'.format(lspl[n-1],lspl[n-2])),myclose) + else: + dd = '{} \'{}\' {}'.format(mystr,lspl[n-2],myclose) + + thedict = eval(dd) + params = recursive_merge(params, thedict) + + return params + + def select(self, thedict, type_sn, model_sn): + """ + Method to select thedict keys according to type_sn and model_sn + + Parameters + --------------- + thedict: dict + original dict + type_sn: str + SN type + model_sn: str + SN model for simulation + """ + + res = dict(thedict) + if 'Ia' not in type_sn: + del res['SN']['x1'] + del res['SN']['color'] + del res['SN']['modelPar'] + + return res + diff --git a/sn_simu_wrapper/sn_object.py b/sn_simu_wrapper/sn_object.py index 4f7d8988..41c2886d 100644 --- a/sn_simu_wrapper/sn_object.py +++ b/sn_simu_wrapper/sn_object.py @@ -1,17 +1,22 @@ import numpy as np -import astropy.units as u -from astropy.table import Table -from collections import OrderedDict as odict +# import astropy.units as u +# from astropy.table import Table +# from collections import OrderedDict as odict class SN_Object: - def __init__(self, name, sn_parameters, gen_parameters, cosmology, - telescope, snid, area, x0_grid, salt2Dir='SALT2_Files', + def __init__(self, name, sn_parameters, simulator_parameters, + gen_parameters, cosmology, + snid, area, x0_grid=None, + salt2Dir='SALT2_Files', mjdCol='mjd', RACol='pixRa', DecCol='pixDec', - filterCol='band', exptimeCol='exptime', nexpCol='numExposures', - m5Col='fiveSigmaDepth', seasonCol='season', + filterCol='band', + exptimeCol='exptime', + nexpCol='numExposures', + nightCol='night', m5Col='fiveSigmaDepth', seasonCol='season', seeingEffCol='seeingFwhmEff', seeingGeomCol='seeingFwhmGeom', - airmassCol='airmass', skyCol='sky', moonCol='moonPhase'): + airmassCol='airmass', skyCol='sky', moonCol='moonPhase', + psf_flux='', frac_flux_seeing=None, ccd_full_well=-1.): """ class SN object handles sn name, parameters, cosmology, snid, telescope... @@ -28,7 +33,6 @@ def __init__(self, name, sn_parameters, gen_parameters, cosmology, simulation parameters cosmology: dict cosmological parameters used for simulation - telescope: sn_tools.Telescope snid: int supernova identifier area: float @@ -66,17 +70,18 @@ def __init__(self, name, sn_parameters, gen_parameters, cosmology, """ self._name = name self._sn_parameters = sn_parameters + self._simulator_parameters = simulator_parameters self._gen_parameters = gen_parameters self._cosmology = cosmology - self._telescope = telescope self._SNID = snid - self.mjdCol = mjdCol self.RACol = RACol self.DecCol = DecCol + self.filterCol = filterCol self.exptimeCol = exptimeCol self.nexpCol = nexpCol + self.nightCol = nightCol self.m5Col = m5Col self.seasonCol = seasonCol self.seeingEffCol = seeingEffCol @@ -89,35 +94,39 @@ def __init__(self, name, sn_parameters, gen_parameters, cosmology, self.salt2Dir = salt2Dir self.x0_grid = x0_grid - @property + self.psf_flux = psf_flux + self.frac_flux_seeing = frac_flux_seeing + self.ccd_full_well = ccd_full_well + + @ property def name(self): return self._name - @property + @ property def sn_parameters(self): """SN parameters """ return self._sn_parameters - @property + @ property + def simulator_parameters(self): + """SN parameters + """ + return self._simulator_parameters + + @ property def gen_parameters(self): """ Simulation parameters """ return self._gen_parameters - @property + @ property def cosmology(self): """ Cosmology """ return self._cosmology - @property - def telescope(self): - """ Telescope - """ - return self._telescope - - @property + @ property def SNID(self): """ SN identifier """ @@ -125,7 +134,10 @@ def SNID(self): def cutoff(self, obs, T0, z, min_rf_phase, max_rf_phase, - blue_cutoff=350., red_cutoff=800.): + blue_cutoffs=dict( + zip('ugrizy', [380., 380., 380., 360., 380., 380.])), + red_cutoffs=dict(zip('ugrizy', + [700., 700., 700., 700., 700., 700.]))): """ select observations depending on phases Parameters @@ -146,18 +158,117 @@ def cutoff(self, obs, T0, z, array of obs passing the selection """ - mean_restframe_wavelength = np.asarray( - [self.telescope.mean_wavelength[obser[self.filterCol][-1]] / - (1. + z) for obser in obs]) + self.blue_cutoffs = blue_cutoffs + self.red_cutoffs = red_cutoffs + + filters = np.array(obs[self.filterCol].tolist()) + # airmass = np.array(obs['airmass'].tolist()) + + # filt_air = list(zip(filters,airmass)) + filters = filters.reshape((len(filters), 1)) + + blue_values = np.apply_along_axis(self.blues, 1, filters) + red_values = np.apply_along_axis(self.reds, 1, filters) + + """ + mean_restframe_wavelength = \ + np.apply_along_axis(self.mean_wave, 1, filters,airmass)/(1.+z) + mean_restframe_wavelength= mean_restframe_wavelength.reshape((len(filters),1)) + """ + # print('jjj',self.mean_wavelength_airmass['g'](1.2)) + + # filt_air = np.array([*map(self.mean_wavelength_airmass.get, filters)]) + # res = list(map(lambda x:self.mean_wavelength_airmass[x[0]](x[1]),filt_air)) + + mean_wavelength = obs['mean_wave']/(1.+z) p = (obs[self.mjdCol]-T0)/(1.+z) idx = (p >= 1.000000001*min_rf_phase) & (p <= 1.00001*max_rf_phase) + """ idx &= (mean_restframe_wavelength > blue_cutoff) idx &= (mean_restframe_wavelength < red_cutoff) - return obs[idx] + """ + idx &= (mean_wavelength - blue_values >= 0.) + idx &= (mean_wavelength - red_values <= 0.) + + selobs = obs[idx] + + return selobs + + def blues(self, band): + """ + Method to return blue_cutoff value + + Parameters + ---------- + band: str + the band to process + + Returns + ------- + the blue cutoff value + + """ + return self.blue_cutoffs[band[0]] + + def reds(self, band): + """ + Method to return the red_cutoff value + + Parameters + ---------- + band: str + the band to process + + Returns + ------- + the red cutoff value + + """ + + return self.red_cutoffs[band[0]] + + def mean_wave(self, band): + """ + Method to return the mean_restframe_wavelength + + Parameters + ---------- + band : str + the band to process + + Returns + ------- + float + the mean restframe wavelength corresponding to band + + """ - def plotLC(self, table, time_display): + # return self.mean_wavelength[band[0]] + return self.mean_wavelength_airmass[band[0]] + + def mean_wave_airmass(self, vv): + """ + Method to return the mean_restframe_wavelength + + Parameters + ---------- + band : str + the band to process + + Returns + ------- + float + the mean restframe wavelength corresponding to band + + """ + + print('hello', vv) + return self.filt_air(vv) + + @ staticmethod + def plotLC(table, time_display, airmass=1.2): """ Light curve plot using sncosmo methods Parameters @@ -170,7 +281,7 @@ def plotLC(self, table, time_display): import pylab as plt import sncosmo - prefix = 'LSST::' + # prefix = 'LSST::' """ _photdata_aliases = odict([ ('time', set(['time', 'date', 'jd', 'mjd', 'mjdobs', 'mjd_obs'])), @@ -182,42 +293,33 @@ def plotLC(self, table, time_display): ('zpsys', set(['zpsys', 'zpmagsys', 'magsys'])) ]) """ - for band in 'grizy': - name_filter = prefix+band - if self.telescope.airmass > 0: - bandpass = sncosmo.Bandpass( - self.telescope.atmosphere[band].wavelen, - self.telescope.atmosphere[band].sb, - name=name_filter, - wave_unit=u.nm) - else: - bandpass = sncosmo.Bandpass( - self.telescope.system[band].wavelen, - self.telescope.system[band].sb, - name=name_filter, - wave_unit=u.nm) - # print('registering',name_filter) - sncosmo.registry.register(bandpass, force=True) - z = table.meta['z'] - x1 = table.meta['x1'] - color = table.meta['color'] + if 'x1' in table.meta.keys(): + x1 = table.meta['x1'] + color = table.meta['color'] + x0 = table.meta['x0'] + else: + x1 = 0. + color = 0. + x0 = 0. daymax = table.meta['daymax'] model = sncosmo.Model('salt2') model.set(z=z, c=color, t0=daymax, - # x0=self.X0, + # x0=x0, x1=x1) """ - print('tests',isinstance(table, np.ndarray),isinstance(table,Table),isinstance(table,dict)) + print('tests',isinstance(table, np.ndarray), + isinstance(table,Table),isinstance(table,dict)) array_tab = np.asarray(table) print(array_tab.dtype) colnames = array_tab.dtype.names # Create mapping from lowercased column names to originals - lower_to_orig = dict([(colname.lower(), colname) for colname in colnames]) - + lower_to_orig = dict([(colname.lower(), colname) + for colname in colnames]) + # Set of lowercase column names lower_colnames = set(lower_to_orig.keys()) orig_colnames_to_use = [] @@ -228,12 +330,19 @@ def plotLC(self, table, time_display): '(case independent)'.format(', '.join(aliases))) orig_colnames_to_use.append(lower_to_orig[i.pop()]) - + new_data = table[orig_colnames_to_use].copy() print('bbbb',orig_colnames_to_use,_photdata_aliases.keys(),new_data.dtype.names) new_data.dtype.names = _photdata_aliases.keys() """ - sncosmo.plot_lc(data=table, model=model) + # display only 1 sigma LC points + table = table[table['flux']/table['fluxerr'] >= 1.] + """ + if 'x1' in table.meta.keys(): + sncosmo.plot_lc(data=table,model=model) + else: + """ + sncosmo.plot_lc(data=table) plt.draw() plt.pause(time_display) diff --git a/sn_simu_wrapper/sn_simu.py b/sn_simu_wrapper/sn_simu.py index 85f7c560..4ece5c36 100644 --- a/sn_simu_wrapper/sn_simu.py +++ b/sn_simu_wrapper/sn_simu.py @@ -1,73 +1,51 @@ import numpy as np -from lsst.sims.maf.metrics import BaseMetric -from sn_stackers.coadd_stacker import CoaddStacker import healpy as hp import os import time -import multiprocessing -from astropy.table import Table -import h5py +# import multiprocessing +import astropy +from astropy.table import Table, vstack, unique from astropy.cosmology import w0waCDM from importlib import import_module -from sn_tools.sn_telescope import Telescope from sn_simu_wrapper.sn_object import SN_Object -from sn_tools.sn_utils import SimuParameters +from sn_tools.sn_utils import SimuParameters, multiproc from sn_tools.sn_obs import season as seasoncalc -from sn_tools.sn_utils import GetReference, LoadGamma, LoadDust -from scipy.interpolate import interp1d -from sn_tools.sn_io import check_get_dir - -class SNSimulation(BaseMetric): - """LC simulation wrapper class - - Parameters - --------------- - - mjdCol: str, opt - mjd col name in observations (default: 'observationStartMJD') - RACol: str, opt - RA col name in observations (default: 'fieldRA') - DecCol:str, opt - Dec col name in observations (default: 'fieldDec') - filterCol: str, opt - filter col name in observations (default: filter') - m5Col: str, opt - 5-sigma depth col name in observations (default: 'fiveSigmaDepth') - exptimeCol: str, opt - exposure time col name in observations (default: 'visitExposureTime) - nightCol: str, opt - night col name in observations (default: 'night') - obsidCol: str, opt - observation id col name in observations (default: 'observationId') - nexpCol: str, opt - number of exposures col name in observations (default: 'numExposures') - visitimeCol: str, opt - visit time col name in observations (default: 'visiTime') - seeingEffCol: str, opt - seeing eff col name in observations (default: 'seeingFwhmEff') - seeingGeomCol: str, opt - seeing geom col name in observations (default: 'seeingFwhmGeom') - coadd: bool, opt - coaddition of obs (per band and per night) if set to True (default: True) - config: dict - configuration dict for simulation (SN parameters, cosmology, telescope, ...) - ex: {'ProductionID': 'DD_baseline2018a_Cosmo', 'SN parameters': {'Id': 100, 'x1_color': {'type': 'fixed', 'min': [-2.0, 0.2], 'max': [0.2, 0.2], 'rate': 'JLA'}, 'z': {'type': 'uniform', 'min': 0.01, 'max': 0.9, 'step': 0.05, 'rate': 'Perrett'}, 'daymax': {'type': 'unique', 'step': 1}, 'min_rf_phase': -20.0, 'max_rf_phase': 60.0, 'absmag': -19.0906, 'band': 'bessellB', 'magsys': 'vega', 'differential_flux': False}, 'Cosmology': {'Model': 'w0waCDM', 'Omega_m': 0.3, 'Omega_l': 0.7, 'H0': 72.0, 'w0': -1.0, 'wa': 0.0}, 'Instrument': { - 'name': 'LSST', 'throughput_dir': 'LSST_THROUGHPUTS_BASELINE', 'atmos_dir': 'THROUGHPUTS_DIR', 'airmass': 1.2, 'atmos': True, 'aerosol': False}, 'Observations': {'filename': '/home/philippe/LSST/DB_Files/kraken_2026.db', 'fieldtype': 'DD', 'coadd': True, 'season': 1}, 'Simulator': {'name': 'sn_simulator.sn_cosmo', 'model': 'salt2-extended', 'version': 1.0, 'Reference File': 'LC_Test_today.hdf5'}, 'Host Parameters': 'None', 'Display_LC': {'display': True, 'time': 1}, 'Output': {'directory': 'Output_Simu', 'save': True}, 'Multiprocessing': {'nproc': 1}, 'Metric': 'sn_mafsim.sn_maf_simulation', 'Pixelisation': {'nside': 64}} - x0_norm: array of float - grid ox (x1,color,x0) values +from sn_tools.sn_calcFast import GetReference, LoadDust +from sn_tools.sn_stacker import CoaddStacker +from sn_telmodel.sn_throughputs import load_throughputs_from_config +import numpy.lib.recfunctions as rf +import pandas as pd +from sn_tools.sn_utils import register_bands_sncosmo +# import sncosmo as sncosmo_emul +from sn_tools.sn_obs import load_season +# import tracemalloc + + +class SNSimu_Params: + def __init__(self, mjdCol, + RACol, DecCol, + filterCol, m5Col, + exptimeCol, + nightCol, obsidCol, + nexpCol, + vistimeCol, seeingEffCol, + airmassCol, + skyCol, moonCol, + seeingGeomCol, config, dust_map, zp_airmass=None): + """ + Class to load simulation parameters - """ + Parameters + ---------- + config : dict + dict of parameters - def __init__(self, metricName='SNSimulation', - mjdCol='observationStartMJD', RACol='fieldRA', DecCol='fieldDec', - filterCol='filter', m5Col='fiveSigmaDepth', exptimeCol='visitExposureTime', - nightCol='night', obsidCol='observationId', nexpCol='numExposures', - vistimeCol='visitTime', seeingEffCol='seeingFwhmEff', - airmassCol='airmass', - skyCol='sky', moonCol='moonPhase', - seeingGeomCol='seeingFwhmGeom', - uniqueBlocks=False, config=None, x0_norm=None,**kwargs): + Returns + ------- + None. + """ + # data columns self.mjdCol = mjdCol self.m5Col = m5Col self.filterCol = filterCol @@ -84,28 +62,14 @@ def __init__(self, metricName='SNSimulation', self.airmassCol = airmassCol self.skyCol = skyCol self.moonCol = moonCol + self.obs_coadd = config['Observations']['coadd'] + self.lc_coadd = config['LC']['coadd'] - cols = [self.RACol, self.DecCol, self.nightCol, self.m5Col, self.filterCol, self.mjdCol, self.obsidCol, - self.nexpCol, self.vistimeCol, self.exptimeCol, self.seeingEffCol, self.seeingGeomCol, self.nightCol, - self.airmassCol, self.moonCol] - - # self.airmassCol, self.skyCol, self.moonCol] - self.stacker = None - - coadd = config['Observations']['coadd'] - - if coadd: - # cols += ['sn_coadd'] - self.stacker = CoaddStacker(mjdCol=self.mjdCol, - RACol=self.RACol, DecCol=self.DecCol, - m5Col=self.m5Col, nightCol=self.nightCol, - filterCol=self.filterCol, numExposuresCol=self.nexpCol, - visitTimeCol=self.vistimeCol, visitExposureTimeCol='visitExposureTime') - super(SNSimulation, self).__init__( - col=cols, metricName=metricName, **kwargs) + # load stacker + self.stacker = self.load_stacker(config['Observations']['coadd']) # bands considered - self.filterNames = 'ugrizy' + self.filterNames = 'grizy' # grab config file self.config = config @@ -115,59 +79,85 @@ def __init__(self, metricName='SNSimulation', self.area = hp.nside2pixarea(self.nside, degrees=True) # prodid - self.prodid = config['ProductionID'] + self.prodid = config['ProductionIDSimu'] # load cosmology - cosmo_par = config['Cosmology'] - self.cosmology = w0waCDM(H0=cosmo_par['H0'], - Om0=cosmo_par['Omega_m'], - Ode0=cosmo_par['Omega_l'], - w0=cosmo_par['w0'], wa=cosmo_par['wa']) - - # load telescope - tel_par = config['Instrument'] - self.telescope = Telescope(name=tel_par['name'], - throughput_dir=tel_par['throughput_dir'], - atmos_dir=tel_par['atmos_dir'], - atmos=tel_par['atmos'], - aerosol=tel_par['aerosol'], - airmass=tel_par['airmass']) + self.cosmology = self.load_cosmology(config['Cosmology']) # sn parameters - self.sn_parameters = config['SN parameters'] - self.gen_par = SimuParameters(self.sn_parameters, cosmo_par, mjdCol=self.mjdCol, area=self.area, - dirFiles=self.sn_parameters['x1_color']['dirFile'], - web_path=config['Web path']) + self.sn_parameters = config['SN'] + + """ + dirFiles = None + if 'modelPar' in self.sn_parameters.keys(): + dirFiles = self.sn_parameters['modelPar']['dirFile'] + """ + self.gen_par = SimuParameters(self.sn_parameters, config['Cosmology'], + mjdCol=self.mjdCol, area=self.area, + web_path=config['WebPathSimu']) + + # simulator parameters + self.simulator_parameters = config['Simulator'] + # simu params from file + """ + self.simuParamsFile = self.simu_params_from_file( + self.sn_parameters['simuFile']) + """ # this is for output - save_status = config['Output']['save'] + save_status = config['OutputSimu']['save'] self.save_status = save_status - self.outdir = config['Output']['directory'] + clean_status = config['OutputSimu']['clean'] + self.clean_status = clean_status + self.outdir = config['OutputSimu']['directory'] + self.throwaway_empty = config['OutputSimu']['throwempty'] + self.throwafterdump = config['OutputSimu']['throwafterdump'] # number of procs to run simu here - self.nprocs = config['Multiprocessing']['nproc'] + self.nprocs = config['MultiprocessingSimu']['nproc'] - # if saving activated, prepare output dirs - """ - if self.save_status: - self.prepareSave(outdir, prodid) - """ # simulator parameter self.simu_config = config['Simulator'] + # reference files + self.reffiles = config['ReferenceFiles'] + + # load zp vs airmass + """ + self.zp_airmass = self.load_zp(config['WebPathSimu'], + self.reffiles['zpDir'], + self.reffiles['zpFile']) + """ + # load flux_pixel if necessary for saturation effects + self.frac_flux_seeing = None + self.ccd_full_well = -1. + self.psf_flux = 'None' + saturation_effect = config['saturation']['effect'] + if saturation_effect: + psf = config['saturation']['psf'] + fName = 'psf_pixel_{}_summary.npy'.format(psf) + dira = config['WebPathSimu'] + dirb = self.reffiles['fluxpixelDir'] + self.frac_flux_seeing = self.load_flux_to_pixel(dira, + dirb, fName) + self.ccd_full_well = config['saturation']['ccdfullwell'] + self.psf_flux = psf # LC display in "real-time" - self.display_lc = config['Display_LC']['display'] - self.time_display = config['Display_LC']['time'] - # fieldtype, season + self.display_lc = config['Display']['LC']['display'] + self.time_display = config['Display']['LC']['time'] + + # fieldtype self.field_type = config['Observations']['fieldtype'] - self.season = config['Observations']['season'] + + # seasons + self.season = load_season(config['Observations']['season']) self.type = 'simulation' # get the x0_norm values to be put on a 2D(x1,color) griddata - self.x0_grid = x0_norm + # self.x0_grid = x0_norm # SALT2DIR self.salt2Dir = self.sn_parameters['salt2Dir'] @@ -177,47 +167,24 @@ def __init__(self, metricName='SNSimulation', # load reference LC if simulator is sn_fast self.reference_lc = None - self.gamma = None - self.mag_to_flux = None + # self.gamma = None + # self.mag_to_flux = None self.dustcorr = None - web_path = config['Web path'] - self.error_model = self.simu_config['error_model'] - - if 'sn_fast' in self.simu_config['name']: - templateDir = self.simu_config['Template Dir'] - gammaDir = self.simu_config['Gamma Dir'] - gammaFile = self.simu_config['Gamma File'] - dustDir = self.simu_config['DustCorr Dir'] - - # x1 and color are unique for this simulator - x1 = self.sn_parameters['x1']['min'] - color = self.sn_parameters['color']['min'] - bluecutoff = self.sn_parameters['blue_cutoff'] - redcutoff = self.sn_parameters['red_cutoff'] - # Loading reference file - lcname = 'LC_{}_{}_{}_{}_ebvofMW_0.0_vstack.hdf5'.format( - x1, color, bluecutoff, redcutoff) - - self.reference_lc = GetReference(templateDir, - lcname, gammaDir, gammaFile, web_path, self.telescope) - - dustFile = 'Dust_{}_{}_{}_{}.hdf5'.format( - x1, color, bluecutoff, redcutoff) - self.dustcorr = LoadDust(dustDir, dustFile, web_path).dustcorr + web_path = config['WebPathSimu'] + self.error_model = self.simu_config['errorModel'] + if 'sn_fast' in self.simu_config['name']: + self.reference_lc, self.dustcorr = self.load_for_snfast(web_path) + """ else: gammas = LoadGamma( - 'grizy', self.simu_config['Gamma Dir'], - self.simu_config['Gamma File'], - web_path, self.telescope) - + 'grizy', self.reffiles['GammaDir'], + self.reffiles['GammaFile'], + web_path) self.gamma = gammas.gamma self.mag_to_flux = gammas.mag_to_flux - - - if self.error_model: - SALT2Dir = 'SALT2.Guy10_UV2IR' - check_get_dir(web_path,SALT2Dir,SALT2Dir) + """ + self.filterNames = ['g', 'r', 'i', 'z', 'y'] self.nprocdict = {} self.simu_out = {} @@ -225,136 +192,1494 @@ def __init__(self, metricName='SNSimulation', self.SNID = {} self.sn_meta = {} - def run(self, obs, slicePoint=None): + # load the instrument(telescope) + self.telescope = load_throughputs_from_config(config['InstrumentSimu']) + + self.config_instr = config['InstrumentSimu'] + self.atmosType = self.config_instr['atmosType'] + + # estimate zp , mean_wave and sigmas vs airmass + if zp_airmass is None: + from sn_telmodel.sn_transtools import zp_from_config + self.zp_atmos = zp_from_config(config['InstrumentSimu']) + else: + self.zp_atmos = zp_airmass + + # load dust_map + self.dust_map = dust_map + + def simu_params_from_file_deprecated(self, simuFile): + """ + Method to grab simu parameters from file + + Parameters + ---------- + simuFile : str + simu file Name. + + Returns + ------- + numpy array + array with simu parameters + + """ + # sn simu parameters from file + df = pd.DataFrame() + simuFile = self.sn_parameters['simuFile'] + + if simuFile != 'None': + df = pd.read_hdf(simuFile) + else: + return df + + # complete df with other simulation parameters + ccols = ['healpixID', 'season', 'z', 'daymax', 'x1', 'color', + 'epsilon_x0', 'epsilon_x1', + 'epsilon_color', 'epsilon_daymax', 'SNID'] + if 'weight' in df.columns: + ccols += ['weight'] + ccolsb = ['minRFphase', 'maxRFphase', + 'minRFphaseQual', 'maxRFphaseQual'] + + for vv in ccolsb: + df[vv] = self.sn_parameters[vv] + + return df[ccols+ccolsb].to_records(index=False) + + def load_for_snfast(self, web_path): + """ + Method to load reference files for sn_fast + + Parameters + ---------- + web_path : str + web dir where files are located. + + Returns + ------- + reference_lc : astropy table + lc reference files + dustcorr : astropy table + dust correction data. + + """ + n_to_load = 1 + templateDir = self.reffiles['TemplateDir'] + gammaDir = self.reffiles['GammaDir'] + gammaFile = self.reffiles['GammaFile'] + dustDir = self.reffiles['DustCorrDir'] + + # x1 and color are unique for this simulator + x1 = self.sn_parameters['x1']['min'] + color = self.sn_parameters['color']['min'] + ebvofMW = self.sn_parameters['ebvofMW'] + sn_model = self.simulator_parameters['model'] + sn_version = self.simulator_parameters['version'] + + cutoff = 'cutoff' + if self.error_model: + cutoff = 'error_model' + lcname = 'LC_{}_{}_{}_{}_{}_ebvofMW_{}_vstack.hdf5'.format( + x1, color, cutoff, sn_model, sn_version, ebvofMW) + lcname = 'LC_{}_{}_{}_{}_{}_ebvofMW_{}_vstack.hdf5'.format( + x1, color, cutoff, sn_model, sn_version, ebvofMW) + dustFile = 'Dust_{}_{}_{}.hdf5'.format( + x1, color, cutoff) + + print('loading reference files') + time_ref = time.time() + + result_queue = multiprocessing.Queue() + p = multiprocessing.Process(name='Subprocess-0', + target=self.loadReference, args=( + templateDir, lcname, gammaDir, + gammaFile, web_path, 0, + result_queue)) + p.start() + + if np.abs(ebvofMW) > 1.e-5: + n_to_load = 2 + print('loading dust files') + pb = multiprocessing.Process( + name='Subprocess-1', target=self.loadDust, + args=(dustDir, dustFile, web_path, 1, result_queue)) + pb.start() + + resultdict = {} + for j in range(n_to_load): + resultdict.update(result_queue.get()) + + for p in multiprocessing.active_children(): + p.join() + + reference_lc = resultdict[0] + dustcorr = None + if n_to_load > 1: + dustcorr = resultdict[1] + + return reference_lc, dustcorr + + def load_stacker(self, coadd=False): + """ + Method to load stacker class + + Parameters + ---------- + coadd : bool, optional + to stack or not. The default is False. + + Returns + ------- + stacker : class + CoaddStacker class instance. + + """ + + stacker = None + + if coadd: + stacker = CoaddStacker(col_sum=[self.nexpCol, self.vistimeCol, + 'visitExposureTime'], + col_mean=[self.mjdCol, + self.RACol, self.DecCol, + 'pixRA', + 'pixDec', 'healpixID', + 'season', 'airmass'], + col_median=[ + 'sky', 'moonPhase', 'seeingFwhmEff', + 'lsst_start'], + col_group=[ + self.filterCol, self.nightCol], + col_coadd=self.m5Col, + col_visit='visitExposureTime', col_atmos=[]) + return stacker + + def load_cosmology(self, cosmo_par): + """ + Method to load cosmology parameters + + Parameters + ---------- + cosmo_par : dict + cosmology parameters + + Returns + ------- + cosmology class + + """ + + cosmology = w0waCDM(H0=cosmo_par['H0'], + Om0=cosmo_par['Om'], + Ode0=cosmo_par['Ol'], + w0=cosmo_par['w0'], wa=cosmo_par['wa']) + return cosmology + + def loadReference(self, templateDir, lcname, gammaDir, + gammaFile, web_path, j=-1, output_q=None): + """ + Method to load reference files (lc and gamma) + + Parameters + ---------- + templateDir : str + template dir. + lcname : str + lc reference name. + gammaDir : str + gamma loc dir. + gammaFile : str + gamma file name. + web_path : str + web path of original files. + j : int, optional + internal tag for multiproc. The default is -1. + output_q : multiprocessing queue, optional + queue for multiprocessing. The default is None. + + Returns + ------- + None or multiprocessing queue. + + """ + + reference_lc = GetReference(templateDir, + lcname, gammaDir, gammaFile, web_path) + + if output_q is not None: + output_q.put({j: reference_lc}) + else: + return None + + def loadDust(self, dustDir, dustFile, web_path, j=-1, output_q=None): + """ + Method to load dust files + + Parameters + ---------- + dustDir : str + dust location dir. + dustFile : str + dust file name. + web_path : str + web path of original files. + j : int, optional + internal tag for multiproc. The default is -1. + output_q : multiprocessing queue, optional + queue for multiprocessing. The default is None. + + Returns + ------- + None or multiprocessing queue. + + """ + + dustcorr = LoadDust(dustDir, dustFile, web_path).dustcorr + + if output_q is not None: + output_q.put({j: dustcorr}) + else: + return None + + def load_zp(self, web_path, templateDir, fName): + """ + Method to load zp_airmass results + + Parameters + ---------- + web_path : str + web path of original files. + templateDir : str + location dir + fName : str + file name + + Returns + ------- + res : record array + array with data. + + """ + + from sn_tools.sn_io import check_get_file + check_get_file(web_path, templateDir, fName) + fullName = '{}/{}'.format(templateDir, fName) + + res = np.load(fullName) + + return res + + def load_flux_to_pixel(self, web_path, templateDir, fName): + """ + Method to load the frac flux distrib (for saturation effects) + + Parameters + ---------- + web_path : str + web path of original file. + templateDir : str + location dir of the file. + fName : str + Name of the file to load. + + Returns + ------- + myinterp : scipy.interpolate.interp1d + interpolator (seeing, frac_flux_median) + + """ + + from sn_tools.sn_io import check_get_file + check_get_file(web_path, templateDir, fName) + fullName = '{}/{}'.format(templateDir, fName) + + res = np.load(fullName) + + from scipy.interpolate import interp1d + myinterp = interp1d(res['seeing'], res['pixel_frac_med'], + bounds_error=False, fill_value=0.) + + return myinterp + + +class SNSimulation(SNSimu_Params): + """LC simulation wrapper class + + Parameters + --------------- + + mjdCol: str, opt + mjd col name in observations (default: 'observationStartMJD') + RACol: str, opt + RA col name in observations (default: 'fieldRA') + DecCol:str, opt + Dec col name in observations (default: 'fieldDec') + filterCol: str, opt + filter col name in observations (default: filter') + m5Col: str, opt + 5-sigma depth col name in observations (default: 'fiveSigmaDepth') + exptimeCol: str, opt + exposure time col name in observations (default: 'visitExposureTime) + nightCol: str, opt + night col name in observations (default: 'night') + obsidCol: str, opt + observation id col name in observations (default: 'observationId') + nexpCol: str, opt + number of exposures col name in observations (default: 'numExposures') + visitimeCol: str, opt + visit time col name in observations (default: 'visiTime') + seeingEffCol: str, opt + seeing eff col name in observations (default: 'seeingFwhmEff') + seeingGeomCol: str, opt + seeing geom col name in observations (default: 'seeingFwhmGeom') + coadd: bool, opt + coaddition of obs (per band and per night) if set to True (default: True) + config: dict + configuration dict for simulation (SN parameters, cosmology, telescope,..) + x0_norm: array of float + grid ox (x1,color,x0) values + zp_airmass: dict(interp1d) + zeropoints vs airmass/filter + + """ + + def __init__(self, metricName='SNSimulation', + mjdCol='observationStartMJD', + RACol='fieldRA', DecCol='fieldDec', + filterCol='filter', m5Col='fiveSigmaDepth', + exptimeCol='visitExposureTime', + nightCol='night', obsidCol='observationId', + nexpCol='numExposures', + vistimeCol='visitTime', seeingEffCol='seeingFwhmEff', + airmassCol='airmass', + skyCol='sky', moonCol='moonPhase', + seeingGeomCol='seeingFwhmGeom', + uniqueBlocks=False, config=None, dust_map=None, + zp_airmass=None, **kwargs): + super().__init__(mjdCol=mjdCol, + RACol=RACol, DecCol=DecCol, + filterCol=filterCol, m5Col=m5Col, + exptimeCol=exptimeCol, + nightCol=nightCol, obsidCol=obsidCol, + nexpCol=nexpCol, + vistimeCol=vistimeCol, seeingEffCol=seeingEffCol, + airmassCol=airmassCol, + skyCol=skyCol, moonCol=moonCol, + seeingGeomCol=seeingGeomCol, + config=config, dust_map=dust_map, zp_airmass=zp_airmass) + + def run(self, obs, gen_simu_params, slicePoint=None, imulti=0): + """ LC simulations + + Parameters + -------------- + obs: array + array of observations + + """ + # np.save('obs_orig.npy', obs) + if self.stacker is not None: + obs = self.stacker._run(obs, atmosType=self.atmosType) + # np.save('obs_coadd.npy', obs) + + # set atmos params and throughputs + obs = self.set_atmos_and_throughput(obs) + + par = {} + par['obs'] = obs + par['nspectra'] = self.sn_parameters['nspectra'] + par['lc_coadd'] = self.lc_coadd + + list_lc = self.simuLoop(gen_simu_params, par, imulti) + + if list_lc: + return list_lc + + return None + + def run_deprecated(self, obs, slicePoint=None, imulti=0): """ LC simulations Parameters - -------------- - obs: array - array of observations + -------------- + obs: array + array of observations + + """ + iproc = 1 + + # select filters + goodFilters = np.in1d(obs[self.filterCol], self.filterNames) + obs = obs[goodFilters] + + # estimate seasons + if 'season' not in obs.dtype.names: + obs = seasoncalc(obs, season_gap=50., force_calc=True) + + # check the number of seasons + # if too low get seasons using clusters + + nseasons = len(np.unique(obs['season'])) + + """ + if nseasons <= 8: + print('seasons from cluster') + obs = self.get_season_from_cluster(obs) + """ + """ + delta_max = self.get_delta_per_season(obs) + # obs = rf.drop_fields(obs, 'season') + print('delta_max', delta_max) + obs = seasoncalc(obs, season_gap=65., force_calc=True) + """ + # save these obs. + # np.save('obs_pixel.npy', obs) + + # plot seasons + # self.plot_seasons(obs) + + if len(obs) == 0: + return None + # stack obs if required + # np.save('obs_orig.npy', obs) + if self.stacker is not None: + obs = self.stacker._run(obs, atmosType=self.atmosType) + # np.save('obs_coadd.npy', obs) + self.fieldname = 'unknown' + self.fieldid = 0 + try: + iter(self.season) + except TypeError: + self.season = [self.season] + + if self.season == [-1]: + seasons = np.unique(obs[self.seasonCol]) + else: + seasons = self.season + + time_ref = time.time() + """ + tracemalloc.start() + start = tracemalloc.take_snapshot() + """ + + """ + print('before hhh', len(obs), obs.dtype.names, + np.unique(obs['season'])) + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + ax.plot(obs['observationStartMJD'], obs['airmass'], 'ko') + plt.show() + """ + # select obs corresponding to seasons + idx = np.in1d(obs['season'], seasons) + obs = obs[idx] + + if len(obs) == 0: + return None + + list_lc = [] + """ + for seas in seasons: + idx = obs['season'] == seas + obs_seas = obs[idx] + if len(obs_seas) >= 5: + ll = self.run_season(obs_seas, seas) + if ll is not None: + list_lc += ll + if ll is not None: + list_lc += ll + """ + list_lc = self.run_season(obs, seasons[0]) + + if list_lc: + return list_lc + + return None + + def run_season(self, obs, seas): + """ + Method to simulate LCs per season + + Parameters + ---------- + obs : numpy array + Observations + seas : int + season number. + + Returns + ------- + list_lc : TYPE + DESCRIPTION. + + """ + + # set atmos params and throughputs + obs = self.set_atmos_and_throughput(obs) + + # get simulation parameters + + gen_params = self.get_all_gen_params(obs, [seas]) + + if gen_params is None: + return None + + list_lc = [] + + if gen_params is not None: + print('NLC to simulate:', len(gen_params), + np.unique(obs['healpixID']), len(obs)) + + # LC simulation using multiprocessing + par = {} + par['obs'] = obs + par['nspectra'] = self.sn_parameters['nspectra'] + par['lc_coadd'] = self.lc_coadd + """ + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + ax.plot(obs['observationStartMJD'], obs['filter'], 'ko') + print(np.max(np.diff(obs['observationStartMJD']))) + plt.show() + """ + # print('running', self.nprocs) + list_lc = multiproc(gen_params, par, self.simuLoop, self.nprocs) + + """ + top_stats = snapshot.statistics('lineno') + + print("[ Top 10 ]") + for stat in top_stats[:10]: + print(stat) + """ + + if list_lc: + return list_lc + + return None + + def set_atmos_and_throughput(self, obs): + """ + Function to set atmos params and register throughput + + Parameters + ---------- + obs : numpy array + Data to process. + + Returns + ------- + obs : numpy array + Output data. + + """ + + # select cols on interest only + ccolsa = obs.dtype.names + ccolsb = ['numExposures', 'visitTime', + 'observationStartMJD', 'fieldRA', 'fieldDec', + 'pixRA', 'pixDec', + 'healpixID', 'season', 'sky', 'moonPhase', + 'seeingFwhmEff', 'seeingFwhmGeom', 'note', + 'filter', 'night', 'visitExposureTime', + 'fiveSigmaDepth', 'visitExposureTime', + 'airmass', 'pwv', 'ozone', 'aerosol', 'lsst_start'] + + ccols_ = list(set(ccolsa) & set(ccolsb)) + obs = obs[ccols_] + + # add atmospheric parameters here + obs = self.set_atmos_params(obs) + + # register throughputs for sn_cosmo + idx = obs['airmass'] <= 3.0 + obs = obs[idx] + + if len(obs) < 2: + return None + + obs = self.register_bands_from_atmos(obs) + + # estimate zp and mean_wavelength (and sigmas) corresponding to obs + if self.atmosType == 'const': + obs = self.add_zp_meanwave_from_interp(obs) + else: + obs = self.add_zp_meanwave_from_obs(obs) + + return obs + + def get_season_from_cluster(self, obs): + """ + Method to estimate seasons from clusters + + Parameters + ---------- + obs : numpy array + Data to process. + + Returns + ------- + obs : numpy array + Original data plus season col. + + """ + + from sn_tools.sn_clusters import makeClusters, anaClusters + nclusters = 10 + obs.sort(order=self.mjdCol) + + nobs = len(obs) + if nobs < nclusters: + nclusters = nobs + points, clus, labels = makeClusters( + nclusters, obs, 'pixRA', self.mjdCol) + + dfcluster = anaClusters( + nclusters, obs, points, clus, labels, 'pixRA', self.mjdCol) + + dfcluster = dfcluster.sort_values(by=['pixRA', self.mjdCol]) + dfclmean = dfcluster.groupby( + 'clusId')[self.mjdCol].mean().reset_index() + dfclmean = dfclmean.sort_values(by=[self.mjdCol]).reset_index() + dfclmean['season'] = dfclmean.index+1 + + dfcluster = dfcluster.merge(dfclmean[['clusId', 'season']], + left_on=['clusId'], + right_on=['clusId'], + suffixes=('', '')) + seasons = dfcluster['season'].to_list() + obs.sort(order=['pixRA', self.mjdCol]) + obs = rf.drop_fields(obs, 'season') + obs = rf.append_fields(obs, 'season', seasons) + + return obs + + def get_delta_per_season(self, obs, vara='observationStartMJD'): + """ + Method to get info on data/season + + Parameters + ---------- + obs : numpy array + Data to process. + vara : str, optional + col to get info from. The default is 'observationStartMJD'. + + Returns + ------- + float + Max delta_MJD (over seasons) + + """ + + delta_max = [] + for seas in np.unique(obs['season']): + idx = obs['season'] == seas + sel = obs[idx] + season_length = np.max(sel[vara])-np.min(sel[vara]) + if len(sel) >= 2: + delta_max.append(np.max(np.diff(sel[vara]))) + print(seas, np.max(np.diff(sel[vara])), season_length) + + return np.max(delta_max) + + def plot_seasons(self, obs): + """ + Method to plot seasons + + Parameters + ---------- + obs : array + Data to process. + + Returns + ------- + None. + + """ + + print('seasons', np.unique(obs['season'])) + vara = 'observationStartMJD' + varb = 'fiveSigmaDepth' + import matplotlib.pyplot as plt + plt.plot(obs[vara], + obs[varb], 'ko', mfc='None') + for seas in np.unique(obs['season']): + idx = obs['season'] == seas + sel = obs[idx] + season_length = np.max(sel[vara])-np.min(sel[vara]) + if len(sel) >= 2: + # print(sel[[vara, self.filterCol, 'fieldRA', 'fieldDec']]) + print(seas, np.max(np.diff(sel[vara])), season_length, + np.min(sel[vara]), np.max(sel[vara])) + + plt.plot(sel[vara], + sel[varb], marker='*', linestyle='None') + + plt.show() + + def get_all_gen_params_deprecated(self, obs, seasons): + """ + Method to get simu parameters for all seasons + + Parameters + ---------- + obs: array + observations. + seasons: list(int) + list of seasons + + Returns + ------- + array + simulation parameters. + + """ + + gp = None + for seas in seasons: + + if len(self.simuParamsFile) == 0: + gen_pars = self.gen_params_from_season(obs, seas) + else: + gen_pars = self.gen_params_from_file(obs, seas) + + if gen_pars is None: + continue + + if gp is None: + gp = gen_pars + else: + gp = np.concatenate((gp, gen_pars)) + + return gp + + def gen_params_from_file_deprecated(self, obs, seas): + """ + Method to grab simu parameters from input file + + Parameters + ---------- + obs : numpy array + array of observations. + seas : int + season to process. + + Returns + ------- + sel : numpy array + array of simu parameters. + + """ + + healpixID = np.unique(obs['healpixID']) + + idx = self.simuParamsFile['healpixID'] == healpixID + idx &= self.simuParamsFile['season'] == seas + + sel = self.simuParamsFile[idx] + + return sel + + def gen_params_from_season_deprecated(self, obs, seas): + """ + Method to grab simu params (estimated from obs) for a season + + Parameters + ---------- + obs : numpy array + Observations + seas : int + season of observation. + + Returns + ------- + gen_pars : numpy array + simulation parameters. + + """ + + idxa = obs[self.seasonCol] == seas + obs_season = obs[idxa] + gen_pars = self.gen_par.simuparams(obs_season) + + if gen_pars is None: + return gen_pars + + gen_pars = rf.append_fields(gen_pars, + 'season', + [seas]*len(gen_pars)) + + return gen_pars + + def simuLoop(self, gen_params, params, j=0, output_q=None): + """ + Method to simulate LC (llop, multiproc) + + Parameters + ---------- + gen_params : array + simulation parameters. + params : dict + parameters to use. + j : int, optional + internal tag for multiproc. The default is 0. + output_q : multiprocessing queue, optional + container for output if multiproc used. The default is None. + + Returns + ------- + list(astropy table) + list of simulated LC. + + """ + import time + time_ref = time.time() + # print('multiproc simu', j, len(gen_params)) + obs = params['obs'] + lc_coadd = params['lc_coadd'] + + lsst_start = -1 + if 'lsst_start' in obs.dtype.names: + lsst_start = np.median(obs['lsst_start']) + + nspectra = params['nspectra'] + simu_out, lc_out = None, None + + if self.save_status: + simu_out, lc_out, sed_out = self.prepareSave( + self.outdir, self.prodid, j) + + if 'sn_fast' not in self.simu_config['name']: + lc_list, lc_list_keep, tab_meta, sed_list = self.loop_gen( + obs, gen_params, lc_coadd, j, lc_out) + else: + lc_list = self.simuLCs(obs, gen_params, lc_coadd) + if not self.throwafterdump: + import copy + lc_list_keep = copy.deepcopy(lc_list) + + if self.save_status: + if nspectra <= 0: + if len(lc_list) > 0: + self.dump(lc_list, lc_out) + lc_list = [] + tab_meta.meta['lc_dir'] = self.outdir + tab_meta.meta['lc_fileName'] = lc_out.split('/')[-1] + tab_meta.meta['lsst_start'] = lsst_start + self.write_meta(tab_meta, simu_out) + else: + self.dump_df(tab_meta, simu_out, lc_list, + lc_out, sed_list, sed_out) + + # print('done', j, time.time()-time_ref) + if output_q is not None: + return output_q.put({j: lc_list_keep}) + else: + return lc_list_keep + + def dump_df(self, tab_meta, simu_out, lc_list, lc_out, sed_list, sed_out): + """ + Method to dump output data in pandas df + using format defined by N. Regnault + to account for Spectra production + + Parameters + ---------- + tab_meta : astropy table + Metadata. + simu_out : str + output path (full) for metadata. + lc_list : list(astropy table) + list of light curves. + lc_out : str + output path (full) for lc. + sed_list : list(astropy table) + list of spectra. + sed_out : str + output path (full) for spectra. + + Returns + ------- + None. + + """ + import pandas as pd + + # meta data or Sn_data + df_meta = tab_meta.to_pandas() + rename_dict = dict(zip(['SNID', 'daymax', 'color', 'ebvofMW'], [ + 'sn', 'tmax', 'col', 'ebv'])) + df_meta = df_meta.rename(columns=rename_dict) + df_meta['valid'] = 1 + df_meta['IAU'] = 0 + cols = ['sn', 'z', 'tmax', 'x0', 'x1', 'col', 'ebv', 'valid', 'IAU'] + df_meta[cols].to_hdf(simu_out, key='sn_data') + # print(df_meta[cols]) + + # light curves + df_lc = pd.DataFrame() + for io, lc in enumerate(lc_list): + lca = lc.to_pandas() + snid = lc.meta['SNID'] + lca['sn'] = snid + lca['lc'] = 'lc_{}'.format(io) + df_lc = pd.concat((df_lc, lca)) + + df_lc['magsys'] = 'AB' + df_lc['valid'] = 1 + rename_dict = dict(zip(['time', 'sky', 'seeingFwhmEff'], + ['mjd', 'mag_sky', 'seeing'])) + + df_lc = df_lc.rename(columns=rename_dict) + cols = ['sn', 'mjd', 'flux', 'fluxerr', 'band', 'magsys', 'exptime', + 'valid', 'lc', 'zp', 'mag_sky', 'seeing'] + df_lc[cols].to_hdf(lc_out, key='lc_data') + + # spectra + df_spectra = pd.DataFrame() + for sed in sed_list: + df_spectra = pd.concat((df_spectra, sed.to_pandas())) + + cols = ['sn', 'mjd', 'wavelength', 'flux', + 'fluxerr', 'valid', 'spec', 'exptime'] + + df_spectra.to_hdf(sed_out, key='spec_data') + + def loop_gen(self, obs, gen_params, lc_coadd, j, lc_out): + """ + Method to generate LCs by looping on genparams + + Parameters + ---------- + obs : array + data to process (observations). + gen_params : array + simulation parameters. + lc_coadd: int + to coadd lc fluxes. + j : int + tag for SNID and output file. + lc_out : str + output file name. + + Returns + ------- + lc_list : TYPE + DESCRIPTION. + lc_list_keep : TYPE + DESCRIPTION. + tab_meta : TYPE + DESCRIPTION. + + """ + + lc_list = [] + lc_list_keep = [] + sed_list = [] + isn = 0 + tab_meta = Table() + + for genpar in gen_params: + isn += 1 + season = genpar['season'] + idx = obs['season'] == season + obs_season = obs[idx] + lc, sed = self.simuLCs(obs_season, genpar, lc_coadd) + if len(lc) == 0: + continue + + lc = lc[0] + + hpix = int(np.mean(obs['healpixID'])) + isn_str = str(isn) + + if 'SNID' in genpar.dtype.names: + sn_id = genpar['SNID'] + else: + sn_id = 'SN_{}_{}_{}_{}'.format( + str(hpix).zfill(7), str(season).zfill(2), isn_str.zfill(5), j) + + lc.meta['SNID'] = sn_id + lc_list += [lc] + if sed: + sed = sed[0] + sed['sn_id'] = sn_id + sed_list += [sed] + + if not self.throwafterdump: + lc_list_keep += [lc] + # every xx SN: dump to file + if self.save_status: + if len(lc_list) >= 100: + self.dump(lc_list, lc_out) + lc_list = [] + tab_meta = vstack([tab_meta, Table(rows=[lc.meta])]) + + return lc_list, lc_list_keep, tab_meta, sed_list + + def prepareSave(self, outdir, prodid, iproc): + """ Prepare output directories for data + + Parameters + -------------- + outdir: str + output directory where to copy the data + prodid: str + production id(label for input files) + iproc: int + internal tag for multiprocessing + + Returns + ---------- + Two output files are open: + - astropy table with (SNID, RA, Dec, X1, Color, z) parameters + -> name: Simu_prodid.hdf5 + - astropy tables with LC + -> name: LC_prodid.hdf5 + + """ + + if not os.path.exists(outdir): + print('Creating output directory', outdir) + os.makedirs(outdir) + # Two files to be opened - tagged by iproc + # One containing a summary of the simulation: + # astropy table with (SNID,RA,Dec,X1,Color,z) parameters + # -> name: Simu_prodid.hdf5 + # A second containing the Light curves (list of astropy tables) + # -> name : LC_prodid.hdf5 + + simu_out = '{}/Simu_{}_{}.hdf5'.format( + outdir, prodid, iproc) + lc_out = '{}/LC_{}_{}.hdf5'.format(outdir, prodid, iproc) + sed_out = '{}/Spectra_{}_{}.hdf5'.format(outdir, prodid, iproc) + + if self.clean_status: + self.check_del(simu_out) + self.check_del(lc_out) + self.check_del(sed_out) + + return simu_out, lc_out, sed_out + + def check_del(self, fileName): + """ + Method to remove a file if already exist + + Parameters + ---------- + fileName: str + file to remove(full path) + + """ + if os.path.exists(fileName): + os.remove(fileName) + + def simuLCs(self, obs, gen_params, lc_coadd): + """ Generate LC for one season and a set of simu parameters + + Parameters + -------------- + obs: array + array of observations + gen_params: dict + generation parameters + lc_coadd: int. + to coadd lc fluxes. + + Returns + ---------- + lc_table: astropy table + table with LC informations(flux, time, ...) + """ + + sn_par = self.sn_parameters.copy() + # simulator_par = self.simulator_parameters.copy() + + # print('oooo',self.sn_parameters) + sn_par['sigmaz'] = self.sn_parameters['z']['sigmaz'] + for name in ['z', 'x1', 'color', 'daymax']: + if name in gen_params.dtype.names: + sn_par[name] = gen_params[name] + + SNID = sn_par['Id'] + + # print('aoooo',self.sn_parameters) + # print(test) + sn_object = SN_Object(self.simu_config['name'], + sn_par, + self.simulator_parameters, + gen_params, + self.cosmology, + SNID, self.area, + x0_grid=None, + salt2Dir=self.salt2Dir, + mjdCol=self.mjdCol, + RACol=self.RACol, + DecCol=self.DecCol, + filterCol=self.filterCol, + exptimeCol=self.exptimeCol, + m5Col=self.m5Col, + psf_flux=self.psf_flux, + frac_flux_seeing=self.frac_flux_seeing, + ccd_full_well=self.ccd_full_well) + + module = import_module(self.simu_config['name']) + + simu = module.SN(sn_object, self.simu_config, + self.reference_lc, self.dustcorr, + dust_map=self.dust_map) + # simulation - this is supposed to be a list of astropytables + lc_table = simu(obs, self.display_lc, self.time_display, lc_coadd) + + seds = [] + nspectra = self.sn_parameters['nspectra'] + if nspectra > 0: + seds = simu.SN_SED(gen_params, nspectra=nspectra) + del simu + del module + return lc_table, seds + + def set_atmos_params(self, obsa): + """ + Method to set atmospheric parameters + + Parameters + ---------- + obs: numpy array + Data to process. + + Returns + ------- + obs : pandas df + output data. + + """ + import numpy.lib.recfunctions as rf + from random import gauss + + obs = np.copy(obsa) + del obsa + + vvt = {} + for atm_param in ['pwv', 'ozone', 'aerosol']: + atm_value = self.config_instr[atm_param] + atm_round_value = self.config_instr['round'][atm_param] + atm_value = np.round(atm_value, atm_round_value) + vvt[atm_param] = [atm_value]*len(obs) + for myval in ['sigma', 'min', 'max', 'round']: + atm_val_value = self.config_instr[myval][atm_param] + vvt['{}_{}'.format(myval, atm_param)] = [ + atm_val_value]*len(obs) + kk = list(vvt.keys()) + vv = list(vvt.values()) + obs = rf.append_fields(obs, kk, vv) + """ + for atm_param in ['pwv', 'ozone', 'aerosol']: + if atm_param not in obs.dtype.names: + atm_value = self.config_instr[atm_param] + atm_round_value = self.config_instr['round'][atm_param] + atm_value = np.round(atm_value, atm_round_value) + obs = rf.append_fields(obs, atm_param, [atm_value]*len(obs)) + + for myval in ['sigma', 'min', 'max', 'round']: + atm_val_value = self.config_instr[myval][atm_param] + obs = rf.append_fields(obs, '{}_{}'.format(myval, + atm_param), [atm_val_value]*len(obs)) + """ + """ + # atm_value = eval('self.{}'.format(atm_param)) + # atm_round_value = eval('self.{}_round'.format(atm_param)) + atm_value = np.round(atm_value, atm_round_value) + obs = rf.append_fields(obs, atm_param, [atm_value]*len(obs)) + # obs[atm_param] = [atm_value]*len(obs) + obs = rf.append_fields(obs, 'sigma_{}'.format( + atm_param), [atm_sigma_value]*len(obs)) + """ + # smear atmospheric parameters + # vv_sigma = [0.0]*len(obs) + """ + if self.atmosType != 'const': + vv_sigma = [ + eval('self.sigma_{}'.format(atm_param))]*len(obs) + obs[atm_param] += np.random.normal(0, vv_sigma) + """ + """ + sigma_val = eval('self.sigma_{}'.format(atm_param)) + obs = rf.append_fields( + obs, 'sigma_{}'.format(atm_param), [sigma_val]*len(obs)) + """ + + # airmass + vvtb = {} + for myval in ['sigma', 'min', 'max', 'round']: + val = self.config_instr[myval]['airmass'] + vvtb['{}_airmass'.format(myval)] = [val]*len(obs) + kk = list(vvtb.keys()) + vv = list(vvtb.values()) + + obs = rf.append_fields(obs, kk, vv) + + """ + for myval in ['sigma', 'min', 'max', 'round']: + val = self.config_instr[myval]['airmass'] + obs = rf.append_fields( + obs, '{}_airmass'.format(myval), [val]*len(obs)) + """ + return obs + + def add_zp_meanwave_from_interp(self, obs): + """ + Method to estimate mean_restframe wavelength + + Parameters + ---------- + obs : TYPE + DESCRIPTION. + + Returns + ------- + TYPE + DESCRIPTION. + + """ + + filters = np.array(obs[self.filterCol].tolist()) + airmass = np.array(obs['airmass'].tolist()) + + filt_airmass = list(zip(filters, airmass)) + mean_wave = list( + map(lambda x: self.zp_atmos['mean_wave'][x[0]](x[1]), filt_airmass)) + sigma_mean_wave = list( + map(lambda x: self.zp_atmos['sigma_mean_wave'][x[0]](x[1]), filt_airmass)) + + zp = list(map(lambda x: self.zp_atmos['zp'][x[0]](x[1]), filt_airmass)) + sigma_zp = list( + map(lambda x: self.zp_atmos['sigma_zp'][x[0]](x[1]), filt_airmass)) + + import numpy.lib.recfunctions as rf + obs = rf.append_fields(obs, 'zp', zp) + obs = rf.append_fields(obs, 'sigma_zp', sigma_zp) + obs = rf.append_fields(obs, 'mean_wave', mean_wave) + obs = rf.append_fields(obs, 'sigma_mean_wave', sigma_mean_wave) + obs = rf.append_fields(obs, 'tel_site_name', [ + self.telescope.site_name]*len(obs)) + + return obs + + def add_zp_meanwave_from_obs(self, obs): + """ + Method to estimate zero points + + Parameters + ---------- + obs : record array + data to process. + + Returns + ------- + obs : numpy array + output result (original array with atmos. params) + + """ + + ra = [] + rb = [] + + for row in obs: + airmass = row['airmass'] + pwv = row['pwv'] + ozone = row['ozone'] + aerosol = row['aerosol'] + b = row['filter'] + exptime = row['visitExposureTime'] + nexp = row['numExposures'] + + self.telescope.new_atmosphere(site_name=self.telescope.site_name, + airmass=airmass, + aerosol=aerosol, + pwv=pwv, ozone=ozone) + self.telescope.reset_data() + self.telescope.mean_wave() + # grab zp + mean_wave = self.telescope.mean_wavelength[b] + zp = self.telescope.zp(b, exptime, nexp) + ra.append((zp)) + rb.append((mean_wave)) + + import numpy.lib.recfunctions as rf + obs = rf.append_fields(obs, 'zp', ra) + obs = rf.append_fields(obs, 'mean_wave', rb) + obs = rf.append_fields(obs, 'tel_site_name', [ + self.telescope.site_name]) + return obs + + def register_bands_from_atmos(self, obs, + cols=['airmass', + 'pwv', 'ozone', 'aerosol']): + """ + Method to register throughputs from atmos params + + Parameters + ---------- + obs : numpy array + observations. + cols : list(str), optional + columns of interest. The default is ['airmass','pwv', 'ozone', 'aerosol']. + + Returns + ------- + obs : numpy array + observations. """ - iproc = 1 - if 'iproc' in obs.dtype.names: - iproc = int(np.mean(obs['iproc'])) + # round atmos parameters + for i, vv in enumerate(cols): + # round_value = eval('self.{}_round'.format(vv)) + round_value = self.config_instr['round'][vv] + obs[vv] = np.round(obs[vv], round_value) - if iproc not in self.nprocdict: - self.nprocdict[iproc] = iproc - if self.save_status: - self.prepareSave(self.outdir, self.prodid, iproc) + # set band_cosmo column + bcols = self.telescope.site_name+'::' + \ + obs[self.filterCol]+'_' +\ + obs[self.airmassCol].astype(str)+'_' +\ + obs['pwv'].astype(str)+'_' +\ + obs['ozone'].astype(str)+'_' +\ + obs['aerosol'].astype(str) - # if 'healpixID' not in obs.dtype.names: - if slicePoint is not None: - import numpy.lib.recfunctions as rf - healpixID = hp.ang2pix( - slicePoint['nside'], np.rad2deg(slicePoint['ra']), np.rad2deg(slicePoint['dec']), nest=True, lonlat=True) - pixRA, pixDec = hp.pix2ang( - self.nside, healpixID, nest=True, lonlat=True) - obs = rf.append_fields(obs, 'healpixID', [healpixID]*len(obs)) - obs = rf.append_fields(obs, 'pixRA', [pixRA]*len(obs)) - obs = rf.append_fields(obs, 'pixDec', [pixDec]*len(obs)) + import numpy.lib.recfunctions as rf + obs = rf.append_fields(obs, 'band_cosmo', bcols.tolist()) - print('processing pixel', healpixID) + # grab unique parameters + all_cols = [self.filterCol]+['band_cosmo']+cols + atmos_table = unique(Table(obs[all_cols])) - # estimate seasons - obs = seasoncalc(obs) - # stack if necessary - if self.stacker is not None: - obs = self.stacker._run(obs) - # obs = Observations(data=tab, names=self.names) - self.fieldname = 'unknown' - self.fieldid = 0 - if self.season == -1: - seasons = np.unique(obs[self.seasonCol]) - else: - seasons = self.season + # register - time_ref = time.time() - for seas in seasons: - self.index_hdf5 += 10000*(seas-1) + self.register_bands_on_the_fly(self.telescope, + atmos_table[all_cols].to_pandas()) - idxa = obs[self.seasonCol] == seas - obs_season = obs[idxa] - # remove the u band - idx = [i for i, val in enumerate( - obs_season[self.filterCol]) if val[-1] != 'u'] + return obs - if len(obs_season[idx]) >= 5: - self.simuSeason(obs_season[idx], seas, iproc) + def register_bands_on_the_fly(self, telescope, data): + """ + Method to register bands on sncosmo - # save metadata - self.save_metadata(np.unique(obs['healpixID']).item()) - # reset metadata dict - self.sn_meta[iproc] = {} + Parameters + ---------- + telescope: Telescope + telescope to use + data: pandas df + data to register - #print('End of simulation', time.time()-time_ref) + Returns + ------- + None. - def prepareSave(self, outdir, prodid, iproc): - """ Prepare output directories for data + """ + import sncosmo + for i, row in data.iterrows(): + bandname = row['band_cosmo'] + band = row[self.filterCol] + airmass = row[self.airmassCol] + pwv = row['pwv'] + ozone = row['ozone'] + aerosol = row['aerosol'] + register_bands_sncosmo(sncosmo, telescope, + bandname, band, + airmass, pwv, ozone, aerosol) + + def dump(self, lc_list, lc_out): + """ + Method to dum lc on file Parameters - -------------- - outdir: str - output directory where to copy the data - prodid: str - production id (label for input files) - iproc: int - internal tag for multiprocessing - + ---------- + list_lc : list(astropytable) + data to dump. + lc_out : str + outputfile name. Returns - ---------- - Two output files are open: - - astropy table with (SNID,RA,Dec,X1,Color,z) parameters - -> name: Simu_prodid.hdf5 - - astropy tables with LC - -> name : LC_prodid.hdf5 + ------- + None. """ + for lc in lc_list: + astropy.io.misc.hdf5.write_table_hdf5( + lc, lc_out, path=lc.meta['SNID'], + append=True, serialize_meta=True) - if not os.path.exists(outdir): - print('Creating output directory', outdir) - os.makedirs(outdir) - # Two files to be opened - tagged by iproc - # One containing a summary of the simulation: - # astropy table with (SNID,RA,Dec,X1,Color,z) parameters - # -> name: Simu_prodid.hdf5 - # A second containing the Light curves (list of astropy tables) - # -> name : LC_prodid.hdf5 + def write_meta(self, meta, out_meta): + """ + Method to dum lc on file - self.simu_out[iproc] = '{}/Simu_{}_{}.hdf5'.format( - outdir, prodid, iproc) - self.lc_out[iproc] = '{}/LC_{}_{}.hdf5'.format(outdir, prodid, iproc) - self.check_del(self.simu_out[iproc]) - self.check_del(self.lc_out[iproc]) - self.SNID[iproc] = 10**iproc - self.sn_meta[iproc] = {} + Parameters + ---------- + list_lc : list(astropytable) + data to dump. + lc_out : str + outputfile name. + + Returns + ------- + None. - """ - # and these files will be removed now (before processing) - # if they exist (to avoid confusions) - if os.path.exists(self.simu_out): - os.remove(self.simu_out) - if os.path.exists(self.lc_out): - os.remove(self.lc_out) """ - def check_del(self, fileName): + path = 'meta_{}'.format(int(np.mean(meta['healpixID']))) + + astropy.io.misc.hdf5.write_table_hdf5( + meta, out_meta, path=path, + append=True, serialize_meta=False) + + def setIndex(self, healpixID, x1, color, z, daymax, + season, epsilon, SNID): """ - Method to remove a file if already exist + Method to set specific index Parameters ---------- - fileName: str - file to remove (full path) + healpixID : int + healpixID. + x1 : float + SN x1. + color : float + SN color. + z : float + SN redshift. + daymax : float + SN T0. + season : int + season num. + epsilon : float + epsilon for SN simu params. + SNID : str + SNID. + + Returns + ------- + index_hdf5 : str + resulting index. """ - if os.path.exists(fileName): - os.remove(fileName) - def simuSeason(self, obs, season, iproc): - """ Generate LC for a season (multiprocessing available) and all simu parameters + index_hdf5 = '{}_{}_{}_{}_{}'.format( + healpixID, z, daymax, season, SNID) + + if x1 != 'undef': + index_hdf5 += '_{}_{}'.format(x1, color) + + # epsilon should be last!! + index_hdf5 += '_{}'.format(epsilon) + + return index_hdf5 + + def simuSeason_deprecated(self, obs, season, iproc): + """ Generate LC for a season(multiprocessing available) + and all simu parameters Parameters -------------- @@ -367,29 +1692,53 @@ def simuSeason(self, obs, season, iproc): """ - gen_params = self.gen_par.Params(obs) + gen_params = self.gen_par.simuparams(obs) if gen_params is None: return - #self.simuLoop(obs, season, gen_params, iproc) + # self.simuLoop(obs, season, gen_params, iproc) npp = self.nprocs if 'sn_fast' in self.simu_config['name']: npp = 1 + list_lc = [] + + if npp == 1: + metadict, list_lc = self.simuLoop(obs, season, gen_params, iproc) + if self.save_status: + if not self.sn_meta[iproc]: + self.sn_meta[iproc] = metadict + else: + # self.sn_meta[iproc]= self.sn_meta[iproc].update(metadict) + for key in metadict.keys(): + self.sn_meta[iproc][key] += metadict[key] + + else: + list_lc = self.multiSeason(obs, season, gen_params, iproc, npp) + + if len(list_lc): + return list_lc + + return None + + def multiSeason_deprecated(self, obs, season, gen_params, iproc, npp): + nlc = len(gen_params) batch = np.linspace(0, nlc, npp+1, dtype='int') - + print('batch for seasons', npp, batch, season) + result_queue = multiprocessing.Queue() for i in range(npp): ida = batch[i] idb = batch[i+1] - - p = multiprocessing.Process(name='Subprocess', target=self.simuLoop, args=( - obs, season, gen_params[ida:idb], iproc, i, result_queue)) + p = multiprocessing.Process(name='Subprocess', + target=self.simuLoop, args=( + obs, season, gen_params[ida:idb], + iproc, i, result_queue)) p.start() resultdict = {} @@ -399,15 +1748,19 @@ def simuSeason(self, obs, season, iproc): for p in multiprocessing.active_children(): p.join() + list_lc = [] for j in range(npp): - metadict = resultdict[j] + metadict = resultdict[j][0] + list_lc += resultdict[j][1] + """ if not self.sn_meta[iproc]: self.sn_meta[iproc] = metadict else: - #self.sn_meta[iproc]= self.sn_meta[iproc].update(metadict) + # self.sn_meta[iproc]= self.sn_meta[iproc].update(metadict) for key in metadict.keys(): self.sn_meta[iproc][key] += metadict[key] - + """ + return list_lc # self.save_metadata() """ SNID = 100 @@ -425,7 +1778,31 @@ def simuSeason(self, obs, season, iproc): self.writeLC(SNID, lc, season) """ - def writeLC(self, SNID, lc, season, iproc, meta_lc): + def complete_meta_deprecated(self, lc, meta_lc): + + epsilon = 0. + if 'epsilon_x0' in lc.meta.keys(): + epsilon = np.int(1000*1.e8*lc.meta['epsilon_x0']) + epsilon += np.int(100*1.e8*lc.meta['epsilon_x1']) + epsilon += np.int(10*1.e8*lc.meta['epsilon_color']) + epsilon += np.int(1*1.e8*lc.meta['epsilon_daymax']) + + if 'x1' not in lc.meta.keys(): + x1 = 'undef' + color = 'undef' + else: + x1 = lc.meta['x1'] + color = lc.meta['color'] + + SNID_tot = '{}_{}'.format( + lc.meta['sn_type'], lc.meta['healpixID']) + + index_hdf5 = SNID_tot + meta_lc['SNID'] = SNID_tot + + return meta_lc + + def writeLC_deprecated(self, SNID, lc, season, iproc, meta_lc): """ Method to save lc on disk and to update metadata @@ -444,18 +1821,43 @@ def writeLC(self, SNID, lc, season, iproc, meta_lc): """ # save LC on disk - epsilon = np.int(1000*1.e8*lc.meta['epsilon_x0']) - epsilon += np.int(100*1.e8*lc.meta['epsilon_x1']) - epsilon += np.int(10*1.e8*lc.meta['epsilon_color']) - epsilon += np.int(1*1.e8*lc.meta['epsilon_daymax']) + epsilon = 0. + if 'epsilon_x0' in lc.meta.keys(): + epsilon = np.int(1000*1.e8*lc.meta['epsilon_x0']) + epsilon += np.int(100*1.e8*lc.meta['epsilon_x1']) + epsilon += np.int(10*1.e8*lc.meta['epsilon_color']) + epsilon += np.int(1*1.e8*lc.meta['epsilon_daymax']) - index_hdf5 = self.setIndex(lc.meta['healpixID'], - lc.meta['x1'], - lc.meta['color'], - np.round(lc.meta['z'], 3), - np.round(lc.meta['daymax'], 3), - season, epsilon) + if 'x1' not in lc.meta.keys(): + x1 = 'undef' + color = 'undef' + else: + x1 = lc.meta['x1'] + color = lc.meta['color'] + SNID_tot = '{}_{}_{}_{}'.format( + lc.meta['sn_type'], lc.meta['healpixID'], iproc, SNID) + """ + index_hdf5 = self.setIndex(lc.meta['healpixID'], + x1, color, + np.round(lc.meta['z'], 4), + np.round(lc.meta['daymax'], 4), + season, epsilon, SNID) + """ + index_hdf5 = SNID_tot + lc.meta['SNID'] = SNID_tot + """ + idx = lc['snr_m5'] > 0. + lc = lc[idx] + lc.meta = {} + """ + # print('writing',lc,lc.meta) + """ + lc.write(self.lc_out[iproc], + path='lc_{}'.format(index_hdf5), + append=True, + compression=True, serialize_meta=True) + """ lc.write(self.lc_out[iproc], path='lc_{}'.format(index_hdf5), append=True, @@ -463,11 +1865,18 @@ def writeLC(self, SNID, lc, season, iproc, meta_lc): # build metadata dict n_lc_points = len(lc) + """ metanames = ['SNID', 'index_hdf5', 'season', 'fieldname', 'fieldid', 'n_lc_points', 'area'] metavals = [SNID, index_hdf5, season, self.fieldname, self.fieldid, n_lc_points, self.area] + """ + metanames = ['index_hdf5', 'season', + 'fieldname', 'fieldid', 'n_lc_points', 'area', 'lcName'] + metavals = [index_hdf5, season, + self.fieldname, self.fieldid, + n_lc_points, self.area, self.lc_out[iproc]] metadict = dict(zip(metanames, metavals)) metadict.update(lc.meta) @@ -488,17 +1897,8 @@ def writeLC(self, SNID, lc, season, iproc, meta_lc): for key in metadict.keys(): meta_lc[key].extend([metadict[key]]) - def setIndex(self, healpixID, x1, color, z, daymax, season, epsilon): - - index_hdf5 = '{}_{}_{}_{}_{}_{}_{}'.format(healpixID, - x1, - color, - z, - daymax, season, epsilon) - - return index_hdf5 - - def simuLoop(self, obs, season, gen_params, iproc, j=0, output_q=None): + def simuLoop_deprecated(self, obs, season, gen_params, + iproc, j=0, output_q=None): """ Method to simulate LC @@ -516,34 +1916,42 @@ def simuLoop(self, obs, season, gen_params, iproc, j=0, output_q=None): internal parameter for multiprocessing output_q: multiprocessing.Queue() - """ time_ref = time.time() list_lc = [] + list_lc_keep = [] meta_lc = {} + if 'sn_fast' not in self.simu_config['name']: for genpar in gen_params: lc = self.simuLCs(obs, season, genpar) if lc: list_lc += lc + if not self.throwafterdump: + list_lc_keep += lc # every 20 SN: dump to file - if len(list_lc) >= 20: - self.dump(list_lc, season, iproc, meta_lc) - list_lc = [] + if self.save_status: + if len(list_lc) >= 20: + self.dump(list_lc, season, iproc, meta_lc) + list_lc = [] else: list_lc = self.simuLCs(obs, season, gen_params) + if not self.throwafterdump: + import copy + list_lc_keep = copy.deepcopy(list_lc) - if len(list_lc) > 0: + if len(list_lc) > 0 and self.save_status: self.dump(list_lc, season, iproc, meta_lc) + list_lc = [] if output_q is not None: - output_q.put({j: meta_lc}) + output_q.put({j: (meta_lc, list_lc_keep)}) else: - return meta_lc + return (meta_lc, list_lc_keep) - def dump(self, list_lc, season, j, meta_lc): + def dump_deprecated(self, list_lc, season, j, meta_lc): """ Method to write a list of lc on disk @@ -557,67 +1965,175 @@ def dump(self, list_lc, season, j, meta_lc): tag for multiprocessing """ + for lc in list_lc: - self.SNID[j] += 1 - self.writeLC(self.SNID[j], lc, season, j, meta_lc) + ido = True + if self.throwaway_empty and len(lc) == 0: + ido = False + if ido: + self.SNID[j] += 1 + self.writeLC(self.SNID[j], lc, season, j, meta_lc) + + def save_metadata_deprecated(self, isav=-1): + """ Copy metadata to disk - def simuLCs(self, obs, season, gen_params): - """ Generate LC for one season and a set of simu parameters + """ + if self.sn_meta: + for key, vals in self.sn_meta.items(): + if vals: + # print('metadata',vals) + Table(vals).write( + self.simu_out[key], 'summary_{}'.format(isav), + append=True, compression=True) + + def prepareSave_deprecated(self, outdir, prodid, iproc): + """ Prepare output directories for data Parameters -------------- - obs: array - array of observations - season: int - season number - gen_params: dict - generation parameters - ex: - index_hdf5: int - SN index in hdf5 output file - output_q: dict - output for multiprocessing + outdir: str + output directory where to copy the data + prodid: str + production id(label for input files) + iproc: int + internal tag for multiprocessing Returns ---------- - lc_table: astropy table - table with LC informations (flux, time, ...) - metadata: dict - metadata of the simulation + Two output files are open: + - astropy table with (SNID, RA, Dec, X1, Color, z) parameters + -> name: Simu_prodid.hdf5 + - astropy tables with LC + -> name: LC_prodid.hdf5 + """ - sn_par = self.sn_parameters.copy() - for name in ['z', 'x1', 'color', 'daymax']: - sn_par[name] = gen_params[name] + if not os.path.exists(outdir): + print('Creating output directory', outdir) + os.makedirs(outdir) + # Two files to be opened - tagged by iproc + # One containing a summary of the simulation: + # astropy table with (SNID,RA,Dec,X1,Color,z) parameters + # -> name: Simu_prodid.hdf5 + # A second containing the Light curves (list of astropy tables) + # -> name : LC_prodid.hdf5 - SNID = sn_par['Id'] - sn_object = SN_Object(self.simu_config['name'], - sn_par, - gen_params, - self.cosmology, - self.telescope, SNID, self.area, - x0_grid=self.x0_grid, - salt2Dir=self.salt2Dir, - mjdCol=self.mjdCol, - RACol=self.RACol, - DecCol=self.DecCol, - filterCol=self.filterCol, - exptimeCol=self.exptimeCol, - m5Col=self.m5Col) + self.simu_out[iproc] = '{}/Simu_{}_{}.hdf5'.format( + outdir, prodid, iproc) + self.lc_out[iproc] = '{}/LC_{}_{}.hdf5'.format(outdir, prodid, iproc) + self.check_del(self.simu_out[iproc]) + self.check_del(self.lc_out[iproc]) + self.SNID[iproc] = 10**iproc + self.sn_meta[iproc] = {} - module = import_module(self.simu_config['name']) - simu = module.SN(sn_object, self.simu_config, - self.reference_lc, self.gamma, self.mag_to_flux, self.dustcorr,error_model=self.error_model) - # simulation - this is supposed to be a list of astropytables - lc_table = simu(obs, self.display_lc, self.time_display) + """ + # and these files will be removed now (before processing) + # if they exist (to avoid confusions) + if os.path.exists(self.simu_out): + os.remove(self.simu_out) + if os.path.exists(self.lc_out): + os.remove(self.lc_out) + """ + + def run_deprecated(self, obs, slicePoint=None, imulti=0): + """ LC simulations - return lc_table + Parameters + -------------- + obs: array + array of observations - def save_metadata(self, isav=-1): - """ Copy metadata to disk + """ + iproc = 1 + + if 'iproc' in obs.dtype.names: + iproc = int(np.mean(obs['iproc'])) + + if iproc not in self.nprocdict: + self.nprocdict[iproc] = iproc + if self.save_status: + self.prepareSave(self.outdir, self.prodid, iproc) + + # if 'healpixID' not in obs.dtype.names: + if slicePoint is not None: + import numpy.lib.recfunctions as rf + healpixID = hp.ang2pix( + slicePoint['nside'], np.rad2deg(slicePoint['ra']), + np.rad2deg(slicePoint['dec']), nest=True, lonlat=True) + pixRA, pixDec = hp.pix2ang( + self.nside, healpixID, nest=True, lonlat=True) + obs = rf.append_fields(obs, 'healpixID', [healpixID]*len(obs)) + obs = rf.append_fields(obs, 'pixRA', [pixRA]*len(obs)) + obs = rf.append_fields(obs, 'pixDec', [pixDec]*len(obs)) + + print('processing pixel', healpixID) + + # estimate seasons + obs = seasoncalc(obs) + + # select filters + goodFilters = np.in1d(obs[self.filterCol], self.filterNames) + obs = obs[goodFilters] + + np.save('obs1.npy', obs) + # stack if necessary + if self.stacker is not None: + obs = self.stacker._run(obs) + + np.save('obs_coadd.npy', obs) + self.fieldname = 'unknown' + self.fieldid = 0 + try: + iter(self.season) + except TypeError: + self.season = [self.season] + + if self.season == [-1]: + seasons = np.unique(obs[self.seasonCol]) + else: + seasons = self.season + + time_ref = time.time() + """ + tracemalloc.start() + start = tracemalloc.take_snapshot() + """ + list_lc = [] + + for seas in seasons: + self.index_hdf5 += 10000*(seas-1) + + idxa = obs[self.seasonCol] == seas + obs_season = obs[idxa] + + if len(obs_season) >= 5: + print('obs in season', len(obs_season), seas, iproc) + simres = self.simuSeason(obs_season, seas, iproc) + if simres is not None: + list_lc += simres + + """ + current = tracemalloc.take_snapshot() + stats = current.compare_to(start, 'filename') + for i, stat in enumerate(stats[:5], 1): + print("since_start", i, str(stat)) + """ + # save metadata + if self.save_status: + self.save_metadata(np.unique(obs['healpixID']).item()) + # reset metadata dict + self.sn_meta[iproc] = {} """ - for key, vals in self.sn_meta.items(): - if vals: - Table(vals).write( - self.simu_out[key], 'summary_{}'.format(isav), append=True, compression=True) + top_stats = snapshot.statistics('lineno') + + print("[ Top 10 ]") + for stat in top_stats[:10]: + print(stat) + """ + # print('End of simulation', time.time()-time_ref) + + if list_lc: + return list_lc + + return None diff --git a/sn_simu_wrapper/sn_wrapper_for_simu.py b/sn_simu_wrapper/sn_wrapper_for_simu.py new file mode 100644 index 00000000..b9ecef17 --- /dev/null +++ b/sn_simu_wrapper/sn_wrapper_for_simu.py @@ -0,0 +1,1489 @@ +from sn_simu_wrapper.sn_simu import SNSimulation +import numpy as np +import yaml +import os +from sn_tools.sn_io import check_get_file +import pandas as pd +import operator +from astropy.table import Table, vstack +from sn_tools.sn_utils import load_config +from sn_tools.sn_obs import load_season +from sn_telmodel.sn_transtools import zp_from_config +from sn_fit_wrapper.sn_wrapper_for_fit import FitWrapper +import warnings +warnings.filterwarnings("ignore", category=FutureWarning) + + +class MakeYaml: + """ + class to generate a yaml file from a generic one + + Parameters + --------------- + dbDir: str + location dir of the database + dbName: str + OS name + db Extens: str + db extension (npy or db) + nside: int + nside for healpix + nproc: int + number of proc for multiprocessing + diffflux: bool + to allow for simulation with differential params (ex: x1+epsilon_x1) + seasnum: list(int) + season numbers + outDir: str + output directory for the production (and also for this yaml file) + fieldType: str + type of the field to process (DD, WFD, fake) + x1Type: str + x1 type for simulation (unique, uniform, random) + x1min: float + x1 min value + x1max: float + x1 max value + x1step: float + x1 step value + colorType: str + color type for simulation (unique, uniform, random) + colormin: float + color min value + colormax: float + color max value + colorstep: float + color step value + zType: str + z type for simulation (unique, uniform, random) + zmin: float + z min value + zmax: float + z max value + zstep: float + z step value + simu: str + simulator type + daymaxType: str + daymax type for simulation (unique, uniform, random) + daymaxstep: float + daymax step value + coadd: bool + to coadd (True) or not (Fals) observations per night + prodid: str + production id ; the resulting yaml file is prodid.yaml + ebvmw: float + to specify an extinction value + bluecutoff: float + blue cutoff for SN + redcutoff: float + redcutoff for SN + error_model: int + error model for flux error estimation + """ + + def __init__(self, dbDir, dbName, dbExtens, nside, nproc, diffflux, + seasnum, outDir, fieldType, + x1Type, x1min, x1max, x1step, + colorType, colormin, colormax, colorstep, + zType, zmin, zmax, zstep, + simu, daymaxType, daymaxstep, + coadd, prodid, + ebvofMW, bluecutoff, redcutoff, error_model): + + self.dbDir = dbDir + self.dbName = dbName + self.dbExtens = dbExtens + self.nside = nside + self.nproc = nproc + self.diffflux = diffflux + self.seasnum = seasnum + self.outDir = outDir + self.fieldType = fieldType + self.x1Type = x1Type + self.x1min = x1min + self.x1max = x1max + self.x1step = x1step + self.colorType = colorType + self.colormin = colormin + self.colormax = colormax + self.colorstep = colorstep + self.zmin = zmin + self.zmax = zmax + self.zstep = zstep + self.simu = simu + self.zType = zType + self.daymaxType = daymaxType + self.daymaxstep = daymaxstep + self.coadd = coadd + self.prodid = prodid + self.ebvofMW = ebvofMW + self.bluecutoff = bluecutoff + self.redcutoff = redcutoff + self.error_model = error_model + + def genYaml(self, input_file): + """ + method to generate a yaml file + with parameters from generic input_file + + Parameters + --------------- + input_file: str + input generic yaml file + + Returns + ----------- + yaml file with parameters + + + """ + with open(input_file, 'r') as file: + filedata = file.read() + + fullDbName = '{}/{}.{}'.format(self.dbDir, self.dbName, self.dbExtens) + filedata = filedata.replace('prodid', self.prodid) + filedata = filedata.replace('fullDbName', fullDbName) + filedata = filedata.replace('nnproc', str(self.nproc)) + filedata = filedata.replace('nnside', str(self.nside)) + filedata = filedata.replace('outputDir', self.outDir) + filedata = filedata.replace('diffflux', str(self.diffflux)) + filedata = filedata.replace('seasval', str(self.seasnum)) + filedata = filedata.replace('ftype', self.fieldType) + filedata = filedata.replace('x1Type', self.x1Type) + filedata = filedata.replace('x1min', str(self.x1min)) + filedata = filedata.replace('x1max', str(self.x1max)) + filedata = filedata.replace('x1step', str(self.x1step)) + filedata = filedata.replace('colorType', self.colorType) + filedata = filedata.replace('colormin', str(self.colormin)) + filedata = filedata.replace('colormax', str(self.colormax)) + filedata = filedata.replace('colorstep', str(self.colorstep)) + filedata = filedata.replace('zmin', str(self.zmin)) + filedata = filedata.replace('zmax', str(self.zmax)) + filedata = filedata.replace('zstep', str(self.zstep)) + filedata = filedata.replace('zType', self.zType) + filedata = filedata.replace('daymaxType', self.daymaxType) + filedata = filedata.replace('daymaxstep', str(self.daymaxstep)) + filedata = filedata.replace('fcoadd', str(self.coadd)) + filedata = filedata.replace('mysimu', self.simu) + filedata = filedata.replace('ebvofMWval', str(self.ebvofMW)) + filedata = filedata.replace('bluecutoffval', str(self.bluecutoff)) + filedata = filedata.replace('redcutoffval', str(self.redcutoff)) + filedata = filedata.replace('errmod', str(self.error_model)) + + return yaml.load(filedata, Loader=yaml.FullLoader) + + +class FitWrapper_deprecated: + def __init__(self, yaml_config_fit): + """ + Class to fit a set of light curves + + Parameters + ---------- + config_fit : dict + parameters fot + + Returns + ------- + None. + + """ + from sn_fit.process_fit import Fitting + + # Fit instance + config = load_config(yaml_config_fit) + + self.fit = Fitting(config) + self.nproc = config['MultiprocessingFit']['nproc'] + + self.saveData = config['OutputFit']['save'] + + self.outDir = config['OutputFit']['directory'] + + self.prodid = config['Simulations']['prodid'] + + if self.saveData: + from sn_tools.sn_io import checkDir + checkDir(self.outDir) + + def __call__(self, lc_list, remove_sat=False): + """ + Main fit method using multiprocessing + + Parameters + ---------- + lc_list : list(lc) + List of light curves to fit. + remove_sat : bool, optional + To remove saturated points. The default is False. + + Returns + ------- + res : pandas df + output results. + + """ + + res = self.fit.fit_multiproc(lc_list, remove_sat, self.nproc) + + return res + + def __call__deprecated(self, lc_list, remove_sat=False): + """ + Method to fit light curves + + Parameters + ---------- + lc_list : list(astropy table) + LC to fit + remove_sat : bool, optional + To remove saturated fluxes. The default is False. + + Returns + ------- + None. + + """ + """ + from astropy.table import Table, vstack + res = Table() + for lc in lc_list: + lc.convert_bytestring_to_unicode() + resfit = self.fit(lc) + if resfit is not None: + res = vstack([res, resfit]) + + return res + """ + # from sn_tools.sn_utils import multiproc + params = {} + params['remove_sat'] = remove_sat + + res = self.multiproc(lc_list, params, self.fit_lcs, self.nproc) + + return res + + def fit_lcs_deprecated(self, lc_list, params, j=0, output_q=None): + """ + Method to fit LCs + + Parameters + ---------- + lc_list : list(astropy table) + light-curves to fit. + params : dict + parameters. + j : int, optional + Tag for multiprocessing. The default is 0. + output_q : multiprocessing queue, optional + queue managing multiprocessing run. The default is None. + + Returns + ------- + astropytable + Result of the fit. + + """ + + from astropy.table import Table, vstack + res = Table() + # print('processing fit', j) + + for lc in lc_list: + lc.convert_bytestring_to_unicode() + resfit = self.fit(lc, params) + if resfit is not None: + resfit = self.check_correct(resfit) + res = vstack([res, resfit]) + + if output_q is not None: + return output_q.put({j: res}) + else: + return res + + def check_correct(self, sn): + """ + Method to correct for Cov_xy col names + + Parameters + ---------- + sn : astropy Table + Data to process. + + Returns + ------- + sn : astropy Table + Processed data + + """ + + varlist = ['z', 't0', 'x0', 'x1', 'color'] + + if 'Cov_zz' not in sn.columns: + varlist = ['t0', 'x0', 'x1', 'color'] + + for i, namea in enumerate(varlist): + for j, nameb in enumerate(varlist): + if j >= i: + vva = 'Cov_{}{}'.format(namea, nameb) + vvb = 'Cov_{}{}'.format(nameb, namea) + if vva not in sn.columns: + sn.rename_column(vvb, vva) + + return sn + + +class InfoWrapper: + def __init__(self, confDict): + """ + class to estimate global parameters of LC + and add a selection flag according to selection values in dict + + Parameters + ---------- + confDict : dict + parameters for selection + + Returns + ------- + None. + + """ + + self.nproc = confDict['nproc_sel'] + from astropy.table import Table + selfile = confDict['selection_params'] + selpars = Table.read(selfile, format='csv', guess=False, comment='#') + + self.snr_min_value = 0 + self.snr_min_op = operator.ge + idx = selpars['selname'] == 'snr_min' + selb = selpars[idx] + if len(selb) > 0: + self.snr_min_value = selb['selval'][0] + self.snrmin_op = selb['selop'][0] + selpars = selpars[~idx] + + self.selparams = selpars + + def __call__(self, light_curves): + """ + Main method to estimate LC shepe params + and add a flag for selection + + Parameters + ---------- + light_curves : list of astropytables + LC curves to process + + Returns + ------- + None. + + """ + + conf_names = ['n_epochs_bef', + 'n_epochs_aft', + 'n_epochs_phase_minus_10', + 'n_epochs_phase_plus_20', + 'n_epochs_m10_p35', + 'n_epochs_m10_p5', + 'n_epochs_p5_p20', + 'n_bands_m8_p10'] + + configs = [('night', 'phase', operator.ge, -20, operator.le, 0), + ('night', 'phase', operator.gt, 0, operator.le, 60), + ('night', 'phase', operator.ge, -100, operator.le, -10.), + ('night', 'phase', operator.gt, 20., operator.le, 100.), + ('night', 'phase', operator.ge, -10., operator.le, 35.), + ('night', 'phase', operator.ge, -10., operator.le, 5.), + ('night', 'phase', operator.ge, 5., operator.le, 20.), + ('band', 'phase', operator.ge, -8, operator.le, 10.)] + + getInfos = dict(zip(conf_names, configs)) + + """ + selParams = [('n_epochs_m10_p35', operator.ge, 4), + ('n_epochs_m10_p5', operator.ge, 1), + ('n_epochs_p5_p20', operator.ge, 1), + ('n_bands_m8_p10', operator.ge, 2)] + """ + from sn_tools.sn_utils import multiproc + params = {} + params['getInfos'] = getInfos + # params['selParams'] = selParams + + lc_list = multiproc(light_curves, params, self.run_list, self.nproc) + + return lc_list + + def run_list(self, light_curves, params, j, output_q=None): + """ + method to estimate general params of the light curve + + Parameters + ---------- + light_curves : list(Table) + Light curve list. + params : dict + Parameters. + j : int + internal parameter (multiprocessing). + output_q : multiprocessing queue, optional + Where to put the results. The default is None. + + Returns + ------- + TYPE + DESCRIPTION. + + """ + + getInfos = params['getInfos'] + # selParams = params['selParams'] + + lc_list = [] + snr_max = [2, 5, 10, 15, 20] + + for lc in light_curves: + T0 = lc.meta['daymax'] + z = lc.meta['z'] + + if len(lc) == 0: + resdict = self.calc_dummy(getInfos, snr_max) + else: + # apply SNR selection + idx = self.snr_min_op(lc['snr'], self.snr_min_value) + lc_sel = lc[idx] + if len(lc_sel) == 0: + resdict = self.calc_dummy(getInfos, snr_max) + else: + resdict = self.calc_infos(lc_sel, T0, z, + getInfos, selParams={}) + for vval in snr_max: + vc = 'Nfilt_{}'.format(vval) + resdict[vc] = self.nfilt_snrmax(lc_sel, snr_max=vval) + + # update meta data + + lc.meta.update(resdict) + lc_list.append(lc) + + if output_q is not None: + return output_q.put({j: lc_list}) + else: + return lc_list + + def calc_dummy(self, getInfos, snr_max): + """ + Method returning dummy infos + + Parameters + ---------- + getInfos : dict + dict of selection criteria to measure. + snr_max: list(int). + snr_max values for Nfilt estimation + + Returns + ------- + None. + + """ + + resdict = {} + for key in getInfos.keys(): + resdict[key] = 0 + + for b in 'ugrizy': + resdict['SNR_{}'.format(b)] = -1.0 + + resdict['SNR'] = -1.0 + resdict['selected'] = 0 + for vv in snr_max: + resdict['Nfilt_{}'.format(vv)] = 0 + + return resdict + + def calc_infos(self, lc_sel, T0, z, getInfos, selParams={}): + """ + Method returning infos related to getInfos + + Parameters + ---------- + lc_sel : astropy table + LC. + T0 : float + SN daymax. + z : float + SN z. + getInfos : dict + dict of selection criteria to measure. + selParams: dict, opt. + selection parameters. The default is {}. + + Returns + ------- + None. + + """ + + resdict = {} + # add phase column + lc_sel['phase'] = (lc_sel['time']-T0)/(1+z) + """ + if 'filter' in lc_sel.columns: + lc_sel.remove_columns(['filter']) + """ + # self.plotLC(lc_sel) + for key, vals in getInfos.items(): + resdict[key] = self.nepochs_phase( + lc_sel, vals[0], vals[1], vals[2], + vals[3], vals[4], vals[5]) + + if selParams: + resdict['selected'] = int(self.select(resdict, selParams)) + # add snr per band + SNRtot = 0. + + lc_sel = lc_sel.to_pandas() + + for b in 'ugrizy': + idx = lc_sel['band'].str.contains('LSST::{}'.format(b)) + sel = lc_sel[idx] + + SNR = 0. + if len(sel) > 0: + SNR = np.sum(sel['snr_m5']**2) + SNRtot += SNR + resdict['SNR_{}'.format(b)] = np.sqrt(SNR) + + resdict['SNR'] = SNRtot + + del lc_sel + return resdict + + def select(self, res, list_sel): + """ + Method to estimate if a LC passes the cut or not + + Parameters + ---------- + dictval : dict + dict of values + + Returns + ------- + bool decision (1= selected, 0=not selected) + + """ + + idx = True + for vals in list_sel: + idx &= vals[1](res[vals[0]], vals[2]) + if not idx: + return idx + + return idx + """ + for key, vals in dictval.items(): + idx = self.selparams['selname'] == key + pp = self.selparams[idx] + if len(pp) > 0: + op = pp['selop'][0] + selval = pp['selval'][0] + selstr = '{}({},{})'.format(op, vals, selval) + resu = eval(selstr) + if not resu: + return False + + return True + """ + + def nepochs_phase(self, tab, colnum='night', + colsel='phase', opa=operator.ge, vala=0, + opb=operator.le, valb=10): + """ + Method to get the number of epochs between two vals vala and valb + + Parameters + ---------- + tab : astropy table + data to process + colnum : str, optional + column to extract the number of epochs from. + The default is 'night'. + colsel : str, optional + selection column name. The default is 'phase'. + opa : operator, optional + operator to apply. The default is operator.ge. + vala : float, optional + selection value. The default is 0. + opb: operator, optional + operator to apply. The default is operator.le. + valb : float, optional + selection value. The default is 10. + + Returns + ------- + TYPE + DESCRIPTION. + + """ + idx = opa(tab[colsel], vala) + idx &= opb(tab[colsel], valb) + tt = tab[idx] + + res = len(np.unique(tt[colnum])) + + return res + + def nfilt_snrmax(self, lc, snr_max=10): + """ + Method to estimate the number of bands with max SNR >= snr_max + + Parameters + ---------- + lc : astropy table + data to process. + snr_max : float, optional + SNR max value. The default is 10. + + Returns + ------- + TYPE + DESCRIPTION. + + """ + + bands = np.unique(lc['band']) + + r = [] + for b in bands: + idx = lc['band'] == b + sel = lc[idx] + lc_snr_max = np.max(sel['snr_m5']) + if lc_snr_max >= snr_max: + r.append(b) + + return len(r) + + def plotLC(self, tab): + """ + Method to plot LC for cross-checks + + Parameters + ---------- + tab : astropy table + data to process + + Returns + ------- + None. + + """ + """ + from sn_simu_wrapper.sn_object import SN_Object + SN_Object.plotLC(tab,time_display) + """ + import matplotlib.pyplot as plt + plt.plot(tab['phase'], tab['flux_e_sec'], 'ko') + plt.show() + + +class SimInfoFitWrapper: + def __init__(self, yaml_config_simu, + infoDict, + yaml_config_fit, + fit_remove_sat): + """ + + + Parameters + ---------- + yaml_config_simu : yaml file + config file for simulation + infoDict : dict + DESCRIPTION. + + Returns + ------- + None. + + """ + + self.name = 'sim_info_fit' + config_simu = load_config(yaml_config_simu) + # self.config_instr = config_simu['InstrumentSimu'] + self.config_simu = config_simu + self.infoDict = infoDict + self.yaml_config_fit = yaml_config_fit + self.fit_remove_sat_str = fit_remove_sat = fit_remove_sat + """ + self.simu_wrapper = SimuWrapper(yaml_config_simu) + self.info_wrapper = InfoWrapper(infoDict) + self.fit_wrapper = FitWrapper(yaml_config_fit) + self.fit_remove_sat = list(map(int, fit_remove_sat.split(','))) + """ + self.outName = '' + + self.ccolref = [] + fw = load_config(self.yaml_config_fit) + + if fw['OutputFit']['save']: + simu_wrapper = load_config(self.config_simu) + outFile = 'SN_{}.hdf5'.format(simu_wrapper['ProductionIDSimu']) + self.outName = '{}/{}'.format(fw['OutputFit'] + ['directory'], outFile) + # check wether this file already exist and remove it + # import os + if os.path.isfile(self.outName): + os.system('rm {}'.format(self.outName)) + del simu_wrapper + del fw + self.outdf = pd.DataFrame() + + # grab seasons + # grab seasons + + self.seasons = load_season(config_simu['Observations']['season']) + + # getting zeropoints vs airmass + self.zp_atmos = zp_from_config(config_simu['InstrumentSimu']) + + # info required for obs_quality + min_rf_phase_qual = self.config_simu['SN']['minRFphaseQual'] + max_rf_phase_qual = self.config_simu['SN']['maxRFphaseQual'] + self.diff_rf = max_rf_phase_qual-min_rf_phase_qual + self.zmin = self.config_simu['SN']['z']['min'] + + def instances(self): + """ + Method to instantiate necessary classes + + Parameters + ---------- + None + + Returns + ------- + None. + + """ + + self.simu_wrapper = SimuWrapper(self.config_simu, self.zp_atmos) + self.info_wrapper = InfoWrapper(self.infoDict) + self.fit_wrapper = FitWrapper(self.yaml_config_fit) + self.fit_remove_sat = list( + map(int, self.fit_remove_sat_str.split(','))) + + def run(self, obs, imulti=0, verbose=False): + """ + Parameters + ---------- + obs : array + array of observations + imulti : int, optional + internal tag. The default is 0. + verbose : bool, optional + To print infos. The default is False. + + Returns + ------- + None. + + """ + if 'season' not in obs.dtype.names: + from sn_tools.sn_obs import season as seasoncalc + obs = seasoncalc(obs, season_gap=50., force_calc=True) + + for i, seas in enumerate(self.seasons): + idx = obs['season'] == seas + obs_seas = obs[idx] + if not self.obs_quality(obs_seas): + print('bad obs quality', seas, len(obs_seas)) + continue + + # update config for simu + self.config_simu['Observations']['season'] = '{}'.format(seas) + + # instances + + self.instances() + + # grab SNe Ia simulation parameters + """ + gen_simu_params = simu_params(obs, [seas], + self.simu_wrapper.simuParamsFile, + self.simu_wrapper.gen_par) + """ + gen_simu_params = self.simu_wrapper.simu_par_gen(obs, seas) + + if gen_simu_params is None: + continue + + print('Number of LC to generate', len(gen_simu_params)) + # run + params = {} + params['obs'] = obs_seas + params['verbose'] = False + params['imulti'] = imulti + # self.run_season_new(obs_seas, gen_simu_params,imulti) + from sn_tools.sn_utils import multiproc + # simulate LCs + import time + time_ref = time.time() + light_curves = multiproc(gen_simu_params, params, + self.run_season_simulc, + self.simu_wrapper.nproc) + + if verbose: + print('finished', len(light_curves), time.time()-time_ref) + print('LC analysis') + + light_curves_ana = self.info_wrapper(light_curves) + + # fitting here + if verbose: + print('lc fitting', len(light_curves_ana)) + + for rr in self.fit_remove_sat: + fitlc = self.fit_wrapper(light_curves_ana, remove_sat=rr) + if len(fitlc) > 0: + fitlc['remove_sat'] = rr + self.myconcat(fitlc) + + if verbose: + print('end of fitting', time.time()-time_ref) + + def run_season_simulc(self, gen_params, params, j=0, output_q=None): + + obs = params['obs'] + verbose = params['verbose'] + imulti = params['imulti'] + if verbose: + import time + time_ref = time.time() + print('simulation') + + light_curves = self.simu_wrapper(obs, gen_params, j) + + # print('done', j, time.time()-time_ref) + + if output_q is not None: + return output_q.put({j: light_curves}) + else: + return light_curves + + def run_season_deprecated(self, obs, imulti=0, verbose=True): + """ + + Parameters + ---------- + obs : array + array of observations + imulti : int, optional + internal tag. The default is 0. + verbose : bool, optional + To print infos. The default is False. + + Returns + ------- + None. + + """ + + # get Light curves from simuWrapper + # print('processing pixel', np.unique( + # obs[['healpixID', 'pixRA', 'pixDec']])) + if verbose: + import time + time_ref = time.time() + print('simulation') + + light_curves = self.simu_wrapper(obs, imulti) + + if verbose: + print('after simulation', light_curves) + # analyze these LC + flag for selection + if light_curves is None: + return None + + # light_curves = self.myanatest(light_curves) + + if verbose: + print('LC analysis') + + light_curves_ana = self.info_wrapper(light_curves) + + """ + ccols = ['n_epochs_phase_minus_10', 'n_epochs_phase_plus_20', + 'n_epochs_m10_p35', 'n_epochs_m10_p5', + 'n_epochs_p5_p20', 'n_bands_m8_p10'] + + for gg in light_curves_ana: + for col in ccols: + diff = gg.meta[col]-gg.meta['{}_n'.format(col)] + diff = int(diff) + if diff != 0: + print(col, diff, gg.meta[col], gg.meta['{}_n'.format(col)]) + print('---') + """ + # print('nlc analyzed', len(light_curves_ana)) + + # fitting here + if verbose: + print('lc fitting', len(light_curves_ana)) + + for rr in self.fit_remove_sat: + fitlc = self.fit_wrapper(light_curves_ana, remove_sat=rr) + if len(fitlc) > 0: + fitlc['remove_sat'] = rr + self.myconcat(fitlc) + + if verbose: + print('end of fitting', time.time()-time_ref) + + if len(self.outdf) > 10000: + self.dump_df() + self.outdf = pd.DataFrame() + # ccol = ['RA', 'Dec', 'sn_type'] + + # print('nsn', len(fitlc), time.time()-time_ref) + del light_curves + del light_curves_ana + del fitlc + + return None + """ + if self.fit_wrapper.saveData: + outFile = 'SN_{}.hdf5'.format(self.simu_wrapper.prodid) + outName = '{}/{}'.format(self.fit_wrapper.outDir,outFile) + import astropy + astropy.io.misc.hdf5.write_table_hdf5(fitlc, outName, + path='SN', overwrite=True, + serialize_meta=True) + """ + + def obs_quality(self, obs, + duration_min_z=20, nobs_min=5, + mjdCol='observationStartMJD'): + """ + Method to estimate obs quality (nobs and season length) + + Parameters + ---------- + obs : pandas df + Data to check. + duration_min_z : float, optional + min season duration req (z dep). The default is 20. + nobs_min : int, optional + min number of observations. The default is 5. + mjdCol : str, optional + mjd colname. The default is 'observationStartMJD'. + + Returns + ------- + bool + True=ok; False=not ok. + + """ + + if len(obs) < nobs_min: + return False + + daymin = np.min(obs[mjdCol]) + daymax = np.max(obs[mjdCol]) + duration = daymax-daymin + + zlim = (duration-duration_min_z)/self.diff_rf + zlim -= 1 + + if zlim < self.zmin: + return False + + return True + + def myanatest(self, lightcurves): + """ + Method to estimate some values as cross-check + + Parameters + ---------- + lightcurves : list(astropytables) + LC list. + + Returns + ------- + rr : list(astropytables) + LC list (metadata updated). + + """ + + rr = [] + for lc in lightcurves: + dd = {} + daymax = lc.meta['daymax'] + z = lc.meta['z'] + lc['phase'] = (lc['time']-daymax)/(1.+z) + idx = lc['snr'] >= 1. + lc = Table(lc[idx]) + + idx = lc['phase'] <= -10 + sel = lc[idx] + dd['n_epochs_phase_minus_10_n'] = len(np.unique(sel['night'])) + idx = lc['phase'] >= 20 + sel = lc[idx] + dd['n_epochs_phase_plus_20_n'] = len(np.unique(sel['night'])) + idx = lc['phase'] >= -10 + idx &= lc['phase'] <= 35 + sel = lc[idx] + dd['n_epochs_m10_p35_n'] = len(np.unique(sel['night'])) + idx = lc['phase'] >= -10 + idx &= lc['phase'] <= 5 + sel = lc[idx] + dd['n_epochs_m10_p5_n'] = len(np.unique(sel['night'])) + idx = lc['phase'] >= 5 + idx &= lc['phase'] <= 20 + sel = lc[idx] + dd['n_epochs_p5_p20_n'] = len(np.unique(sel['night'])) + idx = lc['phase'] >= -8 + idx &= lc['phase'] <= 10 + sel = lc[idx] + nnb = 0 + if len(sel) > 0: + nnb = len(np.unique(sel['filter'])) + dd['n_bands_m8_p10_n'] = nnb + for key, vals in dd.items(): + lc.meta[key] = vals + + rr.append(lc) + return rr + + def dump(self, fitlc): + """ + + + Parameters + ---------- + fitlc : astropyTable + data to dump + + Returns + ------- + None. + + """ + """ + if self.outName != '': + keyhdf = '{}'.format(int(sn['healpixID'].mean())) + sn.write(self.outName, keyhdf, append=True, compression=True) + """ + if self.fit_wrapper.saveData: + fitlc.convert_bytestring_to_unicode() + df = pd.DataFrame(fitlc.to_pandas()) + + if 'selected' in df.columns: + df = df.drop(columns=['selected']) + + if not self.ccolref: + self.ccolref = df.columns.to_list() + else: + df = df.reindex(columns=self.ccolref) + + """ + print('chisq', df['chisq']) + for vv in df.columns: + print(vv, df[vv].dtype) + """ + """ + for vv in self.ccolref: + print(vv, df[vv].dtype) + """ + """ + cols = ['sn_type', 'sn_model', 'sn_version', 'fitstatus'] + + print(df['fitstatus'].unique(), df['SNID'].unique()) + """ + + print('dumping', len(df)) + df.to_hdf(self.outName, key='SN', append=True) + + def dump_df(self): + """ + Method to dum a pandaf df to a file + + Returns + ------- + None. + + """ + + # print('dumping df', len(self.outdf)) + self.outdf.to_hdf(self.outName, key='SN', append=True) + + def myconcat(self, fitlc): + """ + Method to concat a set of astropy tables with a pandas df + + Parameters + ---------- + fitlc : astropy table + Data to concat. + + Returns + ------- + None. + + """ + + fitlc.convert_bytestring_to_unicode() + df = pd.DataFrame(fitlc.to_pandas()) + + if 'selected' in df.columns: + df = df.drop(columns=['selected']) + + if not self.ccolref: + self.ccolref = df.columns.to_list() + else: + df = df.reindex(columns=self.ccolref) + + self.outdf = pd.concat((self.outdf, df)) + + def finish(self): + """ + Method to use at the end of the run + + Returns + ------- + None. + + """ + + if len(self.outdf) > 0: + self.dump_df() + + +class SimuWrapper: + """ + Wrapper class for simulation + + Parameters + --------------- + yaml_config: str + name of the yaml configuration file + + """ + + def __init__(self, config, zp_airmass, mjdCol='observationStartMJD'): + + # config = load_config(yaml_config) + + self.saveData_simu = config['OutputSimu']['savefromwrapper'] + self.lc_out = None + self.meta_table = Table() + self.meta_out = None + if self.saveData_simu: + from sn_tools.sn_io import checkDir + outDir = config['OutputSimu']['directory'] + checkDir(outDir) + prodid = config['ProductionIDSimu'] + lcpath = '{}/{}.hdf5'.format(outDir, prodid) + metapath = '{}/{}.hdf5'.format(outDir, + prodid.replace('LC', 'Simu')) + self.lc_out = lcpath + self.meta_out = metapath + self.outDir = outDir + + if os.path.isfile(lcpath): + os.system('rm {}'.format(lcpath)) + if os.path.isfile(metapath): + os.system('rm {}'.format(metapath)) + self.name = 'simulation' + + # get X0 for SNIa normalization + x0_tab = None + x0_griddata = config['SN']['x0']['griddata'] + if x0_griddata: + x0_tab = self.x0(config) + + # load references if simulator = sn_fast + # reference_lc = self.load_reference(config) + + # now define the metric instance + # self.metric = SNMAFSimulation(config=config, x0_norm=x0_tab, + # reference_lc=reference_lc, + # coadd=config['Observations']['coadd']) + + self.metric = SNSimulation( + config=config, x0_norm=x0_tab, zp_airmass=zp_airmass) + + self.prodid = config['ProductionIDSimu'] + self.outlc = [] + + # gen simu parameters instance + + import healpy as hp + nside = config['Pixelisation']['nside'] + area = hp.nside2pixarea(nside, degrees=True) + from sn_tools.sn_utils import SN_simu_params + self.simu_par_gen = SN_simu_params(config['SN'], config['Cosmology'], + mjdCol=mjdCol, area=area, + web_path=config['WebPathSimu']) + + # check if the dust map is available in reference_files + # if not: load it from web + # if not available: consider producing it! + fName = 'reference_files/dustmap_{}.hdf5'.format(nside) + if not os.path.isfile(fName): + self.getRefFile( + config['WebPathSimu'], 'reference_files', 'dustmap_{}.hdf5'.format(nside)) + + dust_map = pd.DataFrame() + if not os.path.isfile(fName): + print('File', fName, 'not found') + print('You should consider using the following script to gen it') + print('python run_scripts/dust_for_fast/gen_disp_dustmap.py') + else: + dust_map = pd.read_hdf(fName) + self.metric = SNSimulation( + config=config, dust_map=dust_map, zp_airmass=zp_airmass) + """ + # simu params from file + from sn_tools.sn_utils import simu_params_from_file + self.simuParamsFile = simu_params_from_file(config['SN']) + """ + + self.nproc = config['MultiprocessingSimu']['nproc'] + + def getRefFile(self, web_path, refdir, fname): + """ + Method to get a file from the web + + Parameters + ---------- + web_path : str + web path. + refdir : str + Dir of reference (whre the file is supposed to be). + fname : str + File name. + + Returns + ------- + None. + + """ + fullname = '{}/dust_maps/{}'.format(web_path, fname) + + # check whether the file is available; if not-> get it! + if not os.path.isfile(fname): + print('wget path:', fullname) + cmd = 'wget --no-clobber --no-verbose {} --directory-prefix {}'.format( + fullname, refdir) + os.system(cmd) + + def x0(self, config): + """ + Method to load x0 data + + Parameters + --------------- + config: dict + parameters to load and (potentially) regenerate x0s + + Returns + ----------- + + """ + # check whether X0_norm file exist or not + # (and generate it if necessary) + absMag = config['SN']['absmag'] + x0normFile = 'reference_files/X0_norm_{}.npy'.format(absMag) + if not os.path.isfile(x0normFile): + # if this file does not exist, grab it from a web server + check_get_file(config['WebPathSimu'], 'reference_files', + 'X0_norm_{}.npy'.format(absMag)) + + if not os.path.isfile(x0normFile): + # if the file could not be found, then have to generate it! + salt2Dir = config['SN']['salt2Dir'] + model = config['Simulator']['model'] + version = str(config['Simulator']['version']) + + # need the SALT2 dir for this + from sn_tools.sn_io import check_get_dir + check_get_dir(config['Web path'], 'SALT2', salt2Dir) + from sn_tools.sn_utils import X0_norm + X0_norm(salt2Dir=salt2Dir, model=model, version=version, + absmag=absMag, outfile=x0normFile) + + return np.load(x0normFile) + + def run(self, obs, gen_simu_params, imulti=0): + """ + Method to run the metric + + Parameters + --------------- + obs: array + data to process + + """ + + light_curves = self.metric.run(obs, gen_simu_params, imulti=imulti) + + if light_curves is None: + return None + + if len(light_curves) == 0: + return None + + if self.saveData_simu: + self.increment_data(light_curves) + if len(self.outlc) > 100: + self.dump_lc() + self.outlc = [] + + return None + + return light_curves + + __call__ = run + + def increment_data(self, light_curves): + """ + Method to build metadata table + + Parameters + ---------- + light_curves : Table + light curves. + + Returns + ------- + None. + + """ + self.outlc += light_curves + for ll in light_curves: + tt = Table(rows=[list(ll.meta.values())], + names=list(ll.meta.keys())) + tt.meta['lc_dir'] = self.outDir + tt.meta['lc_fileName'] = self.lc_out.split('/')[-1] + self.meta_table = vstack([self.meta_table, tt]) + + def dump_lc(self): + """ + Method to dum a pandaf df to a file + + Returns + ------- + None. + + """ + import astropy + + for lc in self.outlc: + astropy.io.misc.hdf5.write_table_hdf5( + lc, self.lc_out, path=lc.meta['SNID'], + append=True, serialize_meta=True) + + def finish(self): + """ + Method to use at the end of the run + + Returns + ------- + None. + + """ + import astropy + if len(self.outlc) > 0: + self.dump_lc() + + if len(self.meta_table) > 0: + self.meta_table.meta['lc_dir'] = self.outDir + self.meta_table.meta['lc_fileName'] = self.lc_out.split('/')[-1] + astropy.io.misc.hdf5.write_table_hdf5( + self.meta_table, self.meta_out, path='metadata', + append=True, serialize_meta=True) + + +class InfoFitWrapper: + def __init__(self, infoDict, yaml_config_fit): + """ + + + Parameters + ---------- + yaml_config_simu : yaml file + config file for simulation + infoDict : dict + DESCRIPTION. + + Returns + ------- + None. + + """ + self.name = 'info_fit' + self.info_wrapper = InfoWrapper(infoDict) + + self.outName = '' + + if self.fit_wrapper.saveData: + outFile = 'SN_{}.hdf5'.format(self.fit_wrapper.prodid) + self.outName = '{}/{}'.format(self.fit_wrapper.outDir, outFile) + # check wether this file already exist and remove it + import os + if os.path.isfile(self.outName): + os.system('rm {}'.format(self.outName)) + + self.yaml_config_fit = yaml_config_fit + + def run(self, light_curves): + """ + + + Parameters + ---------- + light_curves : list of astropy table + LC to process + + Returns + ------- + None. + + """ + + # analyze these LC + flag for selection + light_curves_ana = self.info_wrapper(light_curves) + print('nlc analyzed bis', len(light_curves_ana)) + + # fitting here + self.fit_wrapper = FitWrapper(self.yaml_config_fit) + fitlc = self.fit_wrapper(light_curves_ana) + + self.dump(fitlc) + + return fitlc + + def dump(self, sn): + """ + + + Parameters + ---------- + sn : astropyTable + data to dump + + Returns + ------- + None. + + """ + if self.outName != '': + keyhdf = '{}'.format(int(sn['healpixID'].mean())) + sn.write(self.outName, keyhdf, append=True, compression=True) diff --git a/sn_simulator/sn_cosmo.py b/sn_simulator/sn_cosmo.py index 3052e38c..10407afc 100644 --- a/sn_simulator/sn_cosmo.py +++ b/sn_simulator/sn_cosmo.py @@ -1,32 +1,43 @@ import sncosmo import numpy as np -from lsst.sims.photUtils import Bandpass, Sed -from lsst.sims.photUtils import SignalToNoise -from lsst.sims.photUtils import PhotometricParameters -from astropy.table import Table, Column -from lsst.sims.catUtils.dust import EBV -from scipy.interpolate import griddata -import h5py +from astropy.table import Table, vstack, unique +from scipy.interpolate import interp1d from sn_simu_wrapper.sn_object import SN_Object import time from sn_tools.sn_utils import SNTimer -from sn_tools.sn_calcFast import srand import pandas as pd -import operator -from astropy import units as u +import os +# from dustmaps.sfd import SFDQuery +# from random import gauss + class SN(SN_Object): - def __init__(self, param, simu_param, reference_lc=None, gamma=None, mag_to_flux=None, dustcorr=None, snr_fluxsec='interp',error_model=True): - super().__init__(param.name, param.sn_parameters, param.gen_parameters, - param.cosmology, param.telescope, param.SNID, param.area, param.x0_grid, - mjdCol=param.mjdCol, RACol=param.RACol, DecCol=param.DecCol, - filterCol=param.filterCol, exptimeCol=param.exptimeCol, + def __init__(self, param, simu_param, + reference_lc=None, dustcorr=None, dust_map=None): + super().__init__(param.name, + param.sn_parameters, + param.simulator_parameters, + param.gen_parameters, + param.cosmology, + param.SNID, + param.area, param.x0_grid, + mjdCol=param.mjdCol, + RACol=param.RACol, + DecCol=param.DecCol, + filterCol=param.filterCol, + exptimeCol=param.exptimeCol, nexpCol=param.nexpCol, - m5Col=param.m5Col, seasonCol=param.seasonCol, - seeingEffCol=param.seeingEffCol, seeingGeomCol=param.seeingGeomCol, - airmassCol=param.airmassCol, skyCol=param.skyCol, moonCol=param.moonCol, - salt2Dir=param.salt2Dir) - + nightCol=param.nightCol, + m5Col=param.m5Col, + seasonCol=param.seasonCol, + seeingEffCol=param.seeingEffCol, + seeingGeomCol=param.seeingGeomCol, + airmassCol=param.airmassCol, + skyCol=param.skyCol, moonCol=param.moonCol, + salt2Dir=param.salt2Dir, + psf_flux=param.psf_flux, + frac_flux_seeing=param.frac_flux_seeing, + ccd_full_well=param.ccd_full_well) """ SN class - inherits from SN_Object Parameters @@ -40,31 +51,283 @@ def __init__(self, param, simu_param, reference_lc=None, gamma=None, mag_to_flux version: version of the model (str) reference_lc : griddata,opt reference_light curves (default: None) - gamma: griddata, opt - reference gamma values (default: None) - mag_to_flux: griddata, opt - reference mag->flux values (default: None) - snr_fluxsec: str, opt - type of method to estimate snr and flux in pe.s-1: - lsstsim: estimated from lsstsims tools - interp: estimated from interpolation (default) - all : estimated from the two above-mentioned methods """ + + # self.sncosmo = sncosmo_emul + + self.error_model = self.simulator_parameters['errorModel'] + # self.error_model_cut = self.simulator_parameters['errorModelCut'] + + # dust map + self.dustmap = sncosmo.OD94Dust() + # self.dustmap = sncosmo.CCM89Dust() + + # load info for sigma_mb shift + self.nsigmamb = self.sn_parameters['modelPar']['mbsigma'] + + if np.abs(self.nsigmamb) > 1.e-5: + df = pd.read_hdf(self.sn_parameters['modelPar']['mbsigmafile']) + self.sigma_mb_z = interp1d( + df['z'], df['sigma_mb'], bounds_error=False, fill_value=0.) + model = simu_param['model'] version = str(simu_param['version']) + + self.sn_model = model + self.sn_version = version + + self.dL = self.cosmology.luminosity_distance( + self.sn_parameters['z']).value*1.e3 # in kpc + + self.sn_type = self.sn_parameters['type'] + self.sn_smearFlux = self.sn_parameters['smearFlux'] + + if self.sn_type == 'SN_Ia': + self.source(model, version) + else: + self.random_source(self.sn_type, model, version) + + if self.sn_type == 'SN_Ia': + self.SN.set(t0=self.sn_parameters['daymax'] + + self.gen_parameters['epsilon_daymax']) + if 'salt' in model: + self.SN_SALT(model) + else: + # self.SN.set(t0=self.sn_parameters['daymax']) + # self.SN.set(t0=-36.1+self.sn_parameters['daymax']) + self.SN.set(t0=-37.) + # self.SN.set(t0=-36.06+2.964+self.sn_parameters['daymax']) + + self.defname = dict(zip(['healpixID', 'pixRA', 'pixDec'], [ + 'observationId', param.RACol, param.DecCol])) + + # names for metadata + self.names_meta = ['RA', 'Dec', 'sn_type', 'sn_model', 'sn_version', + 'daymax', + 'z', 'survey_area', + 'healpixID', 'pixRA', 'pixDec', + 'season', 'season_length', 'dL', 'ptime', + 'status', 'ebvofMW', 'lsst_start', 'mjd_max', + 'psf_flux', 'ccdfullwell', 'zmeas'] + + if self.sn_type == 'SN_Ia': + self.names_meta += ['x0', 'epsilon_x0', 'x1', + 'epsilon_x1', 'color', 'epsilon_color', + 'epsilon_daymax'] + """ + self.names_meta += ['minRFphase', 'maxRFphase', + 'minRFphaseQual', 'maxRFphaseQual', 'weight'] + """ + self.names_meta += ['minRFphase', 'maxRFphase', + 'minRFphaseQual', 'maxRFphaseQual'] + self.mag_inf = 100. # mag values to replace infs + + """ + bands = self.zp_airmass['band'].tolist() + slope = self.zp_airmass['slope'].tolist() + intercept = self.zp_airmass['intercept'].tolist() + self.zp_slope = dict(zip(bands, slope)) + self.zp_intercept = dict(zip(bands, intercept)) + """ + # self.register_bands() + + """ + for band in 'grizy': + name = 'LSST::'+band + throughput = self.telescope.atmosphere[band] + try: + band = sncosmo.get_bandpass(name) + except Exception as err: + print(err) + bandcosmo = sncosmo.Bandpass( + throughput.wavelen, throughput.sb, name=name, + wave_unit=u.nm) + sncosmo.registry.register(bandcosmo) + """ + # get sigma_z + z = self.sn_parameters['z'] + sigmaz = self.sn_parameters['sigmaz'] + from random import gauss + self.zmeas = z+gauss(0, sigmaz*(1.+z)) + + # dustmap + self.dust_map = dust_map + + def source(self, model, version): + """ + method to instantiate a source from sncosmo + + Parameters + -------------- + model: str + built-in from sncosmo + version: str + version number + + """ + source = sncosmo.get_source(model, version) + + if model == 'salt3': + source._wave[0] = 1500. # used to be 1700 + source._wave[-1] = 24990. + + self.SN = sncosmo.Model(source=source, + effects=[self.dustmap, self.dustmap], + effect_names=['host', 'mw'], + effect_frames=['rest', 'obs']) + self.model = model + # set cosmology here + self.SN.cosmo = self.cosmology self.version = version - self.gamma = gamma - self.mag_to_flux = mag_to_flux - self.snr_fluxsec = snr_fluxsec - self.error_model = error_model + # self.X0 = self.x0(self.dL) + # self.SN.set(x0=self.X0) + self.SN.set(z=self.sn_parameters['z']) + # print(self.SN) + """ + self.dL = self.cosmology.luminosity_distance( + self.sn_parameters['z']).value*1.e3 + """ + + def random_source(self, sn_type, sn_model='random', sn_version=''): + """ + Method to choose a random source depending on the type + This will occur for non-Ia models + in that case choose return randomly one among available models. + + Parameters + --------------- + sn_type: str + supernovae type + sn_model: str, opt + specific model to run on (default: random) + sn_version: str, opt + specific version to run on (default : ) + + """ + + # load the possible models for this type of supernova + # get the location of the file + location = os.path.realpath( + os.path.join(os.getcwd(), os.path.dirname(__file__))) + + location = location.replace('sn_simulator', 'sn_simu_input') + + # load parameters + df = pd.read_csv( + '{}/sncosmo_builtins.txt'.format(location), delimiter=' ') + + # sn_model!='random': choose this model + df['version'] = df['version'].astype(str) + # print(df.dtypes) + if sn_model != 'random': + idx = df['name'] == sn_model + idx &= df['version'] == str(sn_version) + selm = df[idx] + sn_type = '{}_{}'.format( + selm['type'].values.item(), selm['subtype'].values.item()) + sn_version = str(selm['version'].values.item()) + else: + sn_type, sn_model, sn_version = self.get_sn_fromlist(sn_type, df) + + if sn_type == 'SN_Ia': + sn_type += 'T' + self.sn_type = sn_type + self.sn_model = sn_model + self.sn_version = sn_version + + # self.source(self.sn_model, self.sn_version) + source = sncosmo.get_source(sn_model, sn_version) + self.SN = sncosmo.Model(source=source, + effects=[self.dustmap, self.dustmap], + effect_names=['host', 'mw'], + effect_frames=['rest', 'obs']) + self.SN.set(z=self.sn_parameters['z']) + + """ + self.SN.set_source_peakabsmag(self.sn_parameters['absmag'], + self.sn_parameters['band'], + self.sn_parameters['magsys']) + """ + # print(self.SN) + # self.SN.set(amplitude=2.e-8) + + def get_sn_fromlist(self, sn_type, df): + """ + Method to get a sn model, type, version from a list + + Parameters + -------------- + sn_type: str + type of sn + df: pandas df + list of sn + + Returns + ---------- + sn_type, sn_model, sn_version + """ + + tt = sn_type.split('_')[0] + st = sn_type.split('_')[1] + + if sn_type == 'Non_Ia': + idx = df['type'] == 'SN' + idx &= df['subtype'] != 'Ia' + else: + if sn_type == 'SN_IaT': + idx = df['type'] == tt + idx &= df['subtype'] == 'Ia' + else: + idx = df['type'] == tt + idx &= df['subtype'] == st + + sel = df[idx] + sela = sel.groupby(['type', 'subtype']).size().to_frame( + 'size').reset_index() + io = np.random.choice( + len(sela), 1, p=sela['size'].values/sela['size'].sum())[0] + sn_type = '{}_{}'.format( + sela.iloc[io]['type'], sela.iloc[io]['subtype']) + + main_type = sn_type.split('_')[0] + sub_type = sn_type.split('_')[1] + + idx = df['type'] == main_type + idx &= df['subtype'] == sub_type + idx &= df['name'] != 'snana-04d4jv' # this crashes + + sel = df[idx] + + # print(sel) + # take a random + io = np.random.randint(0, len(sel), 1)[0] + + return sn_type, sel.iloc[io]['name'], str(sel.iloc[io]['version']) + + def SN_SALT(self, model): + """ + Method to set SALT parameters for SN + + Parameters + -------------- + model: str + model name + + """ if model == 'salt2-extended': model_min = 300. model_max = 180000. - wave_min = 3000. - wave_max = 11501. + wave_min = 2000. + wave_max = 11000. + + if model == 'salt3': + model_min = 300. + model_max = 180000. + wave_min = model_min + wave_max = model_max if model == 'salt2': model_min = 3400. @@ -73,69 +336,82 @@ def __init__(self, param, simu_param, reference_lc=None, gamma=None, mag_to_flux wave_max = model_max self.wave = np.arange(wave_min, wave_max, 1.) + self.wave *= (1.+self.sn_parameters['z']) + """ if not self.error_model: source = sncosmo.get_source(model, version=version) else: SALT2Dir = 'SALT2.Guy10_UV2IR' - self.SALT2Templates(SALT2Dir=SALT2Dir, blue_cutoff=10.*self.sn_parameters['blue_cutoff']) + self.SALT2Templates( + SALT2Dir=SALT2Dir, + blue_cutoff=10.*self.sn_parameters['blue_cutoff']) source = sncosmo.SALT2Source(modeldir=SALT2Dir) - - - self.dustmap = sncosmo.OD94Dust() - - self.lsstmwebv = EBV.EBVbase() - - self.SN = sncosmo.Model(source=source, - effects=[self.dustmap, self.dustmap], - effect_names=['host', 'mw'], - effect_frames=['rest', 'obs']) - - self.SN.set(z=self.sn_parameters['z']) - self.SN.set(t0=self.sn_parameters['daymax'] + - self.gen_parameters['epsilon_daymax']) + """ + # set x1 and c parameters self.SN.set(c=self.sn_parameters['color'] + self.gen_parameters['epsilon_color']) self.SN.set(x1=self.sn_parameters['x1'] + self.gen_parameters['epsilon_x1']) - # need to correct X0 for alpha and beta - lumidist = self.cosmology.luminosity_distance( - self.sn_parameters['z']).value*1.e3 - - self.X0 = self.x0(lumidist) - self.dL = lumidist + + self.SN.set_source_peakabsmag(self.sn_parameters['absmag'], + self.sn_parameters['band'], + self.sn_parameters['magsys']) + + if self.sn_parameters['x0']['griddata']: + self.X0 = self.get_x0_from_grid(self.dL) + else: + self.X0 = self.get_x0() + + # set X0 self.SN.set(x0=self.X0) + + # print('after',self.SN.get('x0'),self.SN.get('x1'),self.SN.get('c')) + + def get_x0(self): """ - self.SN.set_source_peakabsmag(self.sn_parameters['absmag'], - self.sn_parameters['band'], self.sn_parameters['magsys']) - self.X0 = self.SN.get('x0') + Method to estimate x0 from (alpha, beta,sigmaint) + + Returns + ------- + X0 : TYPE + DESCRIPTION. + """ - self.defname = dict(zip(['healpixID', 'pixRA', 'pixDec'], [ - 'observationId', param.RACol, param.DecCol])) + X0 = self.SN.get('x0') + alpha = self.sn_parameters['alpha'] + beta = self.sn_parameters['beta'] - # names for metadata - self.names_meta = ['RA', 'Dec', - 'x0', 'epsilon_x0', - 'x1', 'epsilon_x1', - 'color', 'epsilon_color', - 'daymax', 'epsilon_daymax', - 'z', 'survey_area', - 'healpixID', 'pixRA', 'pixDec', - 'season', 'dL', 'ptime', 'snr_fluxsec_meth', 'status', 'ebvofMW'] + X0 *= np.power(10., 0.4*(alpha * + self.sn_parameters['x1'] - beta * + self.sn_parameters['color'])) - self.mag_inf = 100. # mag values to replace infs + if self.sn_parameters['sigmaInt'] > 0 or np.abs(self.nsigmamb) > 1.e-5: - # band registery in sncosmo - for band in 'grizy': - throughput = self.telescope.atmosphere[band] - bandcosmo = sncosmo.Bandpass( - throughput.wavelen, throughput.sb, name='LSST::'+band, wave_unit=u.nm) - sncosmo.registry.register(bandcosmo, force=True) + # estimate mb + mb = -2.5*np.log10(X0)+10.635 + if self.sn_parameters['sigmaInt'] > 0: + # smear it + from random import gauss - def x0(self, lumidist): + mb = gauss(mb, self.sn_parameters['sigmaInt']) + + if np.abs(self.nsigmamb) > 1.e-5: + # get sigma from z + sigmamb = self.sigma_mb_z(self.sn_parameters['z']) + mb += self.nsigmamb*sigmamb + + # and recalculate X0 + X0 = 10**(-0.4*(mb-10.635)) + + X0 += self.gen_parameters['epsilon_x0'] + + return X0 + + def get_x0_from_grid(self, lumidist): """" Method to estimate x0 from a griddata @@ -145,45 +421,71 @@ def x0(self, lumidist): luminosity distance """ + from scipy.interpolate import griddata + + X0_grid = griddata((self.x0_grid['x1'], self.x0_grid['color']), + self.x0_grid['x0_norm'], ( + self.sn_parameters['x1'], self.sn_parameters['color']), + method='nearest') - X0_grid = griddata((self.x0_grid['x1'], self.x0_grid['color']), self.x0_grid['x0_norm'], ( - self.sn_parameters['x1'], self.sn_parameters['color']), method='nearest') X0 = X0_grid / lumidist ** 2 - alpha = 0.13 - beta = 3. + + alpha = self.sn_parameters['alpha'] + beta = self.sn_parameters['beta'] + X0 *= np.power(10., 0.4*(alpha * self.sn_parameters['x1'] - beta * self.sn_parameters['color'])) + + if self.sn_parameters['sigmaInt'] > 0 or np.abs(self.nsigmamb) > 1.e-5: + + # estimate mb + mb = -2.5*np.log10(X0)+10.635 + + if self.sn_parameters['sigmaInt'] > 0: + # smear it + from random import gauss + + mb = gauss(mb, self.sn_parameters['sigmaInt']) + + if np.abs(self.nsigmamb) > 1.e-5: + # get sigma from z + sigmamb = self.sigma_mb_z(self.sn_parameters['z']) + mb += self.nsigmamb*sigmamb + + # and recalculate X0 + X0 = 10**(-0.4*(mb-10.635)) + X0 += self.gen_parameters['epsilon_x0'] - - return X0 + return X0 - def SALT2Templates(self,SALT2Dir='SALT2.Guy10_UV2IR', blue_cutoff=3800.): + def SALT2Templates(self, SALT2Dir='SALT2.Guy10_UV2IR', blue_cutoff=3800.): """ Method to load SALT2 templates and apply cutoff on SED. - + Parameters -------------- SALT2Dir: str, opt - SALT2 directory (default: SALT2.Guy10_UV2IR) + SALT2 directory(default: SALT2.Guy10_UV2IR) blue_cutoff: float, opt - blue cut off to apply (in nm - default: 3800.) + blue cut off to apply(in nm - default: 3800.) """ for vv in ['salt2_template_0', 'salt2_template_1']: fName = '{}/{}_orig.dat'.format(SALT2Dir, vv) - data = np.loadtxt(fName, dtype={'names': ('phase', 'wavelength', 'flux'), - 'formats': ('f8', 'i4', 'f8')}) - #print(data) + data = np.loadtxt(fName, dtype={'names': ('phase', 'wavelength', + 'flux'), + 'formats': ('f8', 'i4', 'f8')}) + # print(data) data['flux'][data['wavelength'] <= blue_cutoff] = 0.0 - #print(data) + # print(data) np.savetxt('{}/{}.dat'.format(SALT2Dir, vv), data, fmt=['%1.2f', '%4d', '%.7e', ]) - - def __call__(self, obs, display=False, time_display=0.): + + def __call__(self, obs, display=False, time_display=0., lc_coadd=0): """ Simulation of the light curve Parameters @@ -196,6 +498,8 @@ def __call__(self, obs, display=False, time_display=0.): time_display: float duration(sec) for which the display is visible default: 0 + lc_coadd: int, opt + to coadd lc fluxes Returns ----------- @@ -223,7 +527,8 @@ def __call__(self, obs, display=False, time_display=0.): flux: SN flux(Jy) fluxerr: EN error flux(Jy) snr_m5: Signal-to-Noise Ratio(float) - gamma: gamma parameter(see LSST: From Science...data products eq. 5)(float) + gamma: gamma parameter(see LSST: + From Science...data products eq. 5)(float) m5: five-sigma depth(float) seeingFwhmEff: seeing eff(float) seeingFwhmGeom: seeing geom(float) @@ -237,10 +542,25 @@ def __call__(self, obs, display=False, time_display=0.): time: time(days)(float) phase: phase(float) """ + + # start timer + ti = SNTimer(time.time()) + + if len(obs) == 0: + pix = {} + for vv in ['healpixID', 'pixRA', 'pixDec']: + pix[vv] = -1 + return [self.nosim(-1., -1., pix, -1., -1., -1., + ti, -1, -1., -1., -1., + self.psf_flux, self.ccd_full_well, -1.)] + ra = np.mean(obs[self.RACol]) dec = np.mean(obs[self.DecCol]) area = self.area + season = np.unique(obs['season'])[0] + season_length = np.max(obs[self.mjdCol])-np.min(obs[self.mjdCol]) + pix = {} for vv in ['healpixID', 'pixRA', 'pixDec']: if vv in obs.dtype.names: @@ -249,163 +569,575 @@ def __call__(self, obs, display=False, time_display=0.): pix[vv] = np.mean(obs[self.defname[vv]]) ebvofMW = self.sn_parameters['ebvofMW'] - # apply dust here since Ra, Dec is known - if ebvofMW < 0.: - ebvofMW = self.lsstmwebv.calculateEbv( - equatorialCoordinates=np.array( - [[ra], [dec]]))[0] + if ebvofMW < 0: + idx = self.dust_map['healpixID'] == pix['healpixID'] + sel = self.dust_map[idx] + ebvofMW = sel['ebvofMW'].mean() self.SN.set(mwebv=ebvofMW) - # start timer - ti = SNTimer(time.time()) + # add atmospheric parameters here + # obs = self.set_atmos_params(obs) - # Are there observations with the filters? - goodFilters = np.in1d(obs[self.filterCol], - np.array([b for b in 'grizy'])) + # register throughputs for sn_cosmo - if len(obs[goodFilters]) == 0: - return [self.nosim(ra, dec, pix, area, season, ti, self.snr_fluxsec, -1, ebvofMW)] + # obs = self.register_bands_from_atmos(obs) - # Select obs depending on min and max phases - # blue and red cutoffs applied + # estimate zp and mean_wavelength corresponding to obs + """ + if self.atmosType == 'const': + obs = self.add_zp_meanwave_from_interp(obs) + else: + obs = self.add_zp_meanwave_from_obs(obs) + """ - if not self.error_model: - obs = self.cutoff(obs, self.sn_parameters['daymax'], - self.sn_parameters['z'], - self.sn_parameters['min_rf_phase'], - self.sn_parameters['max_rf_phase'], - self.sn_parameters['blue_cutoff'], - self.sn_parameters['red_cutoff']) + # filter cutoffs - if len(obs) == 0: - return [self.nosim(ra, dec, pix, area, season, ti, self.snr_fluxsec, -1, ebvofMW)] + obs = self.select_filter_cutoff(obs, + ra, dec, pix, area, season, + season_length, + ti, ebvofMW) - # Sort data according to mjd - obs.sort(order=self.mjdCol) + lsst_start = -1 + mjd_max = -1.0 + if len(obs) == 0: + return [self.nosim(ra, dec, pix, area, season, season_length, + ti, -1, ebvofMW, lsst_start, mjd_max, + self.psf_flux, self.ccd_full_well, -1.)] # preparing the results : stored in lcdf pandas DataFrame - outvals = [self.m5Col, self.mjdCol, - self.exptimeCol, self.nexpCol, self.filterCol] - for bb in [self.airmassCol, self.skyCol, self.moonCol, self.seeingEffCol, self.seeingGeomCol]: - if bb in obs.dtype.names: - outvals.append(bb) + a = [self.m5Col, self.mjdCol, + self.exptimeCol, self.nexpCol, + self.filterCol, self.nightCol] + a += [self.airmassCol, self.skyCol, self.moonCol, + self.seeingEffCol, self.seeingGeomCol, + 'zp', 'mean_wave', 'pwv', 'ozone', + 'aerosol', 'band_cosmo', + 'sigma_zp', 'tel_site_name'] + for vv in ['airmass', 'pwv', 'ozone', 'aerosol']: + for vvb in ['sigma', 'min', 'max', 'round']: + a.append('{}_{}'.format(vvb, vv)) + + b = obs.dtype.names - lcdf = pd.DataFrame(obs[outvals]) + outvals = list(set(a) & set(b)) - # print(self.fluxSED(lcdf)) + lcdf = pd.DataFrame(np.copy(obs[outvals])) - band_cosmo = '{}_cosmo'.format(self.filterCol) - lcdf[band_cosmo] = 'LSST::'+lcdf[self.filterCol] + # smear zp here + lcdf['zp'] += np.random.normal(0, lcdf['sigma_zp']) + lcdf['zpsys'] = 'ab' + + # get band flux + # lcdf = self.get_register_throughput(lcdf) + + # lcdf['zp'] = lcdf['zp_e_sec'] + # get band flux + lcdf['flux'] = self.SN.bandflux( + lcdf['band_cosmo'], lcdf[self.mjdCol], zpsys=lcdf['zpsys'], + zp=lcdf['zp']) + + """ lcdf['flux'] = self.SN.bandflux( + lcdf[band_cosmo], lcdf[self.mjdCol], zpsys='ab', + zp=2.5*np.log10(3631)) + lcdf['zp'] = 2.5*np.log10(3631) + """ + + """ + # flux in JY + lcdf['flux_old'] = self.SN.bandflux( lcdf[band_cosmo], lcdf[self.mjdCol], zpsys='ab', zp=2.5*np.log10(3631)) + lcdf['mag_old'] = -2.5 * np.log10(lcdf['flux_old'] / 3631.0) + """ # estimate error model (if necessary) - if self.error_model: + + lcdf['fluxerr_model'] = 0. + if self.error_model and self.sn_type == 'SN_Ia': fluxcov_cosmo = self.SN.bandfluxcov( - lcdf[band_cosmo], lcdf[self.mjdCol], zpsys='ab', zp=2.5*np.log10(3631)) - lcdf['variance_model'] = np.diag(fluxcov_cosmo[1]) - + lcdf['band_cosmo'], lcdf[self.mjdCol], zpsys='ab', + zp=2.5*np.log10(3631)) + lcdf['fluxerr_model'] = np.sqrt(np.diag(fluxcov_cosmo[1])) + + """ + lcdf.loc[lcdf.flux <= 0., 'fluxerr_photo'] = -1. + lcdf.loc[lcdf.flux <= 0., 'fluxerr_model'] = -1. + lcdf.loc[lcdf.flux <= 0., 'flux'] = 9999. + lcdf.loc[lcdf.flux_old <= 0., 'flux'] = 9999. + """ + # positive flux only + """ idx = lcdf['flux'] > 0. lcdf = lcdf[idx] - #print('simulating',season,len(lcdf)) + """ + if len(lcdf) == 0: - return [] + return [self.nosim(ra, dec, pix, area, season, season_length, + ti, -1, ebvofMW, -1., -1., + self.psf_flux, self.ccd_full_well, -1.)] # ti(time.time(), 'fluxes_b') - # magnitudes - integrated fluxes are in Jy - lcdf['mag'] = -2.5 * np.log10(lcdf['flux'] / 3631.0) + # magnitudes - fluxes are in ADU/s + lcdf['mag'] = -2.5 * np.log10(lcdf['flux'])+lcdf['zp'] # if mag have inf values -> set to 50. lcdf['mag'] = lcdf['mag'].replace([np.inf, -np.inf], self.mag_inf) - # ti(time.time(), 'mags') - - # SNR and flux in pe.sec estimations + # flux error + flux5 = 10**(-0.4*(lcdf[self.m5Col]-lcdf['zp'])) + sigma_5 = flux5/5. + shot_noise = np.sqrt(lcdf['flux']/lcdf[self.exptimeCol]) + # shot_noise = 0.0 + lcdf['fluxerr'] = np.sqrt(sigma_5**2+shot_noise**2) + lcdf['snr_m5'] = lcdf['flux']/lcdf['fluxerr'] + # lcdf['fluxerr'] = lcdf['flux']/lcdf['snr_m5'] - if self.snr_fluxsec == 'all' or self.snr_fluxsec == 'lsstsim': - lcdf = self.calcSNR_Flux(lcdf, transes) - - if self.snr_fluxsec == 'all' or self.snr_fluxsec == 'interp': - gammaName = 'gamma' - fluxName = 'flux_e_sec' - snrName = 'snr_m5' - - if gammaName in lcdf.columns: - gammaName += '_interp' - fluxName += '_interp' - snrName += '_interp' - - - lcdf = lcdf.groupby([self.filterCol]).apply( - lambda x: self.interp_gamma_flux(x, gammaName, fluxName)).reset_index() - - - lcdf[snrName] = 1./srand( - lcdf[gammaName].values, lcdf['mag'], lcdf[self.m5Col]) - - # ti(time.time(), 'estimate 1') - - # complete the LC - lcdf['magerr'] = (2.5/np.log(10.))/lcdf['snr_m5'] # mag error - lcdf['fluxerr'] = lcdf['flux']/lcdf['snr_m5'] # flux error - - if self.error_model: - #print(z,lcdf[['variance_model','fluxerr']]) - lcdf['fluxerr'] = np.sqrt(lcdf['variance_model']+lcdf['fluxerr']**2) #flux error - lcdf['snr_m5'] = lcdf['flux']/lcdf['fluxerr'] # snr - lcdf['magerr'] = (2.5/np.log(10.))/lcdf['snr_m5'] # mag error - - - lcdf['zp'] = 2.5*np.log10(3631) # zp - lcdf['zpsys'] = 'ab' # zpsys + lcdf['magerr_phot'] = (2.5/np.log(10.))/lcdf['snr_m5'] # mag error + lcdf['fluxerr_photo'] = lcdf['fluxerr'] + + lcdf['fluxerr'] = np.sqrt( + lcdf['fluxerr_model']**2+lcdf['fluxerr_photo']**2) # flux error + lcdf['snr'] = lcdf['flux']/lcdf['fluxerr'] # snr + lcdf['magerr'] = (2.5/np.log(10.))/lcdf['snr'] # mag error + # lcdf['zpsys'] = 'ab' # zpsys lcdf['phase'] = (lcdf[self.mjdCol]-self.sn_parameters['daymax'] )/(1.+self.sn_parameters['z']) # phase - - - # rename some of the columns lcdf = lcdf.rename( - columns={self.mjdCol: 'time', self.filterCol: 'band', self.m5Col: 'm5', self.exptimeCol: 'exptime'}) - lcdf['band'] = 'LSST::'+lcdf['band'] + columns={self.mjdCol: 'time', self.filterCol: 'band', + self.m5Col: 'm5', self.exptimeCol: 'exptime'}) - # remove rows with mag_inf values + lcdf['filter'] = lcdf['band'] + lcdf['band'] = lcdf['band_cosmo'] + # remove rows with mag_inf values + """ idf = lcdf['mag'] < self.mag_inf lcdf = lcdf[idf] + """ + + lcdf.loc[lcdf.fluxerr_model < 0, 'flux'] = 0. + lcdf.loc[lcdf.fluxerr_model < 0, 'fluxerr_photo'] = 10. + lcdf.loc[lcdf.fluxerr_model < 0, 'fluxerr'] = 10. + lcdf.loc[lcdf.fluxerr_model < 0, 'snr_m5'] = 0. + lcdf.loc[lcdf.fluxerr_model < 0, 'fluxerr_model'] = 10. + + # smear lc fluxes + if self.sn_smearFlux: + lcdf['flux'] = lcdf['flux']+np.random.normal(0., lcdf['fluxerr']) + lcdf.loc[lcdf.flux < 0, 'flux'] = 0. + lcdf['snr'] = lcdf['flux']/lcdf['fluxerr'] + + # smear atmos parameters + # print('before', lcdf[['airmass', 'pwv', 'ozone', 'aerosol']]) + lcdf = self.smear_atmos(lcdf) + # print('after', lcdf[['airmass', 'pwv', 'ozone', 'aerosol']]) + + """ + filters = np.array(lcdf['filter']) + filters = filters.reshape((len(filters), 1)) + + lcdf['lambdabar'] = np.apply_along_axis(self.lambdabar, 1, filters) + """ if len(lcdf) == 0: - return [self.nosim(ra, dec, pix, area, season, ti, self.snr_fluxsec, -1, ebvofMW)] + return [self.nosim(ra, dec, pix, area, season, season_length, + ti, -1, ebvofMW, -1., -1., + self.psf_flux, self.ccd_full_well, -1.)] # get the processing time ptime = ti.finish(time.time())['ptime'].item() - # transform pandas df to astropy Table + """ + print('kkkkk', lcdf[['mag', 'mag_old', 'flux', + 'flux_old', 'snr_m5', 'm5', 'band', 'filter']], len(lcdf)) + """ + + # include saturation effects here + if self.frac_flux_seeing is not None: + lcdf = self.estimate_satured_flux(lcdf) + else: + lcdf['sat'] = 0 + # transform pandas df to astropy Table + table_lc = Table.from_pandas(lcdf) + # set metadata + if 'lsst_start' in obs.dtype.names: + lsst_start = np.median(obs['lsst_start']) + + mjd_max = -1 + if len(table_lc) > 0: + mjd_max = np.max(table_lc['time']) + table_lc.meta = self.metadata( - ra, dec, pix, area, season, ptime, self.snr_fluxsec, 1, ebvofMW) + ra, dec, pix, area, season, season_length, ptime, + 1, ebvofMW, lsst_start, mjd_max, + self.psf_flux, self.ccd_full_well, self.zmeas) + + if len(table_lc) == 0: + return [table_lc] # if the user chooses to display the results... if display: self.plotLC(table_lc['time', 'band', - 'flux', 'fluxerr', 'zp', 'zpsys'], time_display) + 'flux', 'fluxerr', 'zp', 'zpsys'], + time_display) + + # remove LC points with too high error model value + """ + if self.error_model: + if self.error_model_cut >0: + idx = table_lc['fluxerr_model']/ \ + table_lc['flux']<= self.error_model_cut + table_lc = table_lc[idx] + """ + toremove = ['m5', 'exptime', 'numExposures', 'filter_cosmo', + 'airmass', 'moonPhase', + 'seeingFwhmEff', 'seeingFwhmGeom', 'gamma', 'mag', + 'magerr', 'flux_e_sec', 'magerr_phot'] + + toremove = ['m5', 'exptime', 'numExposures', 'filter_cosmo', + 'airmass', 'moonPhase', + 'seeingFwhmEff', 'seeingFwhmGeom', 'gamma', 'mag', + 'magerr', 'magerr_phot'] + + toremove = ['filter_cosmo', 'airmass', 'moonPhase', + 'gamma', 'mag', 'magerr', 'magerr_phot'] + + # table_lc.remove_columns(toremove) + + if lc_coadd: + from sn_tools.sn_lcana import coadd_lc + table_lc = coadd_lc(table_lc) return [table_lc] - def calcSNR_Flux(self, df, transm): + def smear_atmos(self, lcdfa): + """ + Method to smear atmospheric parameters + + Parameters + ---------- + lcdf : pandas df + data to smear. + + Returns + ------- + lcdf: pandas df + date with smeared atmos params + + """ + + ccols = ['airmass', 'ozone', 'pwv', 'aerosol'] + + lcdf = pd.DataFrame(lcdfa) + for vv in ccols: + lcdf[vv] += np.random.normal(0., lcdf['sigma_{}'.format(vv)]) + + # check whether values are inside the allowed range (getobsatmo) + + for vv in ccols: + idx = lcdf[vv] < lcdf['min_{}'.format(vv)] + if len(lcdf[idx]) > 0: + lcdf.loc[idx, vv] = lcdf['min_{}'.format(vv)] + idx = lcdf[vv] > lcdf['max_{}'.format(vv)] + if len(lcdf[idx]) > 0: + lcdf.loc[idx, vv] = lcdf['max_{}'.format(vv)] + + res = pd.DataFrame(lcdf) + + # round atmos params """ - Method to estimate SNRs and fluxes (in e.sec) + for vv in ccols: + # round_value = eval('self.{}_round'.format(vv)) + round_value = int(res['round_{}'.format(vv)].mean()) + res[vv] = np.round(res[vv], round_value) + """ + del lcdf + return res + + def register_bands_from_atmos_deprecated(self, obs, + cols=['airmass', + 'pwv', 'ozone', 'aerosol']): + """ + Method to register throughputs from atmos params + + Parameters + ---------- + obs : numpy array + observations. + cols : list(str), optional + columns of interest. The default is ['airmass','pwv', 'ozone', 'aerosol']. + + Returns + ------- + obs : numpy array + observations. + + """ + import time + time_ref = time.time() + # round atmos parameters + for vv in cols: + obs[vv] = np.round(obs[vv], 2) + + # set band_cosmo column + bcols = self.telescope.site_name+'::' + \ + obs[self.filterCol]+'_' + \ + obs[self.airmassCol].astype(str)+'_' + \ + obs['pwv'].astype(str)+'_' + \ + obs['ozone'].astype(str)+'_' + \ + obs['aerosol'].astype(str) + + import numpy.lib.recfunctions as rf + obs = rf.append_fields(obs, 'band_cosmo', bcols.tolist()) + + # grab unique parameters + all_cols = [self.filterCol]+['band_cosmo']+cols + atmos_table = unique(Table(obs[all_cols])) + + # register + + self.register_bands_on_the_fly(self.telescope, + atmos_table[all_cols].to_pandas()) + + return obs + + def get_zp_deprecated(self, lcdf): + """ + Method to estimate zero points + + Parameters + ---------- + lcdf : TYPE + DESCRIPTION. + + Returns + ------- + lcdf : TYPE + DESCRIPTION. + + """ + + ra = [] + rb = [] + from random import gauss + for i, row in lcdf.iterrows(): + airmass = row['airmass'] + pwv = row['pwv'] + ozone = row['ozone'] + aerosol = row['aerosol'] + b = row['filter'] + pwv += gauss(0.2) + # ozone += gauss(5.) + self.telescope.new_atmosphere(site_name=self.telescope.site_name, + airmass=airmass, + aerosol=aerosol, + pwv=pwv, ozone=ozone) + self.telescope.reset_data() + self.telescope.mean_wave() + # grab zp + mean_wave = self.telescope.mean_wavelength[b] + zp = self.telescope.zp(b) + ra.append((zp)) + rb.append((mean_wave)) + + lcdf['zp_new'] = ra + lcdf['mean_wave'] = rb + + return lcdf + + def set_atmos_params_deprecated(self, lcdf): + """ + Method to set atmospheric parameters + + Parameters + ---------- + lcdf : pandas df + Data to process. + + Returns + ------- + lcdf : pandas df + output data. + + """ + + if 'pwv' not in lcdf.columns: + lcdf['pwv'] = self.pwv + if 'ozone' not in lcdf.columns: + lcdf['ozone'] = self.ozone + if 'aerosol' not in lcdf.columns: + lcdf['aerosol'] = self.aerosol + + lcdf = lcdf.round( + {self.airmassCol: 2, 'pwv': 2, 'ozone': 2, 'aerosol': 2}) + + return lcdf + + def get_register_throughput(self, lcdf): + """ + Method to estimate throughput and register in sncosmo + + Parameters + ---------- + lcdf : pandas df + Data to process. + + Returns + ------- + lcdf : pandas df + init df plus additionnal params. + + """ + + bandname = self.telescope.site_name+'::' + \ + lcdf[self.filterCol]+'_' + \ + lcdf[self.airmassCol].astype(str)+'_' + \ + lcdf['pwv'].astype(str)+'_' + \ + lcdf['ozone'].astype(str)+'_' + \ + lcdf['aerosol'].astype(str) + + lcdf['band_cosmo'] = bandname + + self.register_bands_on_the_fly(self.telescope, + lcdf[['band_cosmo', + self.filterCol, + self.airmassCol, + 'pwv', 'ozone', 'aerosol']]) + + return lcdf + + def estimate_satured_flux(self, lcdf): + """ + Method to remove saturating LC points + + Parameters + ---------- + lcdf : pandas df + Data to process. + + Returns + ------- + lcdf : pandas df + Original data with saturated LC points removed.. + + """ + + lcdf['sat'] = 0 + vvar = 'single_exposure_flux' + lcdf[vvar] = lcdf['flux'] * \ + self.frac_flux_seeing(lcdf['seeingFwhmEff']) + lcdf[vvar] *= lcdf['exptime']/lcdf['numExposures'] + lcdf[vvar] -= self.ccd_full_well + idx = lcdf[vvar] >= 0 + lcdf.loc[idx, 'sat'] = 1 + + # lcdf = pd.DataFrame(lcdf[idx]) + lcdf = lcdf.drop(columns=[vvar]) + return lcdf + + def select_filter_cutoff(self, obs, + ra, dec, pix, area, season, + season_length, + ti, ebvofMW): + """ + Function to select obs: filters and blue and red cutoffs + + Parameters + ---------- + obs : array + data to process. + ra : float + pixRA. + dec : float + pixDec. + pix : TYPE + DESCRIPTION. + area : str + pixarea. + season : int + season of obs. + season_length : float + season length + ti : float + timer val. + ebvofMW : float + E(B-V) MW + + Returns + ------- + array + selected data. + + """ + + # Are there observations with the filters? + goodFilters = np.in1d(obs[self.filterCol], + np.array([b for b in 'grizy'])) + + if len(obs[goodFilters]) == 0: + return [self.nosim(ra, dec, pix, area, season, season_length, + ti, -1, ebvofMW, -1.0, -1.0, + self.psf_flux, self.ccd_full_well, -1.)] + + # Select obs depending on min and max phases + # blue and red cutoffs applied + + blue_cutoffs, red_cutoffs = self.get_cutoffs() + + obs = self.cutoff(obs, self.sn_parameters['daymax'], + self.sn_parameters['z'], + self.sn_parameters['minRFphase'], + self.sn_parameters['maxRFphase'], + blue_cutoffs, red_cutoffs) + if len(obs) > 0: + # Sort data according to mjd + obs.sort(order=self.mjdCol) + + return obs + + def get_cutoffs(self): + """ + Method to estimate blue and red cutoffs as dict + + Returns + ------- + blue_cutoffs : dict + blue cutoffs (key=band). + red_cutoffs : dict + red cutoffs (key=band) + + """ + + blue_cutoff = 0. + red_cutoff = 1.e8 + blue_cutoffs = dict(zip('urgizy', [blue_cutoff]*6)) + red_cutoffs = dict(zip('urgizy', [red_cutoff]*6)) + if self.sn_type == 'SN_Ia' and not self.error_model: + for b in 'ugrizy': + blue_cutoffs[b] = self.sn_parameters['blueCutoff{}'.format( + b)] + red_cutoffs[b] = self.sn_parameters['redCutoff{}'.format( + b)] + + return blue_cutoffs, red_cutoffs + + def calcSNR_Flux_deprecated(self, df, transm): + """ + Method to estimate SNRs and fluxes(in e.sec) using lsst sims estimators Parameters --------------- df: pandas df data to process - transm : array + transm: array throughputs Returns @@ -419,11 +1151,14 @@ def calcSNR_Flux(self, df, transm): # estimate SNR # Get photometric parameters to estimate SNR - photParams = [PhotometricParameters(exptime=vv[self.exptimeCol]/vv[self.nexpCol], - nexp=vv[self.nexpCol]) for index, vv in df.iterrows()] + from rubin_sim.phot_utils.photometric_parameters import PhotometricParameters + from rubin_sim.phot_utils import signaltonoise + photParams = [PhotometricParameters( + exptime=vv[self.exptimeCol]/vv[self.nexpCol], + nexp=vv[self.nexpCol]) for index, vv in df.iterrows()] nvals = range(len(df)) - calc = [SignalToNoise.calcSNR_m5( + calc = [signaltonoise.calcSNR_m5( df.iloc[i]['mag'], transm[i], df.iloc[i][self.m5Col], photParams[i]) for i in nvals] @@ -431,11 +1166,13 @@ def calcSNR_Flux(self, df, transm): df['gamma'] = [calc[i][1] for i in nvals] # estimate the flux in elec.sec-1 df['flux_e_sec'] = self.telescope.mag_to_flux_e_sec( - df['mag'].values, df[self.filterCol].values, df[self.exptimeCol]/df[self.nexpCol], df[self.nexpCol])[:, 1] + df['mag'].values, df[self.filterCol].values, + df[self.exptimeCol]/df[self.nexpCol], df[self.nexpCol])[:, 1] return df - def interp_gamma_flux(self, grp, gammaName='gamma_int', fluxName='flux_e_sec_int'): + def interp_gamma_flux_deprecated(self, grp, gammaName='gamma_int', + fluxName='flux_e_sec_int'): """ Method to estimate gamma and mag_to_flux values from interpolation @@ -461,9 +1198,18 @@ def interp_gamma_flux(self, grp, gammaName='gamma_int', fluxName='flux_e_sec_int grp.loc[:, fluxName] = self.mag_to_flux[grp.name]( (grp['mag'], single_exptime, grp[self.nexpCol])) + """ + # SNR_photo_bd + grp.loc[:, 'f5'] = self.mag_to_flux[grp.name]( + (grp['fiveSigmaDepth'], single_exptime, grp[self.nexpCol])) + + grp.loc[:, 'SNR_photo_bd'] = grp[fluxName]/(grp['f5']/5.) + """ return grp - def nosim(self, ra, dec, pix, area, season, ti, snr_fluxsec, status, ebvofMW): + def nosim(self, ra, dec, pix, area, season, season_length, + ti, status, ebvofMW, lsst_start=-1.0, mjd_max=-1, + psf_flux='', ccd_full_well=133000, zmeas=-1.): """ Method to construct an empty table when no simulation was not possible @@ -479,23 +1225,26 @@ def nosim(self, ra, dec, pix, area, season, ti, snr_fluxsec, status, ebvofMW): survey area season: int season of interest + season_length: float + season length ptime: float processing time - snr_fluxsec: str - method used to estimate snr and flux in pe.s-1 status: int - status of the processing (1=ok, -1=no simu) - ebvofMW : float + status of the processing(1=ok, -1=no simu) + ebvofMW: float E(B-V) of MW """ ptime = ti.finish(time.time())['ptime'].item() table_lc = Table() # set metadata table_lc.meta = self.metadata( - ra, dec, pix, area, season, ptime, snr_fluxsec, status, ebvofMW) + ra, dec, pix, area, season, season_length, ptime, + status, ebvofMW, lsst_start, mjd_max, psf_flux, ccd_full_well, zmeas) return table_lc - def metadata(self, ra, dec, pix, area, season, ptime, snr_fluxsec, status, ebvofMW): + def metadata(self, ra, dec, pix, area, season, season_length, + ptime, status, ebvofMW, lsst_start, mjd_max, psf_flux='', + ccd_full_well=133000, zmeas=-1.): """ Method to fill metadata @@ -506,17 +1255,17 @@ def metadata(self, ra, dec, pix, area, season, ptime, snr_fluxsec, status, ebvof dec: float SN dec pix: dict - pixel infos (ID, RA, Dec) + pixel infos(ID, RA, Dec) area: float area of the survey season: float season number + season_length: float + season length ptime: float processing time - snr_fluxsec: str - method used to estimate snr and flux in pe.s-1 status: int - status of the simulation (1=ok, -1=not simulation) + status of the simulation(1=ok, -1=not simulation) Returns ----------- @@ -524,20 +1273,120 @@ def metadata(self, ra, dec, pix, area, season, ptime, snr_fluxsec, status, ebvof """ - val_meta = [ra, dec, - self.X0, self.gen_parameters['epsilon_x0'], - self.sn_parameters['x1'], self.gen_parameters['epsilon_x1'], - self.sn_parameters['color'], self.gen_parameters['epsilon_color'], - self.sn_parameters['daymax'], self.gen_parameters['epsilon_daymax'], + val_meta = [ra, dec, self.sn_type, self.sn_model, self.sn_version, + self.sn_parameters['daymax'], self.sn_parameters['z'], area, pix['healpixID'], pix['pixRA'], pix['pixDec'], - season, self.dL, ptime, snr_fluxsec, status, ebvofMW] + season, season_length, self.dL, ptime, + status, ebvofMW, lsst_start, mjd_max, + psf_flux, ccd_full_well, zmeas] + + if self.sn_type == 'SN_Ia': + val_meta += [self.X0, self.gen_parameters['epsilon_x0'], + self.sn_parameters['x1'], + self.gen_parameters['epsilon_x1'], + self.sn_parameters['color'], + self.gen_parameters['epsilon_color'], + self.gen_parameters['epsilon_daymax'], + self.gen_parameters['minRFphase'], + self.gen_parameters['maxRFphase'], + self.gen_parameters['minRFphaseQual'], + self.gen_parameters['maxRFphaseQual']] + # self.gen_parameters['weight']] return dict(zip(self.names_meta, val_meta)) + def SN_SED(self, gen_params, nspectra=3): + """ + Method to generate SN SEDs + + Parameters + --------- + gen_params : list + Simulation parameters. + + Returns + ------- + sed : astropy table + Generated SEDs + + """ + + daymax = gen_params['daymax'] + + min_mjd = daymax-10 + max_mjd = daymax+10 + + mjds = np.random.uniform(low=min_mjd, high=max_mjd, size=nspectra) + + sed = Table() + for io, mjd in enumerate(mjds): + diff = mjd-daymax + sedm = self.SN_SED_mjd(mjd) + sedm['spec'] = 'spec_{}'.format(io) + # sedm['time'] = int(diff) + sedm['mjd'] = mjd + sedm['exptime'] = 0.0 + sedm['valid'] = 1 + + sed = vstack([sed, Table(sedm)]) + + # self.plot_SED(sed) + + return [sed] + + def plot_SED(self, sed): + """ + Method to plot generated SEDs + + Parameters + ---------- + sed : astropy table + generated seds. + + Returns + ------- + None. + + """ + + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + # this is to plot SEDs + for sedid in np.unique(sed['SED_id']): + idx = sed['SED_id'] == sedid + sel = sed[idx] + ax.plot(sel['wave'], sel['flux'], marker='.', + label='{}'.format(np.mean(sel['time']))) + ax.legend() + plt.show() + + def SN_SED_mjd(self, mjd): + """ + Method to generate SED flux + + Parameters + ---------- + mjd : float + MJD for the flux generation. + + Returns + ------- + sed : astropy table + generated fluxes. + + """ + + fluxes = 10.*self.SN.flux(mjd, self.wave) + sed = Table([fluxes], names=['flux']) + sed['wavelength'] = self.wave + sed['fluxerr'] = 0.0 + + return sed + def fluxSED(self, obs): """ - Method to estimate the fluxes (in Jy) using + Method to estimate the fluxes(in Jy) using integration of the flux over bandwidth Parameters @@ -551,7 +1400,7 @@ def fluxSED(self, obs): the fluxes """ - + from rubin_sim.phot_utils import Sed # Get the fluxes (vs wavelength) for each obs fluxes = 10.*self.SN.flux(obs[self.mjdCol], self.wave) @@ -573,7 +1422,49 @@ def fluxSED(self, obs): int_fluxes = np.asarray( [seds[i].calcFlux(bandpass=transes[i]) for i in nvals]) - # negative fluxes are a pb for mag estimate -> set neg flux to nearly nothing + # negative fluxes are a pb for mag estimate + # -> set neg flux to nearly nothing int_fluxes[int_fluxes < 0.] = 1.e-10 return int_fluxes + + def lambdabar_deprecated(self, band): + + return self.mean_wavelength[band[0]] + + def ebvofMW_calc_deprecated(self, healpixID): + """ + Method to estimate E(B-V) + + Parameters + -------------- + RA: float + RA coord. + Dec: float + Dec coord + + Returns + ---------- + E(B-V) + + """ + + # in that case ebvofMW value is taken from a map + from astropy.coordinates import SkyCoord + + coords = SkyCoord(RA, Dec, unit='deg') + + """ + try: + sfd = SFDQuery() + except Exception: + from dustmaps.config import config + config['data_dir'] = 'dustmaps' + import dustmaps.sfd + dustmaps.sfd.fetch() + # dustmaps('dustmaps') + """ + sfd = SFDQuery() + ebvofMW = sfd(coords) + + return ebvofMW diff --git a/sn_simulator/sn_fast.py b/sn_simulator/sn_fast.py index ac166894..aabbbfee 100644 --- a/sn_simulator/sn_fast.py +++ b/sn_simulator/sn_fast.py @@ -1,11 +1,9 @@ import numpy as np from astropy.table import Table import time -import pandas as pd -from sn_tools.sn_calcFast import LCfast, srand +from sn_tools.sn_calcFast import LCfast from sn_simu_wrapper.sn_object import SN_Object from sn_tools.sn_utils import SNTimer -from sn_tools.sn_io import dustmaps from astropy.coordinates import SkyCoord from dustmaps.sfd import SFDQuery @@ -23,18 +21,30 @@ class SN(SN_Object): """ - def __init__(self, param, simu_param, reference_lc=None, gamma=None, mag_to_flux=None, dustcorr=None, snr_fluxsec='',error_model=False): - super().__init__(param.name, param.sn_parameters, param.gen_parameters, - param.cosmology, param.telescope, param.SNID, param.area, param.x0_grid, + def __init__(self, param, simu_param, reference_lc=None, dustcorr=None): + super().__init__(param.name, + param.sn_parameters, + param.simulator_parameters, + param.gen_parameters, + param.cosmology, + param.zp_airmass, + param.SNID, + param.area, + param.x0_grid, param.salt2Dir, - mjdCol=param.mjdCol, RACol=param.RACol, DecCol=param.DecCol, - filterCol=param.filterCol, exptimeCol=param.exptimeCol, - m5Col=param.m5Col, seasonCol=param.seasonCol) + mjdCol=param.mjdCol, RACol=param.RACol, + DecCol=param.DecCol, + filterCol=param.filterCol, + exptimeCol=param.exptimeCol, + m5Col=param.m5Col, + seasonCol=param.seasonCol) # x1 and color are unique for this simulator x1 = np.unique(self.sn_parameters['x1']).item() color = np.unique(self.sn_parameters['color']).item() - + sn_type = self.sn_parameters['type'] + self.error_model = self.simulator_parameters['errorModel'] + #self.error_model_cut = self.simulator_parameters['errorModelCut'] """ # Loading reference file fname = '{}/LC_{}_{}_vstack.hdf5'.format( @@ -44,16 +54,19 @@ def __init__(self, param, simu_param, reference_lc=None, gamma=None, mag_to_flux fname, self.gammaFile, param.telescope) """ self.reference_lc = reference_lc - self.gamma = gamma - self.mag_to_flux = mag_to_flux + # self.gamma = gamma + # self.mag_to_flux = mag_to_flux # blue and red cutoffs are taken into account in the reference files # SN parameters for Fisher matrix estimation self.param_Fisher = ['x0', 'x1', 'color', 'daymax'] - bluecutoff = self.sn_parameters['blue_cutoff'] - redcutoff = self.sn_parameters['red_cutoff'] - self.lcFast = LCfast(reference_lc, dustcorr, x1, color, param.telescope, + bluecutoff = self.sn_parameters['blueCutoffg'] + redcutoff = self.sn_parameters['redCutoffg'] + + self.lcFast = LCfast(reference_lc, dustcorr, self.zp_slope, + self.zp_intercept, self.mean_wavelength, + x1, color, param.mjdCol, param.RACol, param.DecCol, param.filterCol, param.exptimeCol, param.m5Col, param.seasonCol, @@ -61,7 +74,8 @@ def __init__(self, param, simu_param, reference_lc=None, gamma=None, mag_to_flux bluecutoff=bluecutoff, redcutoff=redcutoff) - self.premeta = dict(zip(['x1', 'color', 'x0', ], [x1, color, -1.])) + self.premeta = dict( + zip(['x1', 'color', 'x0', 'sn_type'], [x1, color, -1., sn_type])) for vv in self.param_Fisher: vvv = 'epsilon_{}'.format(vv) dd = dict(zip([vvv], [np.unique(self.gen_parameters[vvv]).item()])) @@ -94,23 +108,12 @@ def __call__(self, obs, display=False, time_display=0): pixRA = np.mean(obs['pixRA']) pixDec = np.mean(obs['pixDec']) pixID = np.unique(obs['healpixID']).item() + season = np.mean(obs['season']) dL = -1 - + season_length = np.max(obs[self.mjdCol])-np.min(obs[self.mjdCol]) # get ebvofMW from dust maps - ebvofMW = self.sn_parameters['ebvofMW'] - if ebvofMW < 0.: - # in that case ebvofMW value is taken from a map - coords = SkyCoord(pixRA, pixDec, unit='deg') - try: - sfd = SFDQuery() - except Exception as err: - from dustmaps.config import config - config['data_dir'] = 'dustmaps' - import dustmaps.sfd - dustmaps.sfd.fetch() - - sfd = SFDQuery() - ebvofMW = sfd(coords) + ebvofMW = self.get_ebvofMW( + pixRA, pixDec, self.sn_parameters['ebvofMW']) # start timer ti = SNTimer(time.time()) @@ -119,17 +122,28 @@ def __call__(self, obs, display=False, time_display=0): np.array([b for b in 'grizy'])) if len(obs[goodFilters]) == 0: - return [self.nosim(ra, dec, pixRA, pixDec, pixID, season, ti, -1)] + return [self.nosim(RA, Dec, pixRA, pixDec, pixID, + season, season_length, ti, -1, ebvofMW)] tab_tot = self.lcFast(obs, ebvofMW, self.gen_parameters) + # remove LC points with too high error model value """ - # apply dust correction here - tab_tot = self.dust_corrections(tab_tot, pixRA, pixDec) + if self.error_model: + if self.error_model_cut >0: + idx = tab_tot['fluxerr_model']/tab_tot['flux']<= self.error_model_cut + tab_tot = tab_tot[idx] """ + + # apply dust correction here + # tab_tot = self.dust_corrections(tab_tot, pixRA, pixDec) + ptime = ti.finish(time.time())['ptime'].item() - self.premeta.update(dict(zip(['RA', 'Dec', 'pixRA', 'pixDec', 'healpixID', 'dL', 'ptime', 'status'], - [RA, Dec, pixRA, pixDec, pixID, dL, ptime, 1]))) + self.premeta.update(dict(zip(['RA', 'Dec', 'pixRA', 'pixDec', + 'healpixID', 'dL', 'ptime', + 'status', 'ebvofMW'], + [RA, Dec, pixRA, pixDec, pixID, dL, + ptime, 1, ebvofMW]))) """ ii = tab_tot['band'] == 'LSST::z' @@ -145,14 +159,52 @@ def __call__(self, obs, display=False, time_display=0): if display: for table_lc in list_tables: self.plotLC(table_lc['time', 'band', - 'flux', 'fluxerr', 'zp', 'zpsys'], time_display) + 'flux', 'fluxerr', + 'zp', 'zpsys'], time_display) - return list_tables + def get_ebvofMW(self, pixRA, pixDec, ebvofMW): + """ + Method to get E(B-V) MW. + + Parameters + ---------- + pixRA : float + RA position. + pixDec : float + Dec position. + ebvofMW : float + E(B-V). + + Returns + ------- + ebvofMW : float + E(B-V). + + """ + + if np.abs(ebvofMW) > 1.e-5: + # in that case ebvofMW value is taken from a map + coords = SkyCoord(pixRA, pixDec, unit='deg') + try: + sfd = SFDQuery() + except Exception as err: + print('dustmap not found', err) + from dustmaps.config import config + config['data_dir'] = 'dustmaps' + import dustmaps.sfd + dustmaps.sfd.fetch() + + sfd = SFDQuery() + ebvofMW = sfd(coords) + + return ebvofMW + def transform(self, tab): """ - Method to transform a pandas df to a set of astropytables with metadata + Method to transform a pandas df to a set of astropytables + with metadata Parameters --------------- @@ -165,26 +217,28 @@ def transform(self, tab): """ - groups = tab.groupby(['z', 'daymax']) + groups = tab.groupby(['z', 'daymax', 'season']) tab_tot = [] for name, grp in groups: newtab = Table.from_pandas(grp) - newtab.meta = dict(zip(['z', 'daymax'], name)) + newtab.meta = dict( + zip(['z', 'daymax', 'season'], name)) newtab.meta.update(self.premeta) tab_tot.append(newtab) return tab_tot - def nosim(self, RA, Dec, pixRA, pixDec, healpixID, season, ti, status): + def nosim(self, RA, Dec, pixRA, pixDec, healpixID, season, + season_length, ti, status, ebvofMW): """ Method to construct an empty table when no simulation was not possible Parameters --------------- - ra: float + RA: float SN RA - dec: float + Dec: float SN Dec pixRA: float pixel RA @@ -194,15 +248,19 @@ def nosim(self, RA, Dec, pixRA, pixDec, healpixID, season, ti, status): healpixID season: int season of interest + season_length: float + season length ptime: float processing time status: int status of the processing(1=ok, -1=no simu) - + ebvofMW: float + E(B-V) """ ptime = ti.finish(time.time())['ptime'].item() table_lc = Table() # set metadata table_lc.meta = self.metadata( - ra, dec, pix, area, season, ptime, snr_fluxsec, status) + RA, Dec, pixRA, pixDec, season, season_length, + ptime, status, ebvofMW) return table_lc diff --git a/tests/data_tests/ref_lc_config1.hdf5 b/tests/data_tests/ref_lc_config1.hdf5 new file mode 100644 index 00000000..8049a13e Binary files /dev/null and b/tests/data_tests/ref_lc_config1.hdf5 differ diff --git a/tests/data_tests/ref_lc_config2.hdf5 b/tests/data_tests/ref_lc_config2.hdf5 new file mode 100644 index 00000000..e3749f84 Binary files /dev/null and b/tests/data_tests/ref_lc_config2.hdf5 differ diff --git a/tests/data_tests/ref_lc_config3.hdf5 b/tests/data_tests/ref_lc_config3.hdf5 new file mode 100644 index 00000000..34314a40 Binary files /dev/null and b/tests/data_tests/ref_lc_config3.hdf5 differ diff --git a/tests/data_tests/ref_lc_config4.hdf5 b/tests/data_tests/ref_lc_config4.hdf5 new file mode 100644 index 00000000..e35fa80d Binary files /dev/null and b/tests/data_tests/ref_lc_config4.hdf5 differ diff --git a/tests/data_tests/ref_simu_config1.hdf5 b/tests/data_tests/ref_simu_config1.hdf5 new file mode 100644 index 00000000..97a21d21 Binary files /dev/null and b/tests/data_tests/ref_simu_config1.hdf5 differ diff --git a/tests/data_tests/ref_simu_config2.hdf5 b/tests/data_tests/ref_simu_config2.hdf5 new file mode 100644 index 00000000..2997e4aa Binary files /dev/null and b/tests/data_tests/ref_simu_config2.hdf5 differ diff --git a/tests/data_tests/ref_simu_config3.hdf5 b/tests/data_tests/ref_simu_config3.hdf5 new file mode 100644 index 00000000..f592fe2e Binary files /dev/null and b/tests/data_tests/ref_simu_config3.hdf5 differ diff --git a/tests/data_tests/ref_simu_config4.hdf5 b/tests/data_tests/ref_simu_config4.hdf5 new file mode 100644 index 00000000..11a1a388 Binary files /dev/null and b/tests/data_tests/ref_simu_config4.hdf5 differ diff --git a/tests/testSNsimulation.py b/tests/testSNsimulation.py index c1a6f2ef..a01e7164 100644 --- a/tests/testSNsimulation.py +++ b/tests/testSNsimulation.py @@ -1,20 +1,36 @@ +import pytest from builtins import zip import numpy as np import unittest -import lsst.utils.tests import os from numpy.testing import assert_almost_equal, assert_equal import h5py from astropy.table import Table, vstack from sn_simu_wrapper.sn_simu import SNSimulation from sn_tools.sn_io import check_get_file - +import yaml +from sn_simu_wrapper.config_simulation import ConfigSimulation main_repo = 'https://me.lsst.eu/gris/DESC_SN_pipeline' m5_ref = dict(zip('ugrizy', [23.60, 24.83, 24.38, 23.92, 23.35, 22.44])) def getFile(refdir, fname): + """ + Function to get file (web) named fname and located in refdir + + Parameters + ---------- + refdir : str + directory location of the file. + fname : str + file name. + + Returns + ------- + None. + + """ fullname = '{}/{}/{}'.format(main_repo, refdir, fname) # check whether the file is available; if not-> get it! @@ -25,149 +41,53 @@ def getFile(refdir, fname): def getRefDir(dirname): + """ + Function to grab a directory. + + Parameters + ---------- + dirname : str + directory name. + + Returns + ------- + None. + + """ fullname = '{}/{}'.format(main_repo, dirname) if not os.path.exists(dirname): print('wget path:', fullname) - cmd = 'wget --no-verbose --recursive {} --directory-prefix={} --no-clobber --no-parent -nH --cut-dirs=3 -R \'index.html*\''.format( + cmd = 'wget - -no-verbose - -recursive {} - -directory-prefix = {} \ + --no-clobber - -no-parent - nH - -cut-dirs = 3 - R \'index.html*\''.format( fullname+'/', dirname) os.system(cmd) -def getconfig(prodid, - x1Type, x1min, x1max, x1step, - colorType, colormin, colormax, colorstep, - zType, zmin, zmax, zstep, - daymaxType, daymaxstep, diffflux, - fulldbName, fieldType, fcoadd, seasval, - simulator='sn_cosmo', nside=64, nproc=1, outputDir='.'): - - config = {} - config['ProductionID'] = prodid+'_'+simulator - - # -------------- Supernova parameters ---------------------------------------- - config['SN parameters'] = {} - config['SN parameters']['Id'] = 100 # Id of the first SN - # stretch and color - config['SN parameters']['x1'] = {} - config['SN parameters']['x1']['type'] = x1Type # unique, uniform or random - config['SN parameters']['x1']['min'] = x1min - config['SN parameters']['x1']['max'] = x1max - config['SN parameters']['x1']['step'] = 1. - config['SN parameters']['color'] = {} - # unique, uniform or random - config['SN parameters']['color']['type'] = colorType - config['SN parameters']['color']['min'] = colormin - config['SN parameters']['color']['max'] = colormax - config['SN parameters']['color']['step'] = 0.05 - - config['SN parameters']['x1_color'] = {} +def Observations_band(day0=59000, daymin=59000, cadence=3., + season_length=140., band='r'): """ - config['SN parameters']['x1_color']['type'] = x1colorType # random or fixed - config['SN parameters']['x1_color']['min'] = [x1min, colormin] - config['SN parameters']['x1_color']['max'] = [0.2, 0.2] + Function to generate observations per band + + Parameters + ---------- + day0 : float, optional + Min day for obs. The default is 59000. + daymin : float, optional + Min day for obs. The default is 59000. + cadence : float. optional + cadence of observation (in days). The default is 3.. + season_length : float, optional + Season length (in days). The default is 140.. + band : str, optional + filter for obs. The default is 'r'. + + Returns + ------- + data : array + array of obs. + """ - config['SN parameters']['x1_color']['rate'] = 'JLA' - config['SN parameters']['x1_color']['dirFile'] = 'reference_files' - config['SN parameters']['z'] = {} # redshift - config['SN parameters']['z']['type'] = zType - config['SN parameters']['z']['min'] = zmin - config['SN parameters']['z']['max'] = zmax - config['SN parameters']['z']['step'] = zstep - # Type Ia volumetric rate = Perrett, Ripoche, Dilday. - config['SN parameters']['z']['rate'] = 'Perrett' - config['SN parameters']['daymax'] = {} # Tmax (obs. frame) - # unique, uniform or random - config['SN parameters']['daymax']['type'] = daymaxType - # if uniform: step (in day) between Tmin(obs) and Tmax(obs) - config['SN parameters']['daymax']['step'] = daymaxstep - # obs min and max phase (rest frame) for LC points - config['SN parameters']['min_rf_phase'] = - 20. - config['SN parameters']['max_rf_phase'] = 60. - # obs min and max phase (rest frame) for T0 estimation - config['SN parameters']['min_rf_phase_qual'] = -15 - config['SN parameters']['max_rf_phase_qual'] = 45 - config['SN parameters']['absmag'] = -19.0906 # peak abs mag - config['SN parameters']['band'] = 'bessellB' # band for absmag - config['SN parameters']['magsys'] = 'vega' # magsys for absmag - # to estimate differential flux - config['SN parameters']['differential_flux'] = diffflux - # dir where SALT2 files are located - config['SN parameters']['salt2Dir'] = 'SALT2_Files' - config['SN parameters']['blue_cutoff'] = 380. - config['SN parameters']['red_cutoff'] = 800. - # MW dust - config['SN parameters']['ebvofMW'] = 0 - # NSN factor - config['SN parameters']['NSN factor'] = 1 - - # ------------------cosmology ---------------------- - config['Cosmology'] = {} - config['Cosmology']['Model'] = 'w0waCDM' # Cosmological model - config['Cosmology']['Omega_m'] = 0.30 # Omega_m - config['Cosmology']['Omega_l'] = 0.70 # Omega_l - config['Cosmology']['H0'] = 72.0 # H0 - config['Cosmology']['w0'] = -1.0 # w0 - config['Cosmology']['wa'] = 0.0 # wa - - # -------------------instrument ----------------------- - config['Instrument'] = {} - config['Instrument']['name'] = 'LSST' # name of the telescope (internal) - # dir of throughput - config['Instrument']['throughput_dir'] = 'LSST_THROUGHPUTS_BASELINE' - config['Instrument']['atmos_dir'] = 'THROUGHPUTS_DIR' # dir of atmos - config['Instrument']['airmass'] = 1.2 # airmass value - config['Instrument']['atmos'] = True # atmos - config['Instrument']['aerosol'] = False # aerosol - - # -------------------observations ------------------------ - config['Observations'] = {} - # filename: /sps/lsst/cadence/LSST_SN_PhG/cadence_db/opsim_db/kraken_2026.db # Name of db obs file (full path) - config['Observations']['filename'] = fulldbName - config['Observations']['fieldtype'] = fieldType # DD or WFD - # this is the coaddition per night - config['Observations']['coadd'] = fcoadd - # season to simulate (-1 = all seasons) - config['Observations']['season'] = seasval - - # --------------------simulator ------------------------- - config['Simulator'] = {} - config['Simulator']['name'] = 'sn_simulator.{}'.format( - simulator) # Simulator name= sn_cosmo,sn_sim,sn_ana, sn_fast - config['Simulator']['model'] = 'salt2-extended' # spectra model - config['Simulator']['version'] = 1.0 # version - config['Simulator']['error_model'] = 0 # error model for fluxerr estimation - # Reference File= SN_MAF/Reference_Files/LC_Ref_-2.0_0.2.hdf5 - config['Simulator']['Template Dir'] = 'Template_LC' - config['Simulator']['Gamma Dir'] = 'reference_files' - config['Simulator']['Gamma File'] = 'gamma.hdf5' - config['Simulator']['DustCorr Dir'] = 'Template_Dust' - # -------------------------host --------------------- - config['Host Parameters'] = None # Host parameters - - # --------------------miscellanea ---------------- - config['Display_LC'] = {} # display during LC simulations - config['Display_LC']['display'] = False - # display during time (sec) before closing - config['Display_LC']['time'] = 30 - - config['Output'] = {} - config['Output']['directory'] = 'Output_Simu' - config['Output']['save'] = True - - config['Multiprocessing'] = {} - config['Multiprocessing']['nproc'] = nproc - - config['Metric'] = 'sn_mafsim.sn_maf_simulation' - - config['Pixelisation'] = {} - config['Pixelisation']['nside'] = nside - config['Web path'] = 'https://me.lsst.eu/gris/DESC_SN_pipeline' - - return config - - -def Observations_band(day0=59000, daymin=59000, cadence=3., season_length=140., band='r'): # Define fake data names = ['observationStartMJD', 'fieldRA', 'fieldDec', 'fiveSigmaDepth', 'visitExposureTime', 'numExposures', @@ -207,6 +127,24 @@ def Observations_band(day0=59000, daymin=59000, cadence=3., season_length=140., def Observations_season(day0=59000, mjdmin=59000, cadence=3.): + """ + Function to generate observations per season + + Parameters + ---------- + day0 : float, optional + First day of obs. The default is 59000. + mjdmin : float, optional + min day for obs. The default is 59000. + cadence : float, optional + cadence of observations (in days). The default is 3.. + + Returns + ------- + data : array + array of observations. + + """ bands = 'grizy' Nvisits = dict(zip(bands, [10, 20, 20, 26, 20])) rat = 34./3600./24. @@ -227,7 +165,8 @@ def Observations_season(day0=59000, mjdmin=59000, cadence=3.): for i in range(Nvisits[band]): mjd += shift dat = Observations_band( - daymin=mjd, season_length=season_length, cadence=cadence, band=band) + daymin=mjd, season_length=season_length, + cadence=cadence, band=band) if data is None: data = dat else: @@ -236,341 +175,422 @@ def Observations_season(day0=59000, mjdmin=59000, cadence=3.): return data -class TestSNsimulation(unittest.TestCase): +def fake_data(day0=59000, diff_season=280., nseasons=1): + """ + Function to generate fake data + + Parameters + ---------- + day0 : float, optional + First day of obs. The default is 59000. + diff_season : float, optional + delta time (in days) between two seasons. The default is 280.. + nseasons : int, optional + Number of seasons. The default is 1. + + Returns + ------- + data : array + array of obs. - def testSimuSNCosmo(self): - # set simulation parameters - prodid = 'Fake' - # x1colorType = 'unique' - x1Type = 'unique' - x1min = -2.0 - x1max = 2.0 - x1step = 0.1 - colorType = 'unique' - colormin = 0.2 - colormax = 0.3 - colorstep = 0.02 - zType = 'uniform' - zmin = 0.1 - zmax = 0.8 - zstep = 0.1 - daymaxtype = 'unique' - daymaxstep = 1. - difflux = 0 - fulldbName = 'data_from_fake' - fieldType = 'Fake' - fcoadd = 1 - seasval = [1] - - # get the config file from these - conf = getconfig(prodid, - x1Type, x1min, x1max, x1step, - colorType, colormin, colormax, colorstep, - zType, zmin, zmax, zstep, - daymaxtype, daymaxstep, difflux, - fulldbName, fieldType, fcoadd, seasval) - - # SN_Simulation instance - # getRefDir('SALT2_Files') - # getRefDir('reference_files') - - absMag = conf['SN parameters']['absmag'] - x0normFile = 'reference_files/X0_norm_{}.npy'.format(absMag) - - if not os.path.isfile(x0normFile): - # if this file does not exist, grab it from a web server - check_get_file(conf['Web path'], 'reference_files', - 'X0_norm_{}.npy'.format(absMag)) - x0_norm = np.load(x0normFile) - - area = 9.6 # survey area (deg2) - - simu = SNSimulation(mjdCol='observationStartMJD', - filterCol='filter', - nexpCol='numExposures', - exptimeCol='visitExposureTime', - config=conf, x0_norm=x0_norm) - - # Generate fake data - day0 = 59000 - data = None - - diff_season = 280. - nseasons = 1 - for val in np.arange(59000, 59000+nseasons*diff_season, diff_season): - dat = Observations_season(day0, val) - if data is None: - data = dat - else: - data = np.concatenate((data, dat)) + """ - # now simulate LC on this data + # Generate fake data - simu.run(data) + data = None - # save metadata + for val in np.arange(59000, 59000+nseasons*diff_season, diff_season): + dat = Observations_season(day0, val) + if data is None: + data = dat + else: + data = np.concatenate((data, dat)) - simu.save_metadata() + return data - # check what we have inside the data - simu_name = '{}/Simu_{}_1.hdf5'.format( - conf['Output']['directory'], conf['ProductionID']) - lc_name = '{}/LC_{}_1.hdf5'.format( - conf['Output']['directory'], conf['ProductionID']) +def getSimu(config_name): + """ + Function to perform simulations - f = h5py.File(simu_name, 'r') - # reading the simu file - simul = Table() - for i, key in enumerate(f.keys()): - simul = vstack([simul, Table.read(simu_name, path=key)]) + Parameters + ---------- + config_name : yaml file + Simulation parameters. - # first check on simulation parameters + Returns + ------- + simu : class instance + instance of SNSimulation. + conf : dict + suimulation parameters. - """ - for cc in simul.columns: - print('RefDict[\'{}\']='.format(cc), simul[cc].tolist()) - """ - RefDict = {} - RefDict['SNID'] = [11, 12, 13, 14, 15, 16, 17, 18] - RefDict['index_hdf5'] = ['10_-2.0_0.2_0.1_59023.102_1_0', '10_-2.0_0.2_0.2_59025.202_1_0', '10_-2.0_0.2_0.3_59027.302_1_0', '10_-2.0_0.2_0.4_59029.402_1_0', - '10_-2.0_0.2_0.5_59031.502_1_0', '10_-2.0_0.2_0.6_59033.602_1_0', '10_-2.0_0.2_0.7_59035.702_1_0', '10_-2.0_0.2_0.8_59037.802_1_0'] - RefDict['season'] = [1, 1, 1, 1, 1, 1, 1, 1] - RefDict['fieldname'] = ['unknown', 'unknown', 'unknown', - 'unknown', 'unknown', 'unknown', 'unknown', 'unknown'] - RefDict['fieldid'] = [0, 0, 0, 0, 0, 0, 0, 0] - RefDict['area'] = [0.8392936452111668, 0.8392936452111668, 0.8392936452111668, 0.8392936452111668, - 0.8392936452111668, 0.8392936452111668, 0.8392936452111668, 0.8392936452111668] - RefDict['RA'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['Dec'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['x0'] = [0.0001740043427228556, 3.838217653816398e-05, 1.5291775006387726e-05, 7.813694452227428e-06, - 4.593754146993177e-06, 2.958319253169902e-06, 2.031812357572594e-06, 1.4642465313185475e-06] - RefDict['epsilon_x0'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['x1'] = [-2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0] - RefDict['epsilon_x1'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['color'] = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2] - RefDict['epsilon_color'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['daymax'] = [59023.10190972222, 59025.20190972222, 59027.30190972223, 59029.401909722226, - 59031.501909722225, 59033.60190972222, 59035.70190972222, 59037.80190972223] - RefDict['epsilon_daymax'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['z'] = [0.1, 0.2, 0.30000000000000004, - 0.4, 0.5, 0.6, 0.7000000000000001, 0.8] - RefDict['survey_area'] = [0.8392936452111668, 0.8392936452111668, 0.8392936452111668, - 0.8392936452111668, 0.8392936452111668, 0.8392936452111668, 0.8392936452111668, 0.8392936452111668] - RefDict['healpixID'] = [10, 10, 10, 10, 10, 10, 10, 10] - RefDict['pixRA'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['pixDec'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['dL'] = [447513.8270462983, 952843.778897064, 1509584.9325567451, 2111826.765549938, - 2754245.369069551, 3432131.997266213, 4141376.3917280077, 4878422.389153463] - RefDict['snr_fluxsec_meth'] = ['interp', 'interp', 'interp', - 'interp', 'interp', 'interp', 'interp', 'interp'] - RefDict['status'] = [1, 1, 1, 1, 1, 1, 1, 1] - RefDict['ebvofMW'] = [0, 0, 0, 0, 0, 0, 0, 0] - - for key, vv in RefDict.items(): - if key not in ['index_hdf5', 'fieldname', 'snr_fluxsec_meth']: - assert(np.isclose(vv, simul[key].tolist()).all()) - else: - assert((vv == simul[key]).all()) + """ - # now grab LC + # get the config file from these + ffi = open(config_name) + conf = yaml.load(ffi, Loader=yaml.FullLoader) + ffi.close() - vars = ['snr_m5', 'flux_e_sec', 'mag', - 'exptime', 'magerr', 'band', 'phase','snr_tot'] + # print('conf',conf) + absMag = conf['SN']['absmag'] + x0normFile = 'reference_files/X0_norm_{}.npy'.format(absMag) - for simu in simul: - lc = Table.read(lc_name, path='lc_{}'.format(simu['index_hdf5'])) - idx = lc['snr_m5'] >= 5. - lc = lc[idx][:20] - break + if not os.path.isfile(x0normFile): + # if this file does not exist, grab it from a web server + check_get_file(conf['WebPathSimu'], 'reference_files', + 'X0_norm_{}.npy'.format(absMag)) + x0_norm = np.load(x0normFile) - """ - for cc in vars: - print('RefDict[\'{}\']='.format(cc), lc[cc].tolist()) - """ + simu = SNSimulation(mjdCol='observationStartMJD', + filterCol='filter', + nexpCol='numExposures', + exptimeCol='visitExposureTime', + config=conf, x0_norm=x0_norm) - RefDict = {} - - RefDict['snr_m5'] = [27.77258180379027, 18.98377096566505, 72.79590570928633, 152.87684647830162, 156.06666693213177, 29.621529233519308, 285.0437939047469, 367.04307444774975, 359.1500294310009, - 171.59301350412957, 506.58512140687225, 646.2792675349654, 598.4544644527755, 288.72360876534253, 661.3175452293224, 849.9877278958845, 734.1754283838669, 389.3160994765804, 745.0362859344243, 966.9961020792483] - RefDict['flux_e_sec'] = [38.63207940525376, 30.058075022627055, 121.90767409695435, 229.50791298666263, 265.77760354105226, 47.616478031948624, 597.1868291509262, 626.8365327579263, 680.7371214592486, - 292.13548602917916, 1321.4612615401911, 1298.4282257064317, 1283.1179646726964, 515.1858166393667, 1986.2944122972094, 1911.7598017542887, 1684.8650500083472, 723.2044598429111, 2404.500751077161, 2314.677056897036] - RefDict['mag'] = [24.197549880640725, 24.157628846869798, 23.160950966140664, 22.263092425796074, 21.79125631239912, 23.244278569706797, 21.435793513964224, 21.17226381616802, 20.77021232574782, - 21.27488067257211, 20.5734837320662, 20.381589712793957, 20.081984493733383, 20.658809391021087, 20.130996540814564, 19.96148438529608, 19.7862032720865, 20.290578693537636, 19.92355764391396, 19.75371788930898] - RefDict['exptime'] = [600.0, 600.0, 300.0, 600.0, 600.0, 780.0, 300.0, 600.0, 600.0, - 780.0, 300.0, 600.0, 600.0, 780.0, 300.0, 600.0, 600.0, 780.0, 300.0, 600.0] - RefDict['magerr'] = [0.039093816067541594, 0.05719286261522241, 0.014914797668622531, 0.007102031666464496, 0.006956874431298518, 0.03665361758330578, 0.003809015414385587, 0.0029580620922809125, 0.0030230714625814016, 0.0063273916728084035, - 0.0021432453478753125, 0.0016799799394762239, 0.0018142336121611566, 0.0037604690846066266, 0.0016417774072236856, 0.0012773551536394909, 0.0014788511883980504, 0.0027888294530276484, 0.001457293054386471, 0.001122792741794475] - RefDict['band'] = ['LSST::r', 'LSST::i', 'LSST::g', 'LSST::r', 'LSST::i', 'LSST::z', 'LSST::g', 'LSST::r', 'LSST::i', - 'LSST::z', 'LSST::g', 'LSST::r', 'LSST::i', 'LSST::z', 'LSST::g', 'LSST::r', 'LSST::i', 'LSST::z', 'LSST::g', 'LSST::r'] - RefDict['phase'] = [-15.54029882153951, -15.536721380459229, -12.818181818180495, -12.813026094266784, -12.809448653186502, -12.808501683486444, -10.090909090907767, -10.085753366994057, -10.082175925913774, - - 10.081228956213717, -7.36363636363504, -7.35848063972133, -7.354903198641047, -7.353956228940991, -4.636363636362313, -4.631207912448603, -4.62763047136832, -4.626683501668263, -1.909090909089586, -1.9039351851758757] - - for key, vv in RefDict.items(): - #print('testing', key) - if key not in ['band']: - assert(np.isclose(vv, lc[key].tolist()).all()) - else: - assert(set(vv) == set(lc[key].tolist())) + return simu, conf - def testSimuSNFast(self): - # set simulation parameters - prodid = 'Fake' - # x1colorType = 'unique' - x1Type = 'unique' - x1min = -2.0 - x1max = 2.0 - x1step = 0.1 - colorType = 'unique' - colormin = 0.2 - colormax = 0.3 - colorstep = 0.02 - zType = 'uniform' - zmin = 0.1 - zmax = 0.8 - zstep = 0.1 - daymaxtype = 'unique' - daymaxstep = 1. - difflux = 0 - fulldbName = 'data_from_fake' - fieldType = 'Fake' - fcoadd = 1 - seasval = [1] - - # get the config file from these - conf = getconfig(prodid, - x1Type, x1min, x1max, x1step, - colorType, colormin, colormax, colorstep, - zType, zmin, zmax, zstep, - daymaxtype, daymaxstep, difflux, - fulldbName, fieldType, fcoadd, seasval, - simulator='sn_fast') - - # get the reference LC file - - #referenceName = 'LC_{}_{}_vstack.hdf5'.format(x1min, colormin) - #getFile('Templates', referenceName) - - # instance of SNSimulation - simu = SNSimulation(mjdCol='observationStartMJD', - filterCol='filter', - nexpCol='numExposures', - exptimeCol='visitExposureTime', - config=conf) - - # Generate fake data - day0 = 59000 - data = None - - diff_season = 280. - nseasons = 1 - for val in np.arange(59000, 59000+nseasons*diff_season, diff_season): - dat = Observations_season(day0, val) - if data is None: - data = dat + +def dump(fname, thedict): + """ + Function to dump a dict in a yaml file. + + Parameters + ---------- + fname : str + output file name. + thedict : dict + dict to dum in the yaml file. + + Returns + ------- + None. + + """ + + # with open(fname, 'w') as ffile: + ffile = open(fname, 'w') + documents = yaml.dump(thedict, ffile) + ffile.close() + + +def simulation(data, config_name): + """ + Functio to perform LC simulations + + Parameters + ---------- + data : array + observations used to simulate LCs. + config_name : str + yaml config name. + + Returns + ------- + None. + + """ + + simu, conf = getSimu(config_name) + + # now simulate LC on this data + + simu.run(data) + + # save metadata + + simu.save_metadata() + + +def gimeSimu(simu_name): + + f = h5py.File(simu_name, 'r') + # reading the simu file + simul = Table() + for i, key in enumerate(f.keys()): + simul = vstack([simul, Table.read(simu_name, path=key)]) + + return simul + + +def gimerefSimu(config_name): + """ + Functio to get reference simu params + + Parameters + ---------- + config_name : str + config file name. + + Returns + ------- + tab_simu_ref : astropy table + Reference simulation parameters. + + """ + + name_config = config_name.split('/')[-1].split('.')[0] + + ref_simul = 'data_tests/ref_simu_{}.hdf5'.format(name_config) + print('alors:::', ref_simul) + # load reference parameters + tab_simu_ref = Table.read(ref_simul, path='simu_parameters') + if 'ptime' in tab_simu_ref.columns: + tab_simu_ref.remove_column('ptime') + + return tab_simu_ref + + +def simuName(conf): + """ + get the name of the simu parameter file + + Parameters + ---------- + conf : dict + simu parameters. + + Returns + ------- + simu_name : str + simu file name. + + """ + + simu_name = '{}/Simu_{}_1.hdf5'.format( + conf['OutputSimu']['directory'], conf['ProductionIDSimu']) + + return simu_name + + +def lcName(conf): + """ + get the name of the lc file + + Parameters + ---------- + conf : dict + simu params. + + Returns + ------- + lc_name : astropy table + light curve points + + """ + + lc_name = '{}/LC_{}_1.hdf5'.format( + conf['OutputSimu']['directory'], conf['ProductionIDSimu']) + + return lc_name + + +def conftest(data, config, fname, + cols_simpars=['sn_model', 'sn_version', 'daymax', + 'z', 'ebvofMW', 'x0', 'x1', 'color'], + cols_lc=['snr_m5', 'flux', 'mag', + 'exptime', 'magerr', 'band', 'phase']): + """ + Function to perform unit tests + + Parameters + ---------- + data : array + observations used in simulations. + config : dict + simu parameters. + fname : str + yaml output file name. + cols_simpars : list(str), optional + list of columns for comparison of simu params. + The default is ['sn_model', 'sn_version', 'daymax', + 'z', 'ebvofMW', 'x0', 'x1', 'color']. + cols_lc : list(str), optional + list of columns for LC points comparison. + The default is ['snr_m5', 'flux', 'mag', + 'exptime', 'magerr', 'band', 'phase']. + + Returns + ------- + None. + + """ + + dump(fname, config) + name_config = fname.split('/')[-1].split('.')[0] + + ##################### + # Simu parameters # + ##################### + + # load simu data + simulation(data, fname) + simu_name = simuName(config) + simu_data = gimeSimu(simu_name) + print('kkk', simu_data) + ref_simul = 'data_tests/ref_simu_{}.hdf5'.format(name_config) + + if not os.path.exists(ref_simul): + simu_data.write(ref_simul, 'simu_parameters', + append=True, compression=True) + + refSimu = gimerefSimu(fname) + + # check simu parameters here + @ pytest.mark.parametrize("simu_data,expected", refSimu) + def testsimu_data(simu_data, expected): + print('ccol', expected.columns) + print(simu_data['n_lc_points'], expected['n_lc_points']) + print(type(simu_data), type(expected)) + for key in cols_simpars: + if key not in ['index_hdf5', 'fieldname', 'snr_fluxsec_meth', + 'sn_type', 'sn_model', 'sn_version', 'SNID']: + assert(np.isclose( + expected[key].tolist(), simu_data[key].tolist()).all()) else: - data = np.concatenate((data, dat)) + assert((expected[key] == simu_data[key]).all()) + testsimu_data(simu_data, refSimu) + + #################### + # LC data # + #################### + + lc_name = lcName(config) + + for simu in simu_data: + lc = Table.read(lc_name, path='lc_{}'.format(simu['index_hdf5'])) + idx = lc['snr_m5'] >= 5. + lc = lc[idx][: 20] + break + + ref_lc = 'data_tests/ref_lc_{}.hdf5'.format(name_config) + if not os.path.exists(ref_lc): + lc.write(ref_lc, 'lc_points', append=True, compression=True) + + # load lc reference points + tab_lc_ref = Table.read(ref_lc, path='lc_points') + + if 'index' in tab_lc_ref.columns: + tab_lc_ref.remove_column('index') + + # check lc data here - # now simulate LC on this data + @ pytest.mark.parametrize("lc_data,expected", tab_lc_ref) + def testsimu_lc(lc_data, expected): - tab = simu.run(data) + for key in cols_lc: + if key not in ['band', 'filter_cosmo', 'zpsys']: + assert(np.isclose( + expected[key].tolist(), lc_data[key].tolist()).all()) + else: + assert(set(expected[key].tolist()) + == set(lc_data[key].tolist())) + testsimu_lc(lc, tab_lc_ref) - # save metadata - simu.save_metadata() +class TestSNsimulation(unittest.TestCase): + """ + class to perform unit tests. + """ - # check what we have inside the data + def testSimuSNCosmo(self): + """ + Method to test simulated data with sncosmo - simu_name = '{}/Simu_{}_1.hdf5'.format( - conf['Output']['directory'], conf['ProductionID']) - lc_name = '{}/LC_{}_1.hdf5'.format( - conf['Output']['directory'], conf['ProductionID']) + Returns + ------- + None. - f = h5py.File(simu_name, 'r') - # reading the simu file - simul = Table() - for i, key in enumerate(f.keys()): - simul = vstack([simul, Table.read(simu_name, path=key)]) + """ - # first check on simulation + # fake data + data = fake_data() + + # test Ia simulation + # test Ia - no error model + # get configuration file + config = ConfigSimulation( + 'Ia', 'salt3', '../sn_simu_input/config_simulation.txt').conf_dict + config['ProductionIDSimu'] = 'prod_Ia_sncosmo_errormodel_0' + config['Simulator']['errorModel'] = 0 + config['SN']['z']['type'] = 'uniform' + config['SN']['z']['step'] = 0.1 + config['SN']['z']['min'] = 0.1 + config['SN']['z']['max'] = 1.0 + config['SN']['NSNabsolute'] = 1 + config['SN']['ebvofMW'] = 0.0 + config['ReferenceFiles']['GammaFile'] = 'gamma_DDF.hdf5' """ - for cc in simul.columns: - print('RefDict[\'{}\']='.format(cc), simul[cc].tolist()) + config['SN']['x1']['type'] = 'random' + config['SN']['color']['type'] = 'random' + config['Display']['LC']['display'] = 1 """ - RefDict = {} - - RefDict['SNID'] = [11, 12, 13, 14, 15, 16, 17, 18] - RefDict['index_hdf5'] = ['10_-2.0_0.2_0.1_59023.1_1', '10_-2.0_0.2_0.2_59025.2_1', '10_-2.0_0.2_0.3_59027.3_1', '10_-2.0_0.2_0.4_59029.4_1', - '10_-2.0_0.2_0.5_59031.5_1', '10_-2.0_0.2_0.6_59033.6_1', '10_-2.0_0.2_0.7_59035.7_1', '10_-2.0_0.2_0.8_59037.8_1'] - RefDict['season'] = [1, 1, 1, 1, 1, 1, 1, 1] - RefDict['fieldname'] = ['unknown', 'unknown', 'unknown', - 'unknown', 'unknown', 'unknown', 'unknown', 'unknown'] - RefDict['fieldid'] = [0, 0, 0, 0, 0, 0, 0, 0] - RefDict['n_lc_points'] = [88, 96, 104, 112, 120, 128, 102, 108] - RefDict['area'] = [0.8392936452111668, 0.8392936452111668, 0.8392936452111668, 0.8392936452111668, - 0.8392936452111668, 0.8392936452111668, 0.8392936452111668, 0.8392936452111668] - RefDict['z'] = [0.1, 0.2, 0.30000000000000004, - 0.4, 0.5, 0.6, 0.7000000000000001, 0.8] - RefDict['daymax'] = [59023.10190972222, 59025.20190972222, 59027.30190972223, 59029.401909722226, - 59031.501909722225, 59033.60190972222, 59035.70190972222, 59037.80190972223] - RefDict['x1'] = [-2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0, -2.0] - RefDict['color'] = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2] - RefDict['x0'] = [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0] - RefDict['epsilon_x0'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['epsilon_x1'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['epsilon_color'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['epsilon_daymax'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['RA'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['Dec'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['pixRA'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['pixDec'] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - RefDict['healpixID'] = [10, 10, 10, 10, 10, 10, 10, 10] - RefDict['dL'] = [-1, -1, -1, -1, -1, -1, -1, -1] - - for key, vv in RefDict.items(): - if key not in ['index_hdf5', 'fieldname']: - assert(np.isclose(vv, simul[key].tolist()).all()) - - # now grab LC - - vars = ['snr_m5', 'flux_e_sec', 'mag', - 'visitExposureTime', 'magerr', 'band', 'phase'] - - for simu in simul: - lc = Table.read(lc_name, path='lc_{}'.format(simu['index_hdf5'])) - idx = lc['snr_m5'] >= 5. - lc = lc[idx][:20] - break + fname = 'config2.yaml' + conftest(data, config, fname) + + print('running with error model') + # test Ia - with error model + config['Simulator']['errorModel'] = 1 + config['ProductionIDSimu'] = 'prod_Ia_sncosmo_errormodel_1' + fname = 'config1.yaml' + conftest(data, config, fname) + + # test non Ia - error_model=0 + config['ProductionIDSimu'] = 'prod_Ib_sncosmo_errormodel_0' + error_model = 0 + sn_type = 'SN_Ib' + sn_model = 'nugent-sn2p' + sn_model_version = '1.2' + config['Simulator']['model'] = sn_model + config['Simulator']['version'] = sn_model_version + config['SN']['type'] = sn_type + config['SN']['z']['type'] = 'uniform' + config['SN']['z']['step'] = 0.1 + config['SN']['z']['min'] = 0.01 + config['SN']['z']['max'] = 1.0 + config['Simulator']['errorModel'] = error_model + + # test Ib + fname = 'config3.yaml' + conftest(data, config, fname, + cols_simpars=['sn_model', 'sn_version', 'daymax', + 'z', 'ebvofMW']) + def testSimuSNFast(self): """ - for cc in vars: - print('RefDict[\'{}\']='.format(cc), lc[cc].tolist()) + Method to test SN fast simu output + + Returns + ------- + None. + """ - RefDict = {} - - RefDict['snr_m5'] = [72.76621556433962, 284.91155929878005, 506.3382490861352, 661.0244890696117, 744.7374346856301, 750.9560627892503, 701.5801119418779, 617.0362767220539, 509.49758960144334, - 405.6704063879359, 312.1587215082942, 236.20474351258605, 185.20525529387027, 156.7007460981387, 123.82556613113591, 98.86188081087201, 89.07300629269548, 85.39942468306856, 81.6029185267763, 77.71841283547786] - RefDict['flux_e_sec'] = [121.85475114199077, 596.8249003571814, 1320.5101003419977, 1984.8958672777107, 2402.9349000717466, 2435.4188847893047, 2182.188549348978, 1781.5190648980197, 1332.6543244852983, - 959.9063595202888, 672.3614552092635, 470.38517751045663, 349.6450546810812, 287.0414484278593, 219.024995436545, 170.22546308092205, 151.80546033037837, 144.97271315775106, 137.96525588832066, 130.84706269377557] - RefDict['mag'] = [23.161428491773833, 21.436445105957933, 20.574266452603773, 20.131756917960203, 19.9242657576835, 19.90957575410812, 20.028884191267565, 20.248882689684862, 20.564272843036967, 20.920539051651428, - 21.30696280634007, 21.694807732877024, 22.017022106232524, 22.23124270951677, 22.524886922545807, 22.798310177888723, 22.922900267535, 22.972903089242255, 23.02669579575968, 23.084173861687447] - RefDict['visitExposureTime'] = [300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, - 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0, 300.0] - RefDict['magerr'] = [0.014920883219467768, 0.003810783274045906, 0.0021442903172290873, 0.0016425052667659818, 0.00145787784283523, 0.0014458052322334502, 0.0015475584131838627, 0.0017595986584873079, 0.0021309938003974682, - 0.00267639982523117, 0.0034781543168553785, 0.004596589334372434, 0.0058623401535521326, 0.006928723900766582, 0.008768271679923466, 0.010982354329624782, 0.012189284385331971, 0.012713624345684724, 0.013305114870393606, 0.013970128379442394] - RefDict['band'] = ['LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', - 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g', 'LSST::g'] - RefDict['phase'] = [-12.818181818180495, -10.090909090907767, -7.36363636363504, -4.636363636362313, -1.909090909089586, 0.8181818181831411, 3.545454545455868, 6.272727272728595, 9.000000000001322, - 11.72727272727405, 14.454545454546777, 17.181818181819505, 19.90909090909223, 22.63636363636496, 25.363636363637685, 28.09090909091041, 30.818181818183138, 33.54545454545587, 36.272727272728595, 39.00000000000132] - - for key, vv in RefDict.items(): - if key not in ['band']: - assert(np.isclose(vv, lc[key].tolist()).all()) - else: - assert(set(vv) == set(lc[key].tolist())) + # fake data + data = fake_data() + print('data', data) + + config = ConfigSimulation( + 'Ia', 'salt2', '../sn_simu_input/config_simulation.txt').conf_dict + config['ProductionIDSimu'] = 'prod_Ia_snfast_errormodel_0' + config['Simulator']['errorModel'] = 0 + config['SN']['z']['type'] = 'uniform' + config['SN']['z']['step'] = 0.1 + config['SN']['z']['min'] = 0.1 + config['SN']['z']['max'] = 1.0 + config['SN']['ebvofMW'] = 0.0 + config['SN']['NSNabsolute'] = 1 + + config['Simulator']['name'] = 'sn_simulator.sn_fast' + config['ReferenceFiles']['GammaFile'] = 'gamma_DDF.hdf5' + + print('ccc', config) + fname = 'config4.yaml' + conftest(data, config, fname, cols_simpars=['season', 'daymax', + 'z', 'ebvofMW'], + cols_lc=['snr_m5', 'flux', 'mag', + 'magerr', 'band', 'phase']) if __name__ == "__main__": - lsst.utils.tests.init() unittest.main(verbosity=5) diff --git a/version.py b/version.py deleted file mode 100644 index 0be4d125..00000000 --- a/version.py +++ /dev/null @@ -1 +0,0 @@ -__version__='v1.1.0' diff --git a/version_deprecated.py b/version_deprecated.py new file mode 100644 index 00000000..1ab67ec7 --- /dev/null +++ b/version_deprecated.py @@ -0,0 +1 @@ +__version__='1.1.0'