diff --git a/examples/HowToDelightPZ.md b/examples/HowToDelightPZ.md new file mode 100644 index 00000000..0375583d --- /dev/null +++ b/examples/HowToDelightPZ.md @@ -0,0 +1,131 @@ +# DelightPZ README.md + +- sylvielsstfr, + - creation date : Feb 24th 2021 + - update : May 9th 2021 + +## Delight reference package + +- You must have a working Delight package installed. + +- Get it from https://github.com/LSSTDESC/Delight + +- Then install Delight according the instructions here + +- **https://delight.readthedocs.io/en/latest/install.html** + + +Get the Delight data https://github.com/LSSTDESC/Delight/tree/master/data + +in some path defined in RAIL config file (see file **DelightPZ.yaml**). + + + +## LSSTDESC/RAIL + +### Installation of RAIL + +Installation from here +- **https://github.com/LSSTDESC/RAIL** + + git clone git@github.com:LSSTDESC/RAIL.git + cd RAIL + + +For the moment, the following files are under developpment in the branch **issue/49/delight** + + git checkout issue/49/delight + + +- **RAIL/examples/configs/delightPZ.yaml** +- **RAIL/rail/estimation/algos/delightPZ.py** +- **RAIL/rail/estimation/algos/include_delightPZ/** + +Note **include_delightPZ** includes scripts interfaces from Delight. + + + + rail + ├── README.md + ├── __init__.py + ├── creation + │   ├── README.md + │   └── __init__.py + ├── estimation + │   ├── README.md + │   ├── __init__.py + │   ├── algos + │   │   ├── __init__.py + │   │   ├── delightPZ.py + │   │   ├── flexzboost.py + │   │   ├── include_delightPZ + │   │   │   ├── __init__.py + │   │   │   ├── calibrateTemplateMixturePriors.py + │   │   │   ├── convertDESCcat.py + │   │   │   ├── delightApply.py + │   │   │   ├── delightLearn.py + │   │   │   ├── getDelightRedshiftEstimation.py + │   │   │   ├── makeConfigParam.py + │   │   │   ├── processFilters.py + │   │   │   ├── processSEDs.py + │   │   │   ├── simulateWithSEDs.py + │   │   │   └── templateFitting.py + │   │   ├── randomPZ.py + │   │   ├── sklearn_nn.py + │   │   └── trainZ.py + │   ├── estimator.py + │   └── utils.py + └── evaluation + ├── README.md + └── __init__.py + + +### Build RAIL + + + python setup.py install + + +### Run RAIL + + +- delightPZ generates the configuration file **parametersTest.cfg** required by Delight + +The path of the configuration file is defined in **DelightPZ.yaml** (RAIL/example/configs) . + +#### RAIL temporay files + +The recommended temporary structure in **examples/** directory from where one issue the command to run RAILS with Delight: + + python main.py config/delightPZ.yaml + +- Some temporary directories must be created in **examples/** for Delight input and output data. + +- The recommended structure is the following + +#### For DC2 data +tree tmp + + tmp + ├── delight_data + ├── delight_indata + │   ├── BROWN_SEDs + │   ├── CWW_SEDs + │   └── FILTERS + + + +#### For internal mock data (simulation) +tree tmpsim + + tmpsim + ├── delight_data + ├── delight_indata + │   ├── BROWN_SEDs + │   ├── CWW_SEDs + │   └── FILTERS + + +Note the directories BROWN_SEDs, CWW_SEDs, FILTERS must be copied from Delight installation +https://github.com/LSSTDESC/Delight/tree/master/data . + diff --git a/examples/configs/delightPZ.yaml b/examples/configs/delightPZ.yaml new file mode 100644 index 00000000..bc3eaeba --- /dev/null +++ b/examples/configs/delightPZ.yaml @@ -0,0 +1,76 @@ +########################################################################################## +# +# RAIL configuration file for delightPZ module +# +# Steering Delight from RAIL +# Used on Vera C. Rubin LSST only estimation +# +# Author : Sylvie Dagoret-Campagne +# Affiliation : IJCLab/IN2P3/CNRS/France +# Creation date : March 2021 +# Last update : May 8th 2021 +# +############################################################################################ +run_params: + class_name: delightPZ + run_name: test_delightPZ +#------------------------------------------------ +# redshift range and binning for delight +# dlght_ prepend means a parameter used by Delight +#------------------------------------------------ + dlght_redshiftMin: 0.01 + dlght_redshiftMax: 3.01 + dlght_redshiftNumBinsGPpred: 300 + dlght_redshiftBinSize: 0.01 + dlght_redshiftDisBinSize: 0.2 +#----------------------------------------------- +# Delight input data (Filters and SED templates), note Delight will write inside this directory +# example recommended: +# - DC2 mode : dlght_inputdata: "./tmp/delight_indata" +# - tutorial mode : dlght_inputdata: "./tmpsim/delight_indata" +#---------------------------------------------- + dlght_inputdata: "./tmp/delight_indata" +#--------------------------------------- +# temporary directory for delight work +# example recommended: +# - DC2 mode : tempdir: "./tmp" and tempdatadir: "./tmp/delight_data" +# - tutorial mode : tempdir: "./tmpsim" and tempdatadir: "./tmpsim/delight_data" +#--------------------------------------- + tempdir: "./tmp" + tempdatadir: "./tmp/delight_data" +#---------------------------------------- +# delight configuration file filename +#---------------------------------------- + delightparamfile: "parametersTest.cfg" +#------------------------------------------------------------------------ +# Running mode +# tutorial mode : +# - True : activate internal flux simulation +# - False : activate use of DC2 datasets for training and Validation +#----------------------------------------------------------------------- + dlght_tutorialmode: False +#------------------------------------------------------------------------ +# Filtering of training and Validation dataset +# according flux SNR +#------------------------------------------------------------------------- + flag_filter_training: True + snr_cut_training: 5 + flag_filter_validation: True + snr_cut_validation: 3 +#----------------------------------------------------------------------- +# Special run +# Should not be used +#------------------------------------------------------------------------- + dlght_calibrateTemplateMixturePrior: False +#------------------------------------------------------------------------- +# Delight hyper-parameters that must be optimized +#--------------------------------------------------------------------------- + zPriorSigma: 0.2 + ellPriorSigma: 0.5 + fluxLuminosityNorm: 1.0 + alpha_C: 1.0e3 + V_C: 0.1 + alpha_L: 1.0e2 + V_L: 0.1 + lineWidthSigma: 20 + diff --git a/rail/estimation/algos/delightPZ.py b/rail/estimation/algos/delightPZ.py new file mode 100644 index 00000000..f5f4ad9c --- /dev/null +++ b/rail/estimation/algos/delightPZ.py @@ -0,0 +1,321 @@ +########################################################################################## +# +# RAIL interface class to Delight +# +# Steering Delight from RAIL +# Used on Vera C. Runbin LSST only estimation +# +# Author : Sylvie Dagoret-Campagne +# Affiliation : IJCLab/IN2P3/CNRS/France +# Creation date : March 2021 +# Last update : March 21th 2021 +# +############################################################################################ + +import numpy as np +from scipy.stats import norm +from rail.estimation.estimator import Estimator as BaseEstimation + +import os +import errno + +import coloredlogs +import logging +from pkg_resources import resource_filename + +# Delight initialisation + +#from interfaces.rail.processFilters import processFilters # interface added into delight in branch rail +from rail.estimation.algos.include_delightPZ.processFilters import processFilters +#from interfaces.rail.processSEDs import processSEDs # build a redshift -flux grid model +from rail.estimation.algos.include_delightPZ.processSEDs import processSEDs # build a redshift -flux grid model + +# interface with Delight through files +#from interfaces.rail.makeConfigParam import * # build the parameter file required by Delight +from rail.estimation.algos.include_delightPZ.makeConfigParam import * # build the parameter file required by Delight + +#from interfaces.rail.convertDESCcat import convertDESCcat # convert DESC input file into Delight format +#from interfaces.rail.convertDESCcat import * # convert DESC input file into Delight format +#from from rail.estimation.algos.include_delightPZ.convertDESCcat import convertDESCcat # convert DESC input file into Delight format +from rail.estimation.algos.include_delightPZ.convertDESCcat import * # convert DESC input file into Delight format + +# mock data simulation +#from interfaces.rail.simulateWithSEDs import simulateWithSEDs # simulate its own SED in tutorial mode +from rail.estimation.algos.include_delightPZ.simulateWithSEDs import simulateWithSEDs # simulate its own SED in tutorial mode + +# Delight algorithms + +#from interfaces.rail.templateFitting import templateFitting +from rail.estimation.algos.include_delightPZ.templateFitting import templateFitting +#from interfaces.rail.delightLearn import delightLearn +from rail.estimation.algos.include_delightPZ.delightLearn import delightLearn +#from interfaces.rail.delightApply import delightApply +from rail.estimation.algos.include_delightPZ.delightApply import delightApply + +# other + +#from interfaces.rail.calibrateTemplateMixturePriors import * +from rail.estimation.algos.include_delightPZ.calibrateTemplateMixturePriors import * +#from interfaces.rail.getDelightRedshiftEstimation import * +from rail.estimation.algos.include_delightPZ.getDelightRedshiftEstimation import * + +# Create a logger object. +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s %(name)s[%(process)d] %(levelname)s %(message)s') + + +class delightPZ(BaseEstimation): + + def __init__(self, base_config, config_dict): + """ + Parameters: + ----------- + run_dict: dict + dictionary of all variables read in from the run_params + values in the yaml file + """ + super().__init__(base_config=base_config, config_dict=config_dict) + + inputs = self.config_dict['run_params'] + + self.width = inputs['dlght_redshiftBinSize'] + self.zmin = inputs['dlght_redshiftMin'] + self.zmax = inputs['dlght_redshiftMax'] + self.nzbins = int((self.zmax-self.zmin)/self.width) + + # temporary directories for Delight temprary file + self.tempdir = inputs['tempdir'] + self.tempdatadir = inputs['tempdatadir'] + # name of delight configuration file + self.delightparamfile=inputs["delightparamfile"] + self.delightparamfile = os.path.join(self.tempdir, self.delightparamfile) + + self.delightindata=inputs['dlght_inputdata'] + + # Choice of the running mode + self.tutorialmode = inputs["dlght_tutorialmode"] + self.tutorialpasseval = False # onmy one chunk for simulation + # for standard mode with DC2 dataset + self.flag_filter_training = inputs["flag_filter_training"] + self.snr_cut_training = inputs["snr_cut_training"] + self.flag_filter_validation = inputs["flag_filter_validation"] + self.snr_cut_validation = inputs["snr_cut_validation"] + + + # counter on the chunk validation dataset + self.chunknum=0 + + self.dlght_calibrateTemplateMixturePrior = inputs["dlght_calibrateTemplateMixturePrior"] + # all parameter files + self.inputs = inputs + + np.random.seed(87) + + ############################################################################################################# + # + # INFORM + # + ############################################################################################################ + + def inform(self): + """ + this is random, so does nothing + """ + + msg = " INFORM " + logger.info(msg) + + + logger.info("Try to workout filters") + + # create usefull tempory directory + try: + if not os.path.exists(self.tempdir): + os.makedirs(self.tempdir) + except OSError as e: + if e.errno != errno.EEXIST: + msg = "error creating file "+self.tempdir + logger.error(msg) + raise + + try: + if not os.path.exists(self.tempdatadir): + os.makedirs(self.tempdatadir) + except OSError as e: + if e.errno != errno.EEXIST: + msg = "error creating file " + self.tempdatadir + logger.error(msg) + raise + + + if not os.path.exists(self.delightindata): + msg = " No Delight input data in dir " + self.delightindata + logger.error(msg) + exit(-1) + + SUBDIRs = ['BROWN_SEDs', 'CWW_SEDs', 'FILTERS'] + + for subdir in SUBDIRs: + theinpath=os.path.join(self.delightindata,subdir) + if not os.path.exists(theinpath): + msg = " No Delight input data in dir " + theinpath + logger.error(msg) + exit(-1) + + + + + + + + + # Very awfull way to get the internal data path for Delight + # This is the delight setup.py tools that get the data there + #basedelight_datapath = resource_filename('delight', '../data') + basedelight_datapath = self.delightindata + + logger.debug("Guessed basedelight_datapath = " + basedelight_datapath ) + + # gives to Delight where are its datapath and provide yaml arguments + paramfile_txt=makeConfigParam(basedelight_datapath, self.inputs) + logger.debug(paramfile_txt) + + + # save the config parameter file that Delight needs + with open(self.delightparamfile, 'w') as out: + out.write(paramfile_txt) + + + # Later will steer the calls to the Delight functions with RAIL config file + + # Build LSST filter model + processFilters(self.delightparamfile) + + # Build its own LSST-Flux-Redshift Model + processSEDs(self.delightparamfile) + + if self.tutorialmode: + # Delight build its own Mock simulations + simulateWithSEDs(self.delightparamfile) + + else: # convert hdf5 into ascii in desc input mode + convertDESCcat(self.delightparamfile, self.trainfile, self.testfile,\ + flag_filter_training=self.flag_filter_training,\ + flag_filter_validation=self.flag_filter_validation,\ + snr_cut_training=self.snr_cut_training,\ + snr_cut_validation=self.snr_cut_validation) + + + if self.dlght_calibrateTemplateMixturePrior: + calibrateTemplateMixturePriors(self.delightparamfile) + + + + + # Learn with Gaussian processes + delightLearn(self.delightparamfile) + + + ############################################################################################################# + # + # ESTIMATE + # + ############################################################################################################ + + + def estimate(self, test_data): + + self.chunknum += 1 + + msg = " ESTIMATE : chunk number {} ".format(self.chunknum) + logger.info(msg) + + #basedelight_datapath = resource_filename('delight', '../data') + basedelight_datapath = self.delightindata + + msg=" Delight input data file are in dir : {} ".format(self.delightparamfile) + logger.debug(msg) + + # when Delight runs in tutorial mode call only once delightApply + if self.tutorialmode and not self.tutorialpasseval: + + msg = "TUTORIAL MODE : process chunk {}".format(self.chunknum) + logger.info(msg) + + # Template Fitting + templateFitting(self.delightparamfile) + + # Gaussian process fitting + delightApply(self.delightparamfile) + self.tutorialpasseval = True # avoid latter call to delightApply when running in tutorial mode + + elif self.tutorialmode and self.tutorialpasseval: + msg="TUTORIAL MODE : skip chunk {}".format(self.chunknum) + logger.info(msg) + + elif not self.tutorialmode : # let rail split the test data into chunks + msg = "STANDARD MODE : process chunk {}".format(self.chunknum) + logger.info(msg) + + # Generate a new parameter file for delight this chunk + paramfile_txt=makeConfigParamChunk(basedelight_datapath, self.inputs, self.chunknum) + + # generate the config-parameter filename from chunk number + delightparamfile=self.delightparamfile + logger.debug(delightparamfile) + dirn=os.path.dirname(delightparamfile) + basn=os.path.basename(delightparamfile) + basnsplit=basn.split(".") + basnchunk = basnsplit[0] + "_" + str(self.chunknum) + "." + basnsplit[1] + delightparamfilechunk = os.path.join(dirn,basnchunk) + logger.debug("parameter file for delight :" + delightparamfilechunk) + + # save the config parameter file for the data chunk that Delight needs + with open(delightparamfilechunk, 'w') as out: + out.write(paramfile_txt) + + # convert the chunk data into the required flux-redshift validation file for delight + indexes_sel=convertDESCcatChunk(delightparamfilechunk, test_data, self.chunknum,\ + flag_filter_validation=self.flag_filter_validation,\ + snr_cut_validation=self.snr_cut_validation) + + # template fitting for that chunk + templateFitting(delightparamfilechunk) + + # estimation for that chunk + delightApply(delightparamfilechunk) + + + + else: + msg = "STANDARD MODE : pass chunk {}".format(self.chunknum) + logger.info(msg) + + + pdf = [] + # allow for either format for now + try: + d = test_data['i_mag'] + except Exception: + d = test_data['mag_i_lsst'] + + numzs = len(d) + self.zgrid = np.linspace(self.zmin, self.zmax, self.nzbins) + + + if self.tutorialmode: + # fill creazy values (simulation not send to rail) + zmode = np.round(np.random.uniform(self.zmax, self.zmax, numzs), 3) + widths = self.width * (1.0 + zmode) + + else: + zmode,widths = getDelightRedshiftEstimation(delightparamfilechunk,self.chunknum,numzs,indexes_sel) + zmode = np.round(zmode,3) + + + for i in range(numzs): + pdf.append(norm.pdf(self.zgrid, zmode[i], widths[i])) + + pz_dict = {'zmode': zmode, 'pz_pdf': pdf} + + return pz_dict diff --git a/rail/estimation/algos/include_delightPZ/__init__.py b/rail/estimation/algos/include_delightPZ/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/rail/estimation/algos/include_delightPZ/calibrateTemplateMixturePriors.py b/rail/estimation/algos/include_delightPZ/calibrateTemplateMixturePriors.py new file mode 100644 index 00000000..881daeb8 --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/calibrateTemplateMixturePriors.py @@ -0,0 +1,267 @@ + +import sys +from mpi4py import MPI +import numpy as np +from scipy.interpolate import interp1d +from delight.io import * +from delight.utils import * +from delight.photoz_gp import PhotozGP +from delight.photoz_kernels import Photoz_mean_function, Photoz_kernel +import scipy.stats +import matplotlib.pyplot as plt +import emcee +import corner + +import coloredlogs +import logging + + +# Create a logger object. +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s %(name)s[%(process)d] %(levelname)s %(message)s') + + + + + +def calibrateTemplateMixturePriors(configfilename,make_plot=False): + """ + + :param configfilename: + :return: + """ + # Parse parameters file + + logger.info("--- calibrate Template Mixture Priors ---") + + paramFileName = configfilename + params = parseParamFile(paramFileName, verbose=False) + + comm = MPI.COMM_WORLD + threadNum = comm.Get_rank() + numThreads = comm.Get_size() + + DL = approx_DL() + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + numZ = redshiftGrid.size + + # Locate which columns of the catalog correspond to which bands. + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn,refBandColumn = readColumnPositions(params, prefix="training_") + + dir_seds = params['templates_directory'] + dir_filters = params['bands_directory'] + lambdaRef = params['lambdaRef'] + sed_names = params['templates_names'] + numBands = bandIndices.size + nt = len(sed_names) + + + f_mod = np.zeros((numZ, nt, len(params['bandNames']))) + + # model of flux-redshift for each template + for t, sed_name in enumerate(sed_names): + f_mod[:, t, :] = np.loadtxt(dir_seds + '/' + sed_name + '_fluxredshiftmod.txt') + + numObjectsTraining = np.sum(1 for line in open(params['training_catFile'])) + + msg = 'Number of Training Objects ' + str(numObjectsTraining) + logger.info(msg) + + numMetrics = 7 + len(params['confidenceLevels']) + + + allFluxes = np.zeros((numObjectsTraining, numBands)) + allFluxesVar = np.zeros((numObjectsTraining, numBands)) + + redshifts = np.zeros((numObjectsTraining, 1)) + fmod_atZ = np.zeros((numObjectsTraining, nt, numBands)) + + # Now loop over training set to compute likelihood function + loc = - 1 + trainingDataIter = getDataFromFile(params, 0, numObjectsTraining,prefix="training_", getXY=False) + + # loop on traning + for z, ell, bands, fluxes, fluxesVar, bCV, fCV, fvCV in trainingDataIter: + loc += 1 + allFluxes[loc, :] = fluxes + allFluxesVar[loc, :] = fluxesVar + redshifts[loc, 0] = z + + # loop on SED + for t, sed_name in enumerate(sed_names): + for ib, b in enumerate(bands): + fmod_atZ[loc, t, ib] = ell * np.interp(z, redshiftGrid,f_mod[:, t, b]) + + zZmax = redshifts[:, 0] / redshiftGrid[-1] + + + pmin = np.repeat(0., 2*nt) + pmax = np.repeat(200., 2*nt) + + ndim, nwalkers = 2*nt, 100 + + def approx_flux_likelihood_multiobj( + f_obs, # no, nf + f_obs_var, # no, nf + f_mod, # no, nt, nf + ell_hat, # 1 + ell_var, # 1 + returnChi2=False, + normalized=True): + + assert len(f_obs.shape) == 2 + assert len(f_obs_var.shape) == 2 + assert len(f_mod.shape) == 3 + no, nt, nf = f_mod.shape + f_obs_r = f_obs[:, None, :] + var = f_obs_var[:, None, :] + invvar = np.where(f_obs_r / var < 1e-6, 0.0, var ** -1.0) # nz * nt * nf + FOT = np.sum(f_mod * f_obs_r * invvar, axis=2) \ + + ell_hat / ell_var # no * nt + FTT = np.sum(f_mod ** 2 * invvar, axis=2) \ + + 1. / ell_var # no * nt + FOO = np.sum(f_obs_r ** 2 * invvar, axis=2) \ + + ell_hat ** 2 / ell_var # no * nt + sigma_det = np.prod(var, axis=2) + chi2 = FOO - FOT ** 2.0 / FTT # no * nt + denom = np.sqrt(FTT) + if normalized: + denom *= np.sqrt(sigma_det * (2 * np.pi) ** nf * ell_var) + like = np.exp(-0.5 * chi2) / denom # no * nt + if returnChi2: + return chi2 + else: + return like + + def lnprob(params, nt, allFluxes, allFluxesVar, zZmax, fmod_atZ, pmin, pmax): + if np.any(params > pmax) or np.any(params < pmin): + return - np.inf + + alphas0 = dirichlet(params[0:nt], rsize=1).ravel()[None, :] # 1, nt + alphas1 = dirichlet(params[nt:2 * nt], rsize=1).ravel()[None, :] # 1, nt + alphas_atZ = zZmax[:, None] * (alphas1 - alphas0) + alphas0 # no, nt + # fmod_atZ: no, nt, nf + fmod_atZ_t = (fmod_atZ * alphas_atZ[:, :, None]).sum(axis=1)[:, None, :] + # no, 1, nf + sigma_ell = 1e3 + like_grid = approx_flux_likelihood_multiobj(allFluxes, allFluxesVar, fmod_atZ_t, 1, + sigma_ell ** 2.).ravel() # no, + eps = 1e-305 + ind = like_grid > eps + theprob = np.log(like_grid[ind]).sum() + return theprob + + + p0 = [pmin + (pmax-pmin)*np.random.uniform(0, 1, size=ndim) for i in range(nwalkers)] + + for i in range(10): + print(lnprob(p0[i], nt, allFluxes, allFluxesVar, zZmax, fmod_atZ, pmin, pmax)) + + sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob,threads=4,args=[nt, allFluxes, allFluxesVar, zZmax,fmod_atZ, pmin, pmax]) + pos, prob, state = sampler.run_mcmc(p0, 200) + sampler.reset() + sampler.run_mcmc(pos, 2000) + print("Mean acceptance fraction: {0:.3f}" .format(np.mean(sampler.acceptance_fraction))) + samples = sampler.chain.reshape((-1, ndim)) + lnprob = sampler.lnprobability.reshape((-1, 1)) + + params_mean = samples.mean(axis=0) + params_std = samples.std(axis=0) + + if make_plot: + fig, axs = plt.subplots(4, 5, figsize=(16, 8)) + axs = axs.ravel() + for i in range(ndim): + axs[i].hist(samples[:, i], 50, color="k", histtype="step") + axs[i].axvspan(params_mean[i]-params_std[i],params_mean[i]+params_std[i], color='gray', alpha=0.5) + axs[i].axvline(params_mean[i], c='k', lw=2) + fig.tight_layout() + fig.savefig('priormixture_parameters.pdf') + + fig = plot_params(params_mean) + fig.savefig('priormixture_meanparameters.pdf') + + print("params_mean", params_mean) + print("params_std", params_std) + + alphas = params_mean[0:nt] + alpha0 = np.sum(alphas) + print("alpha0:", ', '.join(['%.2g' % x for x in alphas / alpha0])) + print("alpha0 err:", ', '.join(['%.2g' % x for x in np.sqrt(alphas*(alpha0-alphas)/alpha0**2/(alpha0+1))])) + + alphas = params_mean[nt:2*nt] + alpha0 = np.sum(alphas) + print("alpha1:", ', '.join(['%.2g' % x for x in alphas / alpha0])) + print("alpha1 err:", ', '.join(['%.2g' % x for x in np.sqrt(alphas*(alpha0-alphas)/alpha0**2/(alpha0+1))])) + + + + alphas = params_mean[0:nt] + betas = params_mean[nt:2 * nt] + + alpha0 = np.sum(alphas) + print("p_t:", ' '.join(['%.2g' % x for x in alphas / alpha0])) + print("p_z_t:", ' '.join(['%.2g' % x for x in betas])) + print("p_t err:", ' '.join(['%.2g' % x for x in np.sqrt(alphas * (alpha0 - alphas) / alpha0 ** 2 / (alpha0 + 1))])) + + + + + if make_plot: + fig = corner.corner(samples) + fig.savefig("trianglemixture.pdf") + + +def plot_params(params): + fig, axs = plt.subplots(4, 4, figsize=(16, 8)) + axs = axs.ravel() + alphas = params[0:nt] + alpha0 = np.sum(alphas) + dirsamples = dirichlet(alphas, 1000) + for i in range(nt): + mean = alphas[i]/alpha0 + std = np.sqrt(alphas[i] * (alpha0-alphas[i]) / alpha0**2 / (alpha0+1)) + axs[i].axvspan(mean-std, mean+std, color='gray', alpha=0.5) + axs[i].axvline(mean, c='k', lw=2) + axs[i].axvline(1/nt, c='k', lw=1, ls='dashed') + axs[i].set_title('alpha0 = '+str(alphas[i])) + axs[i].set_xlim([0, 1]) + axs[i].hist(dirsamples[:, i], 50, color="k", histtype="step") + alphas = params[nt:2*nt] + alpha0 = np.sum(alphas) + dirsamples = dirichlet(alphas, 1000) + for i in range(nt): + mean = alphas[i]/alpha0 + std = np.sqrt(alphas[i] * (alpha0-alphas[i]) / alpha0**2 / (alpha0+1)) + axs[nt+i].axvspan(mean-std, mean+std, color='gray', alpha=0.5) + axs[nt+i].axvline(mean, c='k', lw=2) + axs[nt+i].axvline(1/nt, c='k', lw=1, ls='dashed') + axs[nt+i].set_title('alpha1 = '+str(alphas[i])) + axs[nt+i].set_xlim([0, 1]) + axs[nt+i].hist(dirsamples[:, i], 50, color="k", histtype="step") + fig.tight_layout() + return fig + +#----------------------------------------------------------------------------------------- +if __name__ == "__main__": + # execute only if run as a script + + + msg=" calibrateTemplateMixturePriors " + logger.info(msg) + + + logger.debug("__name__:"+__name__) + logger.debug("__file__:"+__file__) + + logger.info(" calibrate Template Priors ---") + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + calibrateTemplateMixturePriors(sys.argv[1]) + + + + + diff --git a/rail/estimation/algos/include_delightPZ/convertDESCcat.py b/rail/estimation/algos/include_delightPZ/convertDESCcat.py new file mode 100644 index 00000000..59c8040e --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/convertDESCcat.py @@ -0,0 +1,646 @@ +####################################################################################################### +# +# script : convertDESCcat.py +# +# convert DESC catalog to be injected in Delight +# produce files `galaxies-redshiftpdfs.txt` and `galaxies-redshiftpdfs2.txt` for training and target +# +######################################################################################################### + + +import sys +import os +import numpy as np +from functools import reduce + +import pprint + +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +from delight.io import * +from delight.utils import * + +import coloredlogs +import logging + +import h5py + +#from interfaces.rail.utils import load_training_data, get_input_data_size_hdf5,load_raw_hdf5_data +from rail.estimation.utils import load_training_data, get_input_data_size_hdf5,load_raw_hdf5_data + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + +# option to convert DC2 flux level (in AB units) into internal Delight units +# this option will be removed when optimisation of parameters will be implemented +FLAG_CONVERTFLUX_TODELIGHTUNIT=True + + + +def group_entries(f): + """ + group entries in single numpy array + + """ + galid = f['id'][()][:, np.newaxis] + redshift = f['redshift'][()][:, np.newaxis] + mag_err_g_lsst = f['mag_err_g_lsst'][()][:, np.newaxis] + mag_err_i_lsst = f['mag_err_i_lsst'][()][:, np.newaxis] + mag_err_r_lsst = f['mag_err_r_lsst'][()][:, np.newaxis] + mag_err_u_lsst = f['mag_err_u_lsst'][()][:, np.newaxis] + mag_err_y_lsst = f['mag_err_y_lsst'][()][:, np.newaxis] + mag_err_z_lsst = f['mag_err_z_lsst'][()][:, np.newaxis] + mag_g_lsst = f['mag_g_lsst'][()][:, np.newaxis] + mag_i_lsst = f['mag_i_lsst'][()][:, np.newaxis] + mag_r_lsst = f['mag_r_lsst'][()][:, np.newaxis] + mag_u_lsst = f['mag_u_lsst'][()][:, np.newaxis] + mag_y_lsst = f['mag_y_lsst'][()][:, np.newaxis] + mag_z_lsst = f['mag_z_lsst'][()][:, np.newaxis] + + full_arr = np.hstack((galid, redshift, mag_u_lsst, mag_g_lsst, mag_r_lsst, mag_i_lsst, mag_z_lsst, mag_y_lsst, \ + mag_err_u_lsst, mag_err_g_lsst, mag_err_r_lsst, mag_err_i_lsst, mag_err_z_lsst, + mag_err_y_lsst)) + return full_arr + + +def filter_mag_entries(d,nb=6): + """ + Filter bad data with bad magnitudes + + input + - d: array of magnitudes and errors + - nb : number of bands + output : + - indexes of row to be filtered + + """ + + u = d[:, 2] + idx_u = np.where(u > 31.8)[0] + + return idx_u + + +def mag_to_flux(d,nb=6): + """ + + Convert magnitudes to fluxes + + input: + -d : array of magnitudes with errors + + + :return: + array of fluxes with error + """ + + fluxes = np.zeros_like(d) + + fluxes[:, 0] = d[:, 0] # object index + fluxes[:, 1] = d[:, 1] # redshift + + for idx in np.arange(nb): + fluxes[:, 2 + idx] = np.power(10, -0.4 * d[:, 2 + idx]) # fluxes + fluxes[:, 8 + idx] = fluxes[:, 2 + idx] * d[:, 8 + idx] # errors on fluxes + return fluxes + + + +def filter_flux_entries(d,nb=6,nsig=5): + """ + Filter noisy data on the the number SNR + + input : + - d: flux and errors array + - nb : number of bands + - nsig : number of sigma + + output: + indexes of row to suppress + + """ + + + # collection of indexes + indexes = [] + #indexes = np.array(indexes, dtype=np.int) + indexes = np.array(indexes, dtype=int) + + for idx in np.arange(nb): + ratio = d[:, 2 + idx] / d[:, 8 + idx] # flux divided by sigma-flux + bad_indexes = np.where(ratio < nsig)[0] + indexes = np.concatenate((indexes, bad_indexes)) + + indexes = np.unique(indexes) + return np.sort(indexes) + + +def convertDESCcatChunk(configfilename,data,chunknum,flag_filter_validation = True, snr_cut_validation = 5): + + """ + convertDESCcatChunk(configfilename,data,chunknum,flag_filter_validation = True, snr_cut_validation = 5) + + Convert files in ascii format to be used by Delight + Input data can be filtered by series of filters. But it is necessary to remember which entries are kept, + which are eliminated + + input args: + - configfilename : Delight configuration file containing path for output files (flux variances and redshifts) + - data : the DC2 data + - chunknum : number of the chunk + - filter_validation : Flag to activate quality filter data + - snr_cut_validation : cut on flux SNR + + output : + - the target file of the chunk which path is in configuration file + :return: + - the list of selected (unfiltered DC2 data) + """ + msg="--- Convert DESC catalogs chunk {}---".format(chunknum) + logger.info(msg) + + if FLAG_CONVERTFLUX_TODELIGHTUNIT: + flux_multiplicative_factor = 2.22e10 + else: + flux_multiplicative_factor = 1 + + + + # produce a numpy array + magdata = group_entries(data) + + + # remember the number of entries + Nin = magdata.shape[0] + msg = "Number of objects = {} , in chunk : {}".format(Nin,chunknum) + logger.debug(msg) + + + # keep indexes to filter data with bad magnitudes + if flag_filter_validation: + indexes_bad_mag = filter_mag_entries(magdata) + #magdata_f = np.delete(magdata, indexes_bad_mag, axis=0) + magdata_f = magdata # filtering will be done later + + + else: + indexes_bad_mag=np.array([]) + magdata_f = magdata + + Nbadmag = len(indexes_bad_mag) + msg = "Number of objects with bad magnitudes = {} , in chunk : {}".format(Nbadmag, chunknum) + logger.debug(msg) + + #print("indexes_bad_mag = ",indexes_bad_mag) + + + # convert mag to fluxes + fdata = mag_to_flux(magdata_f) + + # keep indexes to filter data with bad SNR + if flag_filter_validation: + indexes_bad_snr = filter_flux_entries(fdata, nsig = snr_cut_validation) + fdata_f = fdata + #fdata_f = np.delete(fdata, indexes_bad, axis=0) + #magdata_f = np.delete(magdata_f, indexes_bad, axis=0) + else: + fdata_f=fdata + indexes_bad_snr = np.array([]) + + + Nbadsnr = len(indexes_bad_snr) + msg = "Number of objects with bad SNR = {} , in chunk : {}".format(Nbadsnr, chunknum) + logger.debug(msg) + + #print("indexes_bad_snr = ", indexes_bad_snr) + + # make union of indexes (unique id) before removing them for Delight + idxToRemove = reduce(np.union1d,(indexes_bad_mag,indexes_bad_snr)) + NtoRemove=len(idxToRemove) + msg = "Number of objects filtered out = {} , in chunk : {}".format(NtoRemove, chunknum) + logger.debug(msg) + + #print("indexes_to_remove = ", idxToRemove) + + #pprint(idxToRemove) + + # fdata_f contains the fluxes and errors to be send to Delight + + # indexes of full input dataset + idxInitial = np.arange(Nin) + + if NtoRemove>0: + fdata_f = np.delete(fdata_f,idxToRemove, axis=0) + idxFinal=np.delete(idxInitial,idxToRemove, axis=0) + else: + idxFinal = idxInitial + + + Nkept = len(idxFinal) + msg = "Number of objects kept = {} , in chunk : {}".format(Nkept, chunknum) + logger.debug(msg) + + #print("indexes_kept = ", idxFinal) + + + + gid = fdata_f[:, 0] + rs = fdata_f[:, 1] + + # 2) parameter file + + params = parseParamFile(configfilename, verbose=False, catFilesNeeded=False) + + numB = len(params['bandNames']) + numObjects = len(gid) + + msg = "get {} objects ".format(numObjects) + logger.debug(msg) + + logger.debug(params['bandNames']) + + # Generate target data + # ------------------------- + + # what is fluxes and fluxes variance + fluxes, fluxesVar = np.zeros((numObjects, numB)), np.zeros((numObjects, numB)) + + # loop on objects to simulate for the target and save in output trarget file + for k in range(numObjects): + # loop on number of bands + for i in range(numB): + trueFlux = fdata_f[k, 2 + i] + noise = fdata_f[k, 8 + i] + + # put the DC2 data to the internal units of Delight + trueFlux *= flux_multiplicative_factor + noise *= flux_multiplicative_factor + + + # fluxes[k, i] = trueFlux + noise * np.random.randn() # noisy flux + fluxes[k, i] = trueFlux + + if fluxes[k, i] < 0: + # fluxes[k, i]=np.abs(noise)/10. + fluxes[k, i] = trueFlux + + fluxesVar[k, i] = noise ** 2. + + # container for target galaxies output + # at some redshift, provides the flux and its variance inside each band + + + data = np.zeros((numObjects, 1 + len(params['target_bandOrder']))) + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn, refBandColumn = readColumnPositions(params, + prefix="target_") + + for ib, pf, pfv in zip(bandIndices, bandColumns, bandVarColumns): + data[:, pf] = fluxes[:, ib] + data[:, pfv] = fluxesVar[:, ib] + data[:, redshiftColumn] = rs + data[:, -1] = 0 # NO TYPE + + msg = "write file {}".format(os.path.basename(params['targetFile'])) + logger.debug(msg) + + msg = "write target file {}".format(params['targetFile']) + logger.debug(msg) + + outputdir = os.path.dirname(params['targetFile']) + if not os.path.exists(outputdir): + msg = " outputdir not existing {} then create it ".format(outputdir) + logger.info(msg) + os.makedirs(outputdir) + + np.savetxt(params['targetFile'], data) + + # return the index of selected data + return idxFinal + + + +def convertDESCcat(configfilename,desctraincatalogfile,desctargetcatalogfile,\ + flag_filter_training=True,flag_filter_validation=True,snr_cut_training=5,snr_cut_validation=5): + + """ + convertDESCcat(configfilename,desctraincatalogfile,desctargetcatalogfile,\ + flag_filter_training=True,flag_filter_validation=True,snr_cut_training=5,snr_cut_validation=5): + + + Convert files in ascii format to be used by Delight + + input args: + - configfilename : Delight configuration file containingg path for output files (flux variances and redshifts) + - desctraincatalogfile : training file provided by RAIL (hdf5 format) + - desctargetcatalogfile : target file provided by RAIL (hdf5 format) + - flag_filter_training : Activate filtering on training data + - flag_filter_validation : Activate filtering on validation data + - snr_cut_training : Cut on flux SNR in training data + - snr_cut_validation : Cut on flux SNR in validation data + + output : + - the Delight training and target file which path is in configuration file + + :return: nothing + + """ + + + logger.info("--- Convert DESC training and target catalogs ---") + + if FLAG_CONVERTFLUX_TODELIGHTUNIT: + flux_multiplicative_factor = 2.22e10 + else: + flux_multiplicative_factor = 1 + + + + # 1) DESC catalog file + msg="read DESC hdf5 training file {} ".format(desctraincatalogfile) + logger.debug(msg) + + f = load_raw_hdf5_data(desctraincatalogfile, groupname='photometry') + + # produce a numpy array + magdata = group_entries(f) + + # remember the number of entries + Nin = magdata.shape[0] + msg = "Number of objects = {} , in training dataset".format(Nin) + logger.debug(msg) + + + + # keep indexes to filter data with bad magnitudes + if flag_filter_training: + indexes_bad_mag = filter_mag_entries(magdata) + # magdata_f = np.delete(magdata, indexes_bad_mag, axis=0) + magdata_f = magdata # filtering will be done later + else: + indexes_bad_mag = np.array([]) + magdata_f = magdata + + Nbadmag = len(indexes_bad_mag) + msg = "Number of objects with bad magnitudes {} in training dataset".format(Nbadmag) + logger.debug(msg) + + + # convert mag to fluxes + fdata = mag_to_flux(magdata_f) + + # keep indexes to filter data with bad SNR + if flag_filter_training: + indexes_bad_snr = filter_flux_entries(fdata, nsig=snr_cut_training) + fdata_f = fdata + # fdata_f = np.delete(fdata, indexes_bad, axis=0) + # magdata_f = np.delete(magdata_f, indexes_bad, axis=0) + else: + fdata_f = fdata + indexes_bad_snr = np.array([]) + + Nbadsnr = len(indexes_bad_snr) + msg = "Number of objects with bad SNR = {} , in training dataset".format(Nbadsnr) + logger.debug(msg) + + # make union of indexes (unique id) before removing them for Delight + idxToRemove = reduce(np.union1d, (indexes_bad_mag, indexes_bad_snr)) + NtoRemove = len(idxToRemove) + msg = "Number of objects filtered out = {} , in training dataset".format(NtoRemove) + logger.debug(msg) + + + # fdata_f contains the fluxes and errors to be send to Delight + + # indexes of full input dataset + idxInitial = np.arange(Nin) + + if NtoRemove > 0: + fdata_f = np.delete(fdata_f, idxToRemove, axis=0) + idxFinal = np.delete(idxInitial, idxToRemove, axis=0) + else: + idxFinal = idxInitial + + + Nkept = len(idxFinal) + msg = "Number of objects kept = {} , in training dataset".format(Nkept) + logger.debug(msg) + + + + gid = fdata_f[:, 0] + rs = fdata_f[:, 1] + + + # 2) parameter file + + params = parseParamFile(configfilename, verbose=False, catFilesNeeded=False) + + numB = len(params['bandNames']) + numObjects = len(gid) + + msg = "get {} objects ".format(numObjects) + logger.debug(msg) + + logger.debug(params['bandNames']) + + + + # Generate training data + #------------------------- + + + # what is fluxes and fluxes variance + fluxes, fluxesVar = np.zeros((numObjects, numB)), np.zeros((numObjects, numB)) + + # loop on objects to simulate for the training and save in output training file + for k in range(numObjects): + #loop on number of bands + for i in range(numB): + trueFlux = fdata_f[k,2+i] + noise = fdata_f[k,8+i] + + # put the DC2 data to the internal units of Delight + trueFlux *= flux_multiplicative_factor + noise *= flux_multiplicative_factor + + + #fluxes[k, i] = trueFlux + noise * np.random.randn() # noisy flux + fluxes[k, i] = trueFlux + + if fluxes[k, i]<0: + #fluxes[k, i]=np.abs(noise)/10. + fluxes[k, i] = trueFlux + + fluxesVar[k, i] = noise**2. + + # container for training galaxies output + # at some redshift, provides the flux and its variance inside each band + data = np.zeros((numObjects, 1 + len(params['training_bandOrder']))) + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn,refBandColumn = readColumnPositions(params, prefix="training_") + + for ib, pf, pfv in zip(bandIndices, bandColumns, bandVarColumns): + data[:, pf] = fluxes[:, ib] + data[:, pfv] = fluxesVar[:, ib] + data[:, redshiftColumn] = rs + data[:, -1] = 0 # NO type + + + msg="write training file {}".format(params['trainingFile']) + logger.debug(msg) + + outputdir=os.path.dirname(params['trainingFile']) + if not os.path.exists(outputdir): + msg = " outputdir not existing {} then create it ".format(outputdir) + logger.info(msg) + os.makedirs(outputdir) + + + np.savetxt(params['trainingFile'], data) + + + + + # Generate Target data : procedure similar to the training + #----------------------------------------------------------- + + # 1) DESC catalog file + msg = "read DESC hdf5 validation file {} ".format(desctargetcatalogfile) + logger.debug(msg) + + f = load_raw_hdf5_data(desctargetcatalogfile, groupname='photometry') + + # produce a numpy array + magdata = group_entries(f) + + + # remember the number of entries + Nin = magdata.shape[0] + msg = "Number of objects = {} , in validation dataset".format(Nin) + logger.debug(msg) + + + # filter bad data + # keep indexes to filter data with bad magnitudes + if flag_filter_validation: + indexes_bad_mag = filter_mag_entries(magdata) + # magdata_f = np.delete(magdata, indexes_bad_mag, axis=0) + magdata_f = magdata # filtering will be done later + else: + indexes_bad_mag = np.array([]) + magdata_f = magdata + + Nbadmag = len(indexes_bad_mag) + msg = "Number of objects with bad magnitudes = {} , in validation dataset".format(Nbadmag) + logger.debug(msg) + + + + # convert mag to fluxes + fdata = mag_to_flux(magdata_f) + + # keep indexes to filter data with bad SNR + if flag_filter_validation: + indexes_bad_snr = filter_flux_entries(fdata, nsig=snr_cut_validation) + fdata_f = fdata + # fdata_f = np.delete(fdata, indexes_bad, axis=0) + # magdata_f = np.delete(magdata_f, indexes_bad, axis=0) + else: + fdata_f = fdata + indexes_bad_snr = np.array([]) + + Nbadsnr = len(indexes_bad_snr) + msg = "Number of objects with bad SNR = {} , in validation dataset".format(Nbadsnr) + logger.debug(msg) + + # make union of indexes (unique id) before removing them for Delight + idxToRemove = reduce(np.union1d, (indexes_bad_mag, indexes_bad_snr)) + NtoRemove = len(idxToRemove) + msg = "Number of objects filtered out = {} , in validation dataset".format(NtoRemove) + logger.debug(msg) + + # fdata_f contains the fluxes and errors to be send to Delight + + # indexes of full input dataset + idxInitial = np.arange(Nin) + + if NtoRemove > 0: + fdata_f = np.delete(fdata_f, idxToRemove, axis=0) + idxFinal = np.delete(idxInitial, idxToRemove, axis=0) + else: + idxFinal = idxInitial + + + Nkept = len(idxFinal) + msg = "Number of objects kept = {} , in validation dataset".format(Nkept) + logger.debug(msg) + + gid = fdata_f[:, 0] + rs = fdata_f[:, 1] + + numObjects = len(gid) + msg = "get {} objects ".format(numObjects) + logger.debug(msg) + + fluxes, fluxesVar = np.zeros((numObjects, numB)), np.zeros((numObjects, numB)) + + # loop on objects in target files + for k in range(numObjects): + # loop on bands + for i in range(numB): + # compute the flux in that band at the redshift + trueFlux = fdata_f[k, 2 + i] + noise = fdata_f[k, 8 + i] + + # put the DC2 data to the internal units of Delight + trueFlux *= flux_multiplicative_factor + noise *= flux_multiplicative_factor + + #fluxes[k, i] = trueFlux + noise * np.random.randn() + fluxes[k, i] = trueFlux + + if fluxes[k, i]<0: + #fluxes[k, i]=np.abs(noise)/10. + fluxes[k, i] = trueFlux + + fluxesVar[k, i] = noise**2 + + + + data = np.zeros((numObjects, 1 + len(params['target_bandOrder']))) + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn,refBandColumn = readColumnPositions(params, prefix="target_") + + for ib, pf, pfv in zip(bandIndices, bandColumns, bandVarColumns): + data[:, pf] = fluxes[:, ib] + data[:, pfv] = fluxesVar[:, ib] + data[:, redshiftColumn] = rs + data[:, -1] = 0 # NO TYPE + + msg = "write file {}".format(os.path.basename(params['targetFile'])) + logger.debug(msg) + + msg = "write target file {}".format(params['targetFile']) + logger.debug(msg) + + outputdir = os.path.dirname(params['targetFile']) + if not os.path.exists(outputdir): + msg = " outputdir not existing {} then create it ".format(outputdir) + logger.info(msg) + os.makedirs(outputdir) + + np.savetxt(params['targetFile'], data) + + +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start convertDESCcat.py" + logger.info(msg) + logger.info("--- convert DESC catalogs ---") + + + + if len(sys.argv) < 4: + raise Exception('Please provide a parameter file and the training and validation and catalog files') + + convertDESCcat(sys.argv[1],sys.argv[2],sys.argv[3]) diff --git a/rail/estimation/algos/include_delightPZ/delightApply.py b/rail/estimation/algos/include_delightPZ/delightApply.py new file mode 100644 index 00000000..723c3de0 --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/delightApply.py @@ -0,0 +1,262 @@ + +import sys +from mpi4py import MPI +import numpy as np +from delight.io import * +from delight.utils import * +from delight.photoz_gp import PhotozGP +from delight.photoz_kernels import Photoz_mean_function, Photoz_kernel +from delight.utils_cy import approx_flux_likelihood_cy +from time import time + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + + + +def delightApply(configfilename): + """ + + :param configfilename: + :return: + """ + + + comm = MPI.COMM_WORLD + threadNum = comm.Get_rank() + numThreads = comm.Get_size() + + + params = parseParamFile(configfilename, verbose=False) + + if threadNum == 0: + #print("--- DELIGHT-APPLY ---") + logger.info("--- DELIGHT-APPLY ---") + + + # Read filter coefficients, compute normalization of filters + bandCoefAmplitudes, bandCoefPositions, bandCoefWidths, norms = readBandCoefficients(params) + numBands = bandCoefAmplitudes.shape[0] + + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + f_mod_interp = readSEDs(params) + nt = f_mod_interp.shape[0] + nz = redshiftGrid.size + + dir_seds = params['templates_directory'] + dir_filters = params['bands_directory'] + lambdaRef = params['lambdaRef'] + sed_names = params['templates_names'] + f_mod_grid = np.zeros((redshiftGrid.size, len(sed_names),len(params['bandNames']))) + + + for t, sed_name in enumerate(sed_names): + f_mod_grid[:, t, :] = np.loadtxt(dir_seds + '/' + sed_name +'_fluxredshiftmod.txt') + + numZbins = redshiftDistGrid.size - 1 + numZ = redshiftGrid.size + + numObjectsTraining = np.sum(1 for line in open(params['training_catFile'])) + numObjectsTarget = np.sum(1 for line in open(params['target_catFile'])) + redshiftsInTarget = ('redshift' in params['target_bandOrder']) + Ncompress = params['Ncompress'] + + firstLine = int(threadNum * numObjectsTarget / float(numThreads)) + lastLine = int(min(numObjectsTarget,(threadNum + 1) * numObjectsTarget / float(numThreads))) + numLines = lastLine - firstLine + + if threadNum == 0: + msg= 'Number of Training Objects ' + str(numObjectsTraining) + logger.info(msg) + + msg='Number of Target Objects ' + str(numObjectsTarget) + logger.info(msg) + + comm.Barrier() + + msg= 'Thread '+ str(threadNum) + ' , analyzes lines ' + str(firstLine) + ' to ' + str( lastLine) + logger.info(msg) + + DL = approx_DL() + gp = PhotozGP(f_mod_interp, + bandCoefAmplitudes, bandCoefPositions, bandCoefWidths, + params['lines_pos'], params['lines_width'], + params['V_C'], params['V_L'], + params['alpha_C'], params['alpha_L'], + redshiftGridGP, use_interpolators=True) + + # Create local files to store results + numMetrics = 7 + len(params['confidenceLevels']) + localPDFs = np.zeros((numLines, numZ)) + localMetrics = np.zeros((numLines, numMetrics)) + localCompressIndices = np.zeros((numLines, Ncompress), dtype=int) + localCompEvidences = np.zeros((numLines, Ncompress)) + + # Looping over chunks of the training set to prepare model predictions over z + numChunks = params['training_numChunks'] + for chunk in range(numChunks): + TR_firstLine = int(chunk * numObjectsTraining / float(numChunks)) + TR_lastLine = int(min(numObjectsTraining, (chunk + 1) * numObjectsTarget / float(numChunks))) + targetIndices = np.arange(TR_firstLine, TR_lastLine) + numTObjCk = TR_lastLine - TR_firstLine + redshifts = np.zeros((numTObjCk, )) + model_mean = np.zeros((numZ, numTObjCk, numBands)) + model_covar = np.zeros((numZ, numTObjCk, numBands)) + bestTypes = np.zeros((numTObjCk, ), dtype=int) + ells = np.zeros((numTObjCk, ), dtype=int) + + # loop on training data and training GP coefficients produced by delight_learn + # It fills the model_mean and model_covar predicted by GP + loc = TR_firstLine - 1 + trainingDataIter = getDataFromFile(params, TR_firstLine, TR_lastLine,prefix="training_", ftype="gpparams") + + # loop on training data to load the GP parameter + for loc, (z, ell, bands, X, B, flatarray) in enumerate(trainingDataIter): + t1 = time() + redshifts[loc] = z # redshift of all training samples + gp.setCore(X, B, nt,flatarray[0:nt+B+B*(B+1)//2]) + bestTypes[loc] = gp.bestType # retrieve the best-type found by delight-learn + ells[loc] = ell # retrieve the luminosity parameter l + + # here is the model prediction of Gaussian Process for that particular trainning galaxy + model_mean[:, loc, :], model_covar[:, loc, :] = gp.predictAndInterpolate(redshiftGrid, ell=ell) + t2 = time() + # print(loc, t2-t1) + + #Redshift prior on training galaxy + # p_t = params['p_t'][bestTypes][None, :] + # p_z_t = params['p_z_t'][bestTypes][None, :] + # compute the prior for taht training sample + prior = np.exp(-0.5*((redshiftGrid[:, None]-redshifts[None, :]) /params['zPriorSigma'])**2) + # prior[prior < 1e-6] = 0 + # prior *= p_t * redshiftGrid[:, None] * + # np.exp(-0.5 * redshiftGrid[:, None]**2 / p_z_t) / p_z_t + + if params['useCompression'] and params['compressionFilesFound']: + fC = open(params['compressMargLikFile']) + fCI = open(params['compressIndicesFile']) + itCompM = itertools.islice(fC, firstLine, lastLine) + iterCompI = itertools.islice(fCI, firstLine, lastLine) + + targetDataIter = getDataFromFile(params, firstLine, lastLine,prefix="target_", getXY=False, CV=False) + + # loop on target samples + for loc, (z, normedRefFlux, bands, fluxes, fluxesVar, bCV, dCV, dVCV) in enumerate(targetDataIter): + t1 = time() + ell_hat_z = normedRefFlux * 4 * np.pi * params['fluxLuminosityNorm'] * (DL(redshiftGrid)**2. * (1+redshiftGrid)) + ell_hat_z[:] = 1 + if params['useCompression'] and params['compressionFilesFound']: + indices = np.array(next(iterCompI).split(' '), dtype=int) + sel = np.in1d(targetIndices, indices, assume_unique=True) + # same likelihood as for template fitting + like_grid2 = approx_flux_likelihood(fluxes,fluxesVar,model_mean[:, sel, :][:, :, bands], + f_mod_covar=model_covar[:, sel, :][:, :, bands], + marginalizeEll=True, normalized=False, + ell_hat=ell_hat_z, + ell_var=(ell_hat_z*params['ellPriorSigma'])**2) + like_grid *= prior[:, sel] + else: + like_grid = np.zeros((nz, model_mean.shape[1])) + # same likelihood as for template fitting, but cython + approx_flux_likelihood_cy( + like_grid, nz, model_mean.shape[1], bands.size, + fluxes, fluxesVar, # target galaxy fluxes and variance + model_mean[:, :, bands], # prediction with Gaussian process + model_covar[:, :, bands], + ell_hat=ell_hat_z, # it will find internally the ell + ell_var=(ell_hat_z*params['ellPriorSigma'])**2) + like_grid *= prior[:, :] #likelihood multiplied by redshift training galaxies priors + t2 = time() + localPDFs[loc, :] += like_grid.sum(axis=1) # the final redshift posterior is sum over training galaxies posteriors + + # compute the evidence for each model + evidences = np.trapz(like_grid, x=redshiftGrid, axis=0) + t3 = time() + + if params['useCompression'] and not params['compressionFilesFound']: + if localCompressIndices[loc, :].sum() == 0: + sortind = np.argsort(evidences)[::-1][0:Ncompress] + localCompressIndices[loc, :] = targetIndices[sortind] + localCompEvidences[loc, :] = evidences[sortind] + else: + dind = np.concatenate((targetIndices,localCompressIndices[loc, :])) + devi = np.concatenate((evidences,localCompEvidences[loc, :])) + sortind = np.argsort(devi)[::-1][0:Ncompress] + localCompressIndices[loc, :] = dind[sortind] + localCompEvidences[loc, :] = devi[sortind] + + if chunk == numChunks - 1\ + and redshiftsInTarget\ + and localPDFs[loc, :].sum() > 0: + localMetrics[loc, :] = computeMetrics(z, redshiftGrid,localPDFs[loc, :],params['confidenceLevels']) + t4 = time() + if loc % 100 == 0: + print(loc, t2-t1, t3-t2, t4-t3) + + if params['useCompression'] and params['compressionFilesFound']: + fC.close() + fCI.close() + + comm.Barrier() + + if threadNum == 0: + globalPDFs = np.zeros((numObjectsTarget, numZ)) + globalCompressIndices = np.zeros((numObjectsTarget, Ncompress), dtype=int) + globalCompEvidences = np.zeros((numObjectsTarget, Ncompress)) + globalMetrics = np.zeros((numObjectsTarget, numMetrics)) + else: + globalPDFs = None + globalCompressIndices = None + globalCompEvidences = None + globalMetrics = None + + firstLines = [int(k*numObjectsTarget/numThreads) for k in range(numThreads)] + lastLines = [int(min(numObjectsTarget, (k+1)*numObjectsTarget/numThreads)) for k in range(numThreads)] + numLines = [lastLines[k] - firstLines[k] for k in range(numThreads)] + + sendcounts = tuple([numLines[k] * numZ for k in range(numThreads)]) + displacements = tuple([firstLines[k] * numZ for k in range(numThreads)]) + comm.Gatherv(localPDFs,[globalPDFs, sendcounts, displacements, MPI.DOUBLE]) + + sendcounts = tuple([numLines[k] * Ncompress for k in range(numThreads)]) + displacements = tuple([firstLines[k] * Ncompress for k in range(numThreads)]) + comm.Gatherv(localCompressIndices,[globalCompressIndices, sendcounts, displacements, MPI.LONG]) + comm.Gatherv(localCompEvidences,[globalCompEvidences, sendcounts, displacements, MPI.DOUBLE]) + comm.Barrier() + + sendcounts = tuple([numLines[k] * numMetrics for k in range(numThreads)]) + displacements = tuple([firstLines[k] * numMetrics for k in range(numThreads)]) + comm.Gatherv(localMetrics,[globalMetrics, sendcounts, displacements, MPI.DOUBLE]) + + comm.Barrier() + + if threadNum == 0: + fmt = '%.2e' + fname = params['redshiftpdfFileComp'] if params['compressionFilesFound']\ + else params['redshiftpdfFile'] + np.savetxt(fname, globalPDFs, fmt=fmt) + if redshiftsInTarget: + np.savetxt(params['metricsFile'], globalMetrics, fmt=fmt) + if params['useCompression'] and not params['compressionFilesFound']: + np.savetxt(params['compressMargLikFile'],globalCompEvidences, fmt=fmt) + np.savetxt(params['compressIndicesFile'],globalCompressIndices, fmt="%i") + + +#----------------------------------------------------------------------------------------- +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start Delight Learn.py" + logger.info(msg) + logger.info("--- Process Delight Learn ---") + + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + delightApply(sys.argv[1]) diff --git a/rail/estimation/algos/include_delightPZ/delightLearn.py b/rail/estimation/algos/include_delightPZ/delightLearn.py new file mode 100644 index 00000000..53cd7697 --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/delightLearn.py @@ -0,0 +1,165 @@ +################################################################################################################################## +# +# script : delight-learn.py +# +# input : 'training_catFile' +# output : localData or reducedData usefull for Gaussian Process in 'training_paramFile' +# - find the normalisation of the flux and the best galaxy type +############################################################################################################################ +import sys +from mpi4py import MPI +import numpy as np +from delight.io import * +from delight.utils import * +from delight.photoz_gp import PhotozGP +from delight.photoz_kernels import Photoz_mean_function, Photoz_kernel + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + +def delightLearn(configfilename): + """ + + :param configfilename: + :return: + """ + + + comm = MPI.COMM_WORLD + threadNum = comm.Get_rank() + numThreads = comm.Get_size() + + #parse arguments + + params = parseParamFile(configfilename, verbose=False) + + if threadNum == 0: + logger.info("--- DELIGHT-LEARN ---") + + # Read filter coefficients, compute normalization of filters + bandCoefAmplitudes, bandCoefPositions, bandCoefWidths, norms = readBandCoefficients(params) + numBands = bandCoefAmplitudes.shape[0] + + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + + f_mod = readSEDs(params) + + numObjectsTraining = np.sum(1 for line in open(params['training_catFile'])) + + msg= 'Number of Training Objects ' + str(numObjectsTraining) + logger.info(msg) + + + firstLine = int(threadNum * numObjectsTraining / numThreads) + lastLine = int(min(numObjectsTraining,(threadNum + 1) * numObjectsTraining / numThreads)) + numLines = lastLine - firstLine + + comm.Barrier() + msg ='Thread ' + str(threadNum) + ' , analyzes lines ' + str(firstLine) + ' , to ' + str(lastLine) + logger.info(msg) + + DL = approx_DL() + gp = PhotozGP(f_mod, bandCoefAmplitudes, bandCoefPositions, bandCoefWidths, + params['lines_pos'], params['lines_width'], + params['V_C'], params['V_L'], + params['alpha_C'], params['alpha_L'], + redshiftGridGP, use_interpolators=True) + + B = numBands + numCol = 3 + B + B*(B+1)//2 + B + f_mod.shape[0] + localData = np.zeros((numLines, numCol)) + fmt = '%i ' + '%.12e ' * (localData.shape[1] - 1) + + loc = - 1 + crossValidate = params['training_crossValidate'] + trainingDataIter1 = getDataFromFile(params, firstLine, lastLine,prefix="training_", getXY=True,CV=crossValidate) + + + if crossValidate: + chi2sLocal = None + bandIndicesCV, bandNamesCV, bandColumnsCV,bandVarColumnsCV, redshiftColumnCV = readColumnPositions(params, prefix="training_CV_", refFlux=False) + + for z, normedRefFlux,\ + bands, fluxes, fluxesVar,\ + bandsCV, fluxesCV, fluxesVarCV,\ + X, Y, Yvar in trainingDataIter1: + + loc += 1 + + themod = np.zeros((1, f_mod.shape[0], bands.size)) + for it in range(f_mod.shape[0]): + for ib, band in enumerate(bands): + themod[0, it, ib] = f_mod[it, band](z) + + # really calibrate the luminosity parameter l compared to the model + # according the best type of galaxy + chi2_grid, ellMLs = scalefree_flux_likelihood(fluxes,fluxesVar,themod,returnChi2=True) + + bestType = np.argmin(chi2_grid) # best type + ell = ellMLs[0, bestType] # the luminosity factor + X[:, 2] = ell + + gp.setData(X, Y, Yvar, bestType) + lB = bands.size + localData[loc, 0] = lB + localData[loc, 1] = z + localData[loc, 2] = ell + localData[loc, 3:3+lB] = bands + localData[loc, 3+lB:3+f_mod.shape[0]+lB+lB*(lB+1)//2+lB] = gp.getCore() + + if crossValidate: + model_mean, model_covar = gp.predictAndInterpolate(np.array([z]), ell=ell) + if chi2sLocal is None: + chi2sLocal = np.zeros((numObjectsTraining, bandIndicesCV.size)) + + ind = np.array([list(bandIndicesCV).index(b) for b in bandsCV]) + + chi2sLocal[firstLine + loc, ind] = - 0.5 * (model_mean[0, bandsCV] - fluxesCV)**2 /(model_covar[0, bandsCV] + fluxesVarCV) + + + # use MPI to get the totals + comm.Barrier() + if threadNum == 0: + reducedData = np.zeros((numObjectsTraining, numCol)) + else: + reducedData = None + + if crossValidate: + chi2sGlobal = np.zeros_like(chi2sLocal) + comm.Allreduce(chi2sLocal, chi2sGlobal, op=MPI.SUM) + comm.Barrier() + + firstLines = [int(k*numObjectsTraining/numThreads) for k in range(numThreads)] + lastLines = [int(min(numObjectsTraining, (k+1)*numObjectsTraining/numThreads)) for k in range(numThreads)] + sendcounts = tuple([(lastLines[k] - firstLines[k]) * numCol for k in range(numThreads)]) + displacements = tuple([firstLines[k] * numCol for k in range(numThreads)]) + + comm.Gatherv(localData, [reducedData, sendcounts, displacements, MPI.DOUBLE]) + comm.Barrier() + + # parameters for the GP process on traniing data are transfered to reduced data and saved in file + #'training_paramFile' + if threadNum == 0: + np.savetxt(params['training_paramFile'], reducedData, fmt=fmt) + if crossValidate: + np.savetxt(params['training_CVfile'], chi2sGlobal) + + +#----------------------------------------------------------------------------------------- +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start Delight Learn.py" + logger.info(msg) + logger.info("--- Process Delight Learn ---") + + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + delightLearn(sys.argv[1]) diff --git a/rail/estimation/algos/include_delightPZ/getDelightRedshiftEstimation.py b/rail/estimation/algos/include_delightPZ/getDelightRedshiftEstimation.py new file mode 100644 index 00000000..90c8613a --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/getDelightRedshiftEstimation.py @@ -0,0 +1,73 @@ +import sys +import os +import numpy as np +from functools import reduce + +import pprint + +from delight.io import * +from delight.utils import * +import h5py + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + + + +def getDelightRedshiftEstimation(configfilename,chunknum,nsize,index_sel): + """ + zmode,widths = getDelightRedshiftEstimation(delightparamfilechunk,self.chunknum,nsize,indexes_sel) + + input args: + - nsize : size of arrays to return + - index_sel : indexes in final arays of processed redshits by delight + + + :return: + """ + + msg = "--- getDelightRedshiftEstimation({}) for chunk {}---".format(nsize,chunknum) + logger.info(msg) + + # initialize arrays to be returned + zmode = np.full(nsize, fill_value=-1,dtype=np.float) + widths = np.full(nsize, fill_value=-1, dtype=np.float) + + params = parseParamFile(configfilename, verbose=False) + + # redshiftGrid has nz size + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + + # the pdfs have (m x nz) size + # where m is the number of redshifts calculated by delight + # nz is the number of redshifts + pdfs = np.loadtxt(params['redshiftpdfFile']) + pdfs /= np.trapz(pdfs, x=redshiftGrid, axis=1)[:, None] + + + # find the index of the redshift where there is the mode + # the following arrays have size m + indexes_of_zmode = np.argmax(pdfs,axis=1) + + redshifts_of_zmode = redshiftGrid[indexes_of_zmode] + + + # array of zshift (z-zmode) : of size (m x nz) + zshifts_of_mode = redshiftGrid[np.newaxis,:]-redshifts_of_zmode[:,np.newaxis] + + # array of variances + variance_of_z = np.average(zshifts_of_mode**2,axis=1,weights=pdfs) + + + # copy only the processed redshifts and widths into the final arrays of size nsize + # for RAIL + zmode[index_sel] = redshifts_of_zmode + widths[index_sel] = np.sqrt(variance_of_z)*(1.0 + redshifts_of_zmode) + + + return zmode,widths + diff --git a/rail/estimation/algos/include_delightPZ/libPriorPZ.py b/rail/estimation/algos/include_delightPZ/libPriorPZ.py new file mode 100644 index 00000000..a7a33908 --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/libPriorPZ.py @@ -0,0 +1,159 @@ +####################################################################################### +# +# script : libpriorPZ +# +# Provide a library of priors on photoZ +# +# author : Sylvie Dagoret-Campagne +# affiliation : IJCLab/IN2P3/CNRS +# +# from https://github.com/ixkael/Photoz-tools +# +###################################################################################### +import sys +import numpy as np +from scipy.interpolate import interp1d +from pprint import pprint + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + + +def mknames(nt): + return ['Elliptical ' + str(i + 1) for i in range(nt[0])] \ + + ['Spiral ' + str(i + 1) for i in range(nt[1])] \ + + ['Starburst ' + str(i + 1) for i in range(nt[2])] + + + +# This is the prior HDFN prior from Benitez 2000, adapted from the BPZ code. +# This could be replaced with any redshift, magnitude, and type distribution. +def bpz_prior(z, m, nt): + """ + bpz_prior(z, m, nt): + + - z grid of redshift + - m maximum magnitude + - nt : number of types + + """ + nz = len(z) + momin_hdf = 20. + if m > 32.: m = 32. + if m < 20.: m = 20. + # nt Templates = nell Elliptical + nsp Spiral + nSB starburst + try: # nt is a list of 3 values + nell, nsp, nsb = nt + except: # nt is a single value + nell = 1 # 1 Elliptical in default template set + nsp = 2 # 2 Spirals in default template set + nsb = nt - nell - nsp # rest Irr/SB + nn = nell, nsp, nsb + nt = sum(nn) + # See Table 1 of Benitez00 + a = 2.465, 1.806, 0.906 + zo = 0.431, 0.390, 0.0626 + km = 0.0913, 0.0636, 0.123 + k_t = 0.450, 0.147 + a = np.repeat(a, nn) + zo = np.repeat(zo, nn) + km = np.repeat(km, nn) + k_t = np.repeat(k_t, nn[:2]) + + # Fractions expected at m = 20: 35% E/S0, 50% Spiral, 15% Irr + fo_t = 0.35, 0.5 + fo_t = fo_t / np.array(nn[:2]) + fo_t = np.repeat(fo_t, nn[:2]) + + dm = m - momin_hdf + zmt = np.clip(zo + km * dm, 0.01, 15.) + zmt_at_a = zmt ** (a) + zt_at_a = np.power.outer(z, a) + + # Morphological fractions + nellsp = nell + nsp + f_t = np.zeros((len(a),), float) + f_t[:nellsp] = fo_t * np.exp(-k_t * dm) + f_t[nellsp:] = (1. - np.add.reduce(f_t[:nellsp])) / float(nsb) + + # Formula: zm=zo+km*(m_m_min) and p(z|T,m)=(z**a)*exp(-(z/zm)**a) + p_i = zt_at_a[:nz, :nt] * np.exp(-np.clip(zt_at_a[:nz, :nt] / zmt_at_a[:nt], 0., 700.)) + + # This eliminates the very low level tails of the priors + norm = np.add.reduce(p_i[:nz, :nt], 0) + p_i[:nz, :nt] = np.where(np.less(p_i[:nz, :nt] / norm[:nt], 1e-2 / float(nz)), + 0., p_i[:nz, :nt] / norm[:nt]) + norm = np.add.reduce(p_i[:nz, :nt], 0) + p_i[:nz, :nt] = p_i[:nz, :nt] / norm[:nt] * f_t[:nt] + return p_i # return 2D template nz x nt + + +def libPriorPZ(z_grid,maglim,nt=8): + """ + + :return: + """ + + msg = "--- libPriorPZ" + #logger.info(msg) + + # Just some boolean indexing of templates used. Needed later for some BPZ fcts. + selectedtemplates = np.repeat(False, nt) + + # Using all templates + templatetypesnb = (1, 2, 5) # nb of ellipticals, spirals, and starburst used in the 8-template library. + selectedtemplates[:] = True + + # Uncomment that to use three templates using + # templatetypesnb = (1,1,1) #(1,2,8-3) + # selectedtemplates[0:1] = True + nt = sum(templatetypesnb) + + ellipticals = ['El_B2004a.sed'][0:templatetypesnb[0]] + spirals = ['Sbc_B2004a.sed', 'Scd_B2004a.sed'][0:templatetypesnb[1]] + irregulars = ['Im_B2004a.sed', 'SB3_B2004a.sed', 'SB2_B2004a.sed', + 'ssp_25Myr_z008.sed', 'ssp_5Myr_z008.sed'][0:templatetypesnb[2]] + template_names = [nm.replace('.sed', '') for nm in ellipticals + spirals + irregulars] + + # Use the p(z,t,m) distribution defined above + m = maglim # some reference magnitude + p_z__t_m = bpz_prior(z_grid, m, templatetypesnb) # 2D template nz x nt + + # Convenient function for template names + def mknames(nt): + return ['Elliptical ' + str(i + 1) for i in range(nt[0])] \ + + ['Spiral ' + str(i + 1) for i in range(nt[1])] \ + + ['Starburst ' + str(i + 1) for i in range(nt[2])] + + names = mknames(templatetypesnb) + + return p_z__t_m # return 2D template nz x nt + + + + + +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start libpriorPZ.py" + logger.info(msg) + logger.info("--- libPriorPZ ---") + + z_grid_binsize = 0.001 + z_grid_edges = np.arange(0.0, 3.0, z_grid_binsize) + z_grid = (z_grid_edges[1:] + z_grid_edges[:-1]) / 2. + + m = 26.0 # some reference magnitude + nt=8 + + p_z__t_m = libPriorPZ(z_grid,maglim=m,nt=nt) + + np.set_printoptions(threshold=20, edgeitems=10, linewidth=140, + formatter=dict(float=lambda x: "%.3e" % x)) # float arrays %.3g + print(p_z__t_m ) \ No newline at end of file diff --git a/rail/estimation/algos/include_delightPZ/makeConfigParam.py b/rail/estimation/algos/include_delightPZ/makeConfigParam.py new file mode 100644 index 00000000..a9d6b402 --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/makeConfigParam.py @@ -0,0 +1,699 @@ +#################################################################################################### +# Script name : makeConfigParam.py +# +# Generate Config parameter required by Delight +# +# Some parameters are read from the from the rail configuration file +# Some other parameter are hardcoded in this file +# The fina goal is to retrieve those parameters from RAIL config file +##################################################################################################### +import sys +import numpy as np +from scipy.interpolate import interp1d +from scipy.optimize import leastsq + +from delight.utils import * +from delight.io import * + +import coloredlogs +import logging + + +import os +import yaml + +from pkg_resources import resource_filename + +# Create a logger object. +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s %(name)s[%(process)d] %(levelname)s %(message)s') + + +def makeConfigParam(path,inputs_rail): + """ + makeConfigParam(path,inputs_rail) + + generate Configuration parameter file in ascii. This file is decoded by Delight functions with argparse + + : inputs: + - path : where the FILTERS and SEDs datafiles used by Delight initialisation are stored, + - inputs_rail : RAIL parameter files + + Either the parameters used by Delight are hardcoded here of the can be setup by RAIL config strcture (yaml) in inputs_rail + + :return: paramfile_txt , the string for the configuration file. RAIL will write itself this file. + """ + + logger.debug("__name__:"+__name__) + logger.debug("__file__"+__file__) + + msg = "----- makeConfigParam ------" + logger.info(msg) + + logger.debug(" received path = "+ path) + #logger.debug(" received input_rail = " + inputs_rail) + + # 1) Let 's create a parameter file from scratch. + + #paramfile_txt = "\n" + #paramfile_txt += \ + paramfile_txt = \ +""" +# DELIGHT parameter file +# Syntactic rules: +# - You can set parameters with : or = +# - Lines starting with # or ; will be ignored +# - Multiple values (band names, band orders, confidence levels) +# must beb separated by spaces +# - The input files should contain numbers separated with spaces. +# - underscores mean unused column +""" + + # 2) Filter Section + paramfile_txt += "\n" + paramfile_txt += \ +""" +[Bands] +names: lsst_u lsst_g lsst_r lsst_i lsst_z lsst_y +""" + + paramfile_txt += "directory: " + os.path.join(path, 'FILTERS') + + paramfile_txt += \ +""" +numCoefs: 15 +bands_verbose: True +bands_debug: True +bands_makeplots: False +""" + + # 3) Template Section + + paramfile_txt += \ +""" +[Templates] +""" + paramfile_txt += "directory: " + os.path.join(path, 'CWW_SEDs') + + paramfile_txt += \ +""" +names: El_B2004a Sbc_B2004a Scd_B2004a SB3_B2004a SB2_B2004a Im_B2004a ssp_25Myr_z008 ssp_5Myr_z008 +p_t: 0.27 0.26 0.25 0.069 0.021 0.11 0.0061 0.0079 +p_z_t:0.23 0.39 0.33 0.31 1.1 0.34 1.2 0.14 +lambdaRef: 4.5e3 +""" + + # 4) Simulation Section + + paramfile_txt += \ +""" +[Simulation] +numObjects: 1000 +noiseLevel: 0.03 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +trainingFile: data_lsst/galaxies-fluxredshifts.txt +targetFile: data_lsst/galaxies-fluxredshifts2.txt +""" + else: + thepath=inputs_rail["tempdatadir"] + paramfile_txt += "trainingFile: " + os.path.join(thepath, 'galaxies-fluxredshifts.txt') + paramfile_txt += "\n" + paramfile_txt += "targetFile: " + os.path.join(thepath, 'galaxies-fluxredshifts2.txt') + paramfile_txt += "\n" + + + + # 5) Training Section + + paramfile_txt += \ +""" +[Training] +""" + if inputs_rail == None: + paramfile_txt += \ +""" +catFile: data_lsst/galaxies-fluxredshifts.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "catFile: " + os.path.join(thepath, 'galaxies-fluxredshifts.txt') + + paramfile_txt += \ +""" +bandOrder: lsst_u lsst_u_var lsst_g lsst_g_var lsst_r lsst_r_var lsst_i lsst_i_var lsst_z lsst_z_var lsst_y lsst_y_var redshift +referenceBand: lsst_i +extraFracFluxError: 1e-4 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +paramFile: data_lsst/galaxies-gpparams.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "paramFile: " + os.path.join(thepath, 'galaxies-gpparams.txt') + + + paramfile_txt += \ +""" +crossValidate: False +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +CVfile: data_lsst/galaxies-gpCV.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "CVfile: " + os.path.join(thepath, 'galaxies-gpCV.txt') + + paramfile_txt += \ +""" +crossValidationBandOrder: _ _ _ _ lsst_r lsst_r_var _ _ _ _ _ _ +numChunks: 1 +""" + + # 6) Estimation Section + + + paramfile_txt += \ +""" +[Target] +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +catFile: data_lsst/galaxies-fluxredshifts2.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "catFile: " + os.path.join(thepath, 'galaxies-fluxredshifts2.txt') + + paramfile_txt += \ +""" +bandOrder: lsst_u lsst_u_var lsst_g lsst_g_var lsst_r lsst_r_var lsst_i lsst_i_var lsst_z lsst_z_var lsst_y lsst_y_var redshift +referenceBand: lsst_r +extraFracFluxError: 1e-4 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +redshiftpdfFile: data_lsst/galaxies-redshiftpdfs.txt +redshiftpdfFileTemp: data_lsst/galaxies-redshiftpdfs-cww.txt +metricsFile: data_lsst/galaxies-redshiftmetrics.txt +metricsFileTemp: data_lsst/galaxies-redshiftmetrics-cww.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "redshiftpdfFile: " + os.path.join(thepath, 'galaxies-redshiftpdfs.txt') + paramfile_txt += "\n" + paramfile_txt += "redshiftpdfFileTemp: " + os.path.join(thepath, 'galaxies-redshiftpdfs-cww.txt') + paramfile_txt += "\n" + paramfile_txt += "metricsFile: " + os.path.join(thepath, 'galaxies-redshiftmetrics.txt') + paramfile_txt += "\n" + paramfile_txt += "metricsFileTemp: " + os.path.join(thepath, 'galaxies-redshiftmetrics-cww.txt') + + + paramfile_txt += \ +""" +useCompression: False +Ncompress: 10 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +compressIndicesFile: data_lsst/galaxies-compressionIndices.txt +compressMargLikFile: data_lsst/galaxies-compressionMargLikes.txt +redshiftpdfFileComp: data_lsst/galaxies-redshiftpdfs-comp.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "compressIndicesFile: " + os.path.join(thepath, 'galaxies-compressionIndices.txt') + paramfile_txt += "\n" + paramfile_txt += "compressMargLikFile: " + os.path.join(thepath, 'galaxies-compressionMargLikes.txt') + paramfile_txt += "\n" + paramfile_txt += "redshiftpdfFileComp: " + os.path.join(thepath, 'galaxies-redshiftpdfs-comp.txt') + paramfile_txt += "\n" + + # 7) Other Section + + if inputs_rail == None: + paramfile_txt += \ +""" +[Other] +rootDir: ./ +zPriorSigma: 0.2 +ellPriorSigma: 0.5 +fluxLuminosityNorm: 1.0 +alpha_C: 1.0e3 +V_C: 0.1 +alpha_L: 1.0e2 +V_L: 0.1 +lines_pos: 6500 5002.26 3732.22 +lines_width: 20.0 20.0 20.0 +""" + else: + zPriorSigma = inputs_rail["zPriorSigma"] + ellPriorSigma = inputs_rail["ellPriorSigma"] + fluxLuminosityNorm = inputs_rail["fluxLuminosityNorm"] + alpha_C = inputs_rail["alpha_C"] + V_C = inputs_rail["V_C"] + alpha_L = inputs_rail["alpha_L"] + V_L = inputs_rail["V_L"] + lineWidthSigma = inputs_rail["lineWidthSigma"] + + paramfile_txt += \ +""" +[Other] +rootDir: ./ +""" + + paramfile_txt += "zPriorSigma: " + str(zPriorSigma) + paramfile_txt += "\n" + paramfile_txt += "ellPriorSigma: " + str(ellPriorSigma) + paramfile_txt += "\n" + paramfile_txt += "fluxLuminosityNorm: " + str(fluxLuminosityNorm) + paramfile_txt += "\n" + paramfile_txt += "alpha_C: " + str(alpha_C) + paramfile_txt += "\n" + paramfile_txt += "V_C: " + str(V_C) + paramfile_txt += "\n" + paramfile_txt += "alpha_L: " + str(alpha_L) + paramfile_txt += "\n" + paramfile_txt += "V_L: " + str(V_L) + paramfile_txt += "\n" + paramfile_txt += "lines_pos: 6500 5002.26 3732.22 \n" + paramfile_txt += "\n" + paramfile_txt += "lines_width: " + str(lineWidthSigma) + " " + \ + str(lineWidthSigma) + " " + \ + str(lineWidthSigma) + " " + \ + str(lineWidthSigma) + " " + "\n" + + + if inputs_rail == None: + paramfile_txt += \ +""" +redshiftMin: 0.1 +redshiftMax: 1.101 +redshiftNumBinsGPpred: 100 +redshiftBinSize: 0.001 +redshiftDisBinSize: 0.2 +""" + else: + + msg = "Decode redshift parameter from RAIL config file" + logger.debug(msg) + + dlght_redshiftMin = inputs_rail["dlght_redshiftMin"] + dlght_redshiftMax = inputs_rail["dlght_redshiftMax"] + dlght_redshiftNumBinsGPpred = inputs_rail["dlght_redshiftNumBinsGPpred"] + dlght_redshiftBinSize = inputs_rail["dlght_redshiftBinSize"] + dlght_redshiftDisBinSize = inputs_rail["dlght_redshiftDisBinSize"] + + # will check later what to do with these parameters + + paramfile_txt += "redshiftMin: " + str(dlght_redshiftMin) + paramfile_txt += "\n" + paramfile_txt += "redshiftMax: " + str(dlght_redshiftMax) + paramfile_txt += "\n" + paramfile_txt += "redshiftNumBinsGPpred: " + str(dlght_redshiftNumBinsGPpred) + paramfile_txt += "\n" + paramfile_txt += "redshiftBinSize: " + str(dlght_redshiftBinSize) + paramfile_txt += "\n" + paramfile_txt += "redshiftDisBinSize: " + str(dlght_redshiftDisBinSize) + paramfile_txt += "\n" + + + + + paramfile_txt += \ +""" +confidenceLevels: 0.1 0.50 0.68 0.95 +""" + + + return paramfile_txt + + + + + + +def makeConfigParamChunk(path,inputs_rail,chunknum): + """ + makeConfigParam(path,inputs_rail) + + generate Configuration parameter file in ascii. This file is decoded by Delight functions with argparse + + : inputs: + - path : where the FILTERS and SEDs datafiles used by Delight initialisation are stored, + - inputs_rail : RAIL parameter files + - chunknum : chunknumber + + Either the parameters used by Delight are hardcoded here of the can be setup by RAIL config strcture (yaml) in inputs_rail + + :return: paramfile_txt , the string for the configuration file. RAIL will write itself this file. + """ + + logger.debug("__name__:"+__name__) + logger.debug("__file__"+__file__) + + msg = "----- makeConfigParamChunk ------" + logger.info(msg) + + logger.debug(" received path = "+ path) + #logger.debug(" received input_rail = " + inputs_rail) + + # 1) Let 's create a parameter file from scratch. + + #paramfile_txt = "\n" + #paramfile_txt += \ + paramfile_txt = \ +""" +# DELIGHT parameter file +# Syntactic rules: +# - You can set parameters with : or = +# - Lines starting with # or ; will be ignored +# - Multiple values (band names, band orders, confidence levels) +# must beb separated by spaces +# - The input files should contain numbers separated with spaces. +# - underscores mean unused column +""" + + # 2) Filter Section + paramfile_txt += "\n" + paramfile_txt += \ +""" +[Bands] +names: lsst_u lsst_g lsst_r lsst_i lsst_z lsst_y +""" + + paramfile_txt += "directory: " + os.path.join(path, 'FILTERS') + + paramfile_txt += \ +""" +numCoefs: 15 +bands_verbose: True +bands_debug: True +bands_makeplots: False +""" + + # 3) Template Section + + paramfile_txt += \ +""" +[Templates] +""" + paramfile_txt += "directory: " + os.path.join(path, 'CWW_SEDs') + + paramfile_txt += \ +""" +names: El_B2004a Sbc_B2004a Scd_B2004a SB3_B2004a SB2_B2004a Im_B2004a ssp_25Myr_z008 ssp_5Myr_z008 +p_t: 0.27 0.26 0.25 0.069 0.021 0.11 0.0061 0.0079 +p_z_t:0.23 0.39 0.33 0.31 1.1 0.34 1.2 0.14 +lambdaRef: 4.5e3 +""" + + # 4) Simulation Section + + paramfile_txt += \ +""" +[Simulation] +numObjects: 1000 +noiseLevel: 0.03 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +trainingFile: data_lsst/galaxies-fluxredshifts.txt +targetFile: data_lsst/galaxies-fluxredshifts2.txt +""" + else: + thepath=inputs_rail["tempdatadir"] + paramfile_txt += "trainingFile: " + os.path.join(thepath, 'galaxies-fluxredshifts.txt') + paramfile_txt += "\n" + #### + thepath = inputs_rail["tempdatadir"] + filename = 'galaxies-fluxredshifts2_{}.txt'.format(chunknum) + paramfile_txt += "targetFile: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + + # 5) Training Section + + paramfile_txt += \ +""" +[Training] +""" + if inputs_rail == None: + paramfile_txt += \ +""" +catFile: data_lsst/galaxies-fluxredshifts.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "catFile: " + os.path.join(thepath, 'galaxies-fluxredshifts.txt') + + paramfile_txt += \ +""" +bandOrder: lsst_u lsst_u_var lsst_g lsst_g_var lsst_r lsst_r_var lsst_i lsst_i_var lsst_z lsst_z_var lsst_y lsst_y_var redshift +referenceBand: lsst_i +extraFracFluxError: 1e-4 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +paramFile: data_lsst/galaxies-gpparams.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "paramFile: " + os.path.join(thepath, 'galaxies-gpparams.txt') + + + paramfile_txt += \ +""" +crossValidate: False +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +CVfile: data_lsst/galaxies-gpCV.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + paramfile_txt += "CVfile: " + os.path.join(thepath, 'galaxies-gpCV.txt') + + paramfile_txt += \ +""" +crossValidationBandOrder: _ _ _ _ lsst_r lsst_r_var _ _ _ _ _ _ +numChunks: 1 +""" + + # 6) Estimation Section + + + paramfile_txt += \ +""" +[Target] +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +catFile: data_lsst/galaxies-fluxredshifts2.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + filename = 'galaxies-fluxredshifts2_{}.txt'.format(chunknum) + paramfile_txt += "catFile: " + os.path.join(thepath, filename) + + paramfile_txt += \ +""" +bandOrder: lsst_u lsst_u_var lsst_g lsst_g_var lsst_r lsst_r_var lsst_i lsst_i_var lsst_z lsst_z_var lsst_y lsst_y_var redshift +referenceBand: lsst_r +extraFracFluxError: 1e-4 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +redshiftpdfFile: data_lsst/galaxies-redshiftpdfs.txt +redshiftpdfFileTemp: data_lsst/galaxies-redshiftpdfs-cww.txt +metricsFile: data_lsst/galaxies-redshiftmetrics.txt +metricsFileTemp: data_lsst/galaxies-redshiftmetrics-cww.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + filename = 'galaxies-redshiftpdfs_{}.txt'.format(chunknum) + paramfile_txt += "redshiftpdfFile: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + filename = 'galaxies-redshiftpdfs-cww_{}.txt'.format(chunknum) + paramfile_txt += "redshiftpdfFileTemp: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + filename = 'galaxies-redshiftmetrics_{}.txt'.format(chunknum) + paramfile_txt += "metricsFile: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + filename = 'galaxies-redshiftmetrics-cww_{}.txt'.format(chunknum) + paramfile_txt += "metricsFileTemp: " + os.path.join(thepath, filename) + + + paramfile_txt += \ +""" +useCompression: False +Ncompress: 10 +""" + + if inputs_rail == None: + paramfile_txt += \ +""" +compressIndicesFile: data_lsst/galaxies-compressionIndices.txt +compressMargLikFile: data_lsst/galaxies-compressionMargLikes.txt +redshiftpdfFileComp: data_lsst/galaxies-redshiftpdfs-comp.txt +""" + else: + thepath = inputs_rail["tempdatadir"] + filename = 'galaxies-compressionIndices_{}.txt'.format(chunknum) + paramfile_txt += "compressIndicesFile: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + filename = 'galaxies-compressionMargLikes_{}.txt'.format(chunknum) + paramfile_txt += "compressMargLikFile: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + filename = 'galaxies-redshiftpdfs-comp_{}.txt'.format(chunknum) + paramfile_txt += "redshiftpdfFileComp: " + os.path.join(thepath, filename) + paramfile_txt += "\n" + + # 7) Other Section + + if inputs_rail == None: + paramfile_txt += \ +""" +[Other] +rootDir: ./ +zPriorSigma: 0.2 +ellPriorSigma: 0.5 +fluxLuminosityNorm: 1.0 +alpha_C: 1.0e3 +V_C: 0.1 +alpha_L: 1.0e2 +V_L: 0.1 +lines_pos: 6500 5002.26 3732.22 +lines_width: 20.0 20.0 20.0 +""" + else: + zPriorSigma = inputs_rail["zPriorSigma"] + ellPriorSigma = inputs_rail["ellPriorSigma"] + fluxLuminosityNorm = inputs_rail["fluxLuminosityNorm"] + alpha_C = inputs_rail["alpha_C"] + V_C = inputs_rail["V_C"] + alpha_L = inputs_rail["alpha_L"] + V_L = inputs_rail["V_L"] + lineWidthSigma = inputs_rail["lineWidthSigma"] + + paramfile_txt += \ +""" +[Other] +rootDir: ./ +""" + + paramfile_txt += "zPriorSigma: "+ str(zPriorSigma) + paramfile_txt += "\n" + paramfile_txt += "ellPriorSigma: " + str(ellPriorSigma) + paramfile_txt += "\n" + paramfile_txt += "fluxLuminosityNorm: " + str(fluxLuminosityNorm) + paramfile_txt += "\n" + paramfile_txt += "alpha_C: " + str(alpha_C) + paramfile_txt += "\n" + paramfile_txt += "V_C: " + str(V_C) + paramfile_txt += "\n" + paramfile_txt += "alpha_L: " + str(alpha_L) + paramfile_txt += "\n" + paramfile_txt += "V_L: " + str(V_L) + paramfile_txt += "\n" + paramfile_txt += "lines_pos: 6500 5002.26 3732.22 \n" + paramfile_txt += "\n" + paramfile_txt += "lines_width: " + str(lineWidthSigma) + " " +\ + str(lineWidthSigma) + " " +\ + str(lineWidthSigma) + " " + \ + str(lineWidthSigma) + " " + "\n" + + + + + + + if inputs_rail == None: + paramfile_txt += \ +""" +redshiftMin: 0.1 +redshiftMax: 1.101 +redshiftNumBinsGPpred: 100 +redshiftBinSize: 0.001 +redshiftDisBinSize: 0.2 +""" + else: + + msg = "Decode redshift parameter from RAI config file" + logger.debug(msg) + + dlght_redshiftMin = inputs_rail["dlght_redshiftMin"] + dlght_redshiftMax = inputs_rail["dlght_redshiftMax"] + dlght_redshiftNumBinsGPpred = inputs_rail["dlght_redshiftNumBinsGPpred"] + dlght_redshiftBinSize = inputs_rail["dlght_redshiftBinSize"] + dlght_redshiftDisBinSize = inputs_rail["dlght_redshiftDisBinSize"] + + # will check later what to do with these parameters + + paramfile_txt += "redshiftMin: " + str(dlght_redshiftMin) + paramfile_txt += "\n" + paramfile_txt += "redshiftMax: " + str(dlght_redshiftMax) + paramfile_txt += "\n" + paramfile_txt += "redshiftNumBinsGPpred: " + str(dlght_redshiftNumBinsGPpred) + paramfile_txt += "\n" + paramfile_txt += "redshiftBinSize: " + str(dlght_redshiftBinSize) + paramfile_txt += "\n" + paramfile_txt += "redshiftDisBinSize: " + str(dlght_redshiftDisBinSize) + paramfile_txt += "\n" + + + + + paramfile_txt += \ +""" +confidenceLevels: 0.1 0.50 0.68 0.95 +""" + + + return paramfile_txt + +#----------------------------------------------------------------------------------------- +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start makeConfigParam." + logger.info(msg) + logger.info("--- Make configuration parameter ---") + + logger.debug("__name__:"+__name__) + logger.debug("__file__:"+__file__) + + datapath=resource_filename('delight', '../data') + + logger.debug("datapath = " + datapath) + + + + param_txt=makeConfigParam(datapath,None) + + logger.info(param_txt) diff --git a/rail/estimation/algos/include_delightPZ/processFilters.py b/rail/estimation/algos/include_delightPZ/processFilters.py new file mode 100644 index 00000000..fdf5c2f0 --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/processFilters.py @@ -0,0 +1,165 @@ +#################################################################################################### +# Script name : processFilters.py +# +# fit the band filters with a gaussian mixture +# if make_plot, save images +# +# output file : band + '_gaussian_coefficients.txt' +##################################################################################################### +import sys +import numpy as np +from scipy.interpolate import interp1d +from scipy.optimize import leastsq + +from delight.utils import * +from delight.io import * + +import coloredlogs +import logging + +# Create a logger object. +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s %(name)s[%(process)d] %(levelname)s %(message)s') + + +def processFilters(configfilename): + """ + processFilters(configfilename) + + Develop filter transmission functions as a Gaussian Kernel regression + + : input file : the configuration file + :return: + """ + + msg="----- processFilters ------" + logger.info(msg) + + + params = parseParamFile(configfilename, verbose=False, catFilesNeeded=False) + + numCoefs = params["numCoefs"] + bandNames = params['bandNames'] + make_plots= params['bands_makeplots'] + + fmt = '.res' + max_redshift = params['redshiftMax'] # for plotting purposes + root = params['bands_directory'] + + if make_plots: + import matplotlib.pyplot as plt + cm = plt.get_cmap('brg') + num = len(bandNames) + cols = [cm(i/num) for i in range(num)] + + + # Function we will optimize + # Gaussian function representing filter + def dfunc(p, x, yd): + y = 0*x + n = p.size//2 + for i in range(n): + y += np.abs(p[i]) * np.exp(-0.5*((mus[i]-x)/np.abs(p[n+i]))**2.0) + return yd - y + + if make_plots: + fig0, ax0 = plt.subplots(1, 1, figsize=(8.2, 4)) + + # Loop over bands + for iband, band in enumerate(bandNames): + + fname_in = root + '/' + band + fmt + data = np.genfromtxt(fname_in) + coefs = np.zeros((numCoefs, 3)) + # wavelength - transmission function + x, y = data[:, 0], data[:, 1] + #y /= x # divide by lambda + # Only consider range where >1% max + ind = np.where(y > 0.01*np.max(y))[0] + lambdaMin, lambdaMax = x[ind[0]], x[ind[-1]] + + # Initialize values for amplitude and width of the components + sig0 = np.repeat((lambdaMax-lambdaMin)/numCoefs/4, numCoefs) + # Components uniformly distributed in the range + mus = np.linspace(lambdaMin+sig0[0], lambdaMax-sig0[-1], num=numCoefs) + amp0 = interp1d(x, y)(mus) + p0 = np.concatenate((amp0, sig0)) + print(band, end=" ") + + # fit + popt, pcov = leastsq(dfunc, p0, args=(x, y)) + coefs[:, 0] = np.abs(popt[0:numCoefs]) # amplitudes + coefs[:, 1] = mus # positions + coefs[:, 2] = np.abs(popt[numCoefs:2*numCoefs]) # widths + + # output for gaussian regression fit coefficients + fname_out = root + '/' + band + '_gaussian_coefficients.txt' + np.savetxt(fname_out, coefs, header=fname_in) + + xf = np.linspace(lambdaMin, lambdaMax, num=1000) + yy = 0*xf + for i in range(numCoefs): + yy += coefs[i, 0] * np.exp(-0.5*((coefs[i, 1] - xf)/coefs[i, 2])**2.0) + + if make_plots: + fig, ax = plt.subplots(figsize=(8, 4)) + ax.plot(x[ind], y[ind], lw=3, label='True filter', c='k') + ax.plot(xf, yy, lw=2, c='r', label='Gaussian fit') + # ax0.plot(x[ind], y[ind], lw=3, label=band, color=cols[iband]) + ax0.plot(xf, yy, lw=3, label=band, color=cols[iband]) + + coefs_redshifted = 1*coefs + coefs_redshifted[:, 1] /= (1. + max_redshift) + coefs_redshifted[:, 2] /= (1. + max_redshift) + lambdaMin_redshifted, lambdaMax_redshifted\ + = lambdaMin / (1. + max_redshift), lambdaMax / (1. + max_redshift) + xf = np.linspace(lambdaMin_redshifted, lambdaMax_redshifted, num=1000) + yy = 0*xf + for i in range(numCoefs): + yy += coefs_redshifted[i, 0] *\ + np.exp(-0.5*((coefs_redshifted[i, 1] - xf) / + coefs_redshifted[i, 2])**2.0) + + if make_plots: + ax.plot(xf, yy, lw=2, c='b', label='G fit at z='+str(max_redshift)) + title = band + ' band (' + fname_in +\ + ') with %i' % numCoefs+' components' + ax.set_title(title) + ax.set_ylim([0, data[:, 1].max()*1.2]) + ax.set_yticks([]) + ax.set_xlabel('$\lambda$') + ax.legend(loc='upper center', frameon=False, ncol=3) + + fig.tight_layout() + fname_fig = root + '/' + band + '_gaussian_approximation.png' + fig.savefig(fname_fig) + + if make_plots: + ax0.legend(loc='upper center', frameon=False, ncol=4) + ylims = ax0.get_ylim() + ax0.set_ylim([0, 1.4*ylims[1]]) + ax0.set_yticks([]) + ax0.set_xlabel(r'$\lambda$') + fig0.tight_layout() + fname_fig = root + '/allbands.pdf' + fig0.savefig(fname_fig) + + + +#----------------------------------------------------------------------------------------- +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start processFilters.py" + logger.info(msg) + logger.info("--- Process FILTERS ---") + + #numCoefs = 7 # number of components for the fit + #numCoefs = 21 # for lsst the transmission is too wavy ,number of components for the fit + #make_plots = True + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + processFilters(sys.argv[1]) diff --git a/rail/estimation/algos/include_delightPZ/processSEDs.py b/rail/estimation/algos/include_delightPZ/processSEDs.py new file mode 100644 index 00000000..bbebc3cb --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/processSEDs.py @@ -0,0 +1,113 @@ +#################################################################################################### +# +# script : processSED.py +# +# process the library of SEDs and project them onto the filters, (for the mean fct of the GP) +# (which may take a few minutes depending on the settings you set): +# +# output file : sed_name + '_fluxredshiftmod.txt' +###################################################################################################### + +import sys +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d + + +from delight.io import * +from delight.utils import * + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + + + +def processSEDs(configfilename): + """ + + processSEDs(configfilename) + + Compute the The Flux expected in each band for redshifts in the grid + : input file : the configuration file + + :return: produce the file of flux-redshift in bands + """ + + + + logger.info("--- Process SED ---") + + # decode the parameters + params = parseParamFile(configfilename, verbose=False, catFilesNeeded=False) + bandNames = params['bandNames'] + dir_seds = params['templates_directory'] + dir_filters = params['bands_directory'] + lambdaRef = params['lambdaRef'] + sed_names = params['templates_names'] + fmt = '.dat' + + # Luminosity Distnace + DL = approx_DL() + + #redshift grid + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + numZ = redshiftGrid.size + + # Loop over SEDs + # create a file per SED of all possible flux in band + for sed_name in sed_names: + seddata = np.genfromtxt(dir_seds + '/' + sed_name + fmt) + seddata[:, 1] *= seddata[:, 0] # SDC : multiply luminosity by wl ? + # SDC: OK if luminosity is in wl bins ! To be checked !!!! + ref = np.interp(lambdaRef, seddata[:, 0], seddata[:, 1]) + seddata[:, 1] /= ref # normalisation at lambdaRef + sed_interp = interp1d(seddata[:, 0], seddata[:, 1]) # interpolation + + # container of redshift/ flux : matrix n_z x n_b for each template + # each column correspond to fluxes in the different bands at a a fixed redshift + # redshift along row, fluxes along column + # model of flux as a function of redshift for each template + f_mod = np.zeros((redshiftGrid.size, len(bandNames))) + + # Loop over bands + # jf index on bands + for jf, band in enumerate(bandNames): + fname_in = dir_filters + '/' + band + '.res' + data = np.genfromtxt(fname_in) + xf, yf = data[:, 0], data[:, 1] + #yf /= xf # divide by lambda + # Only consider range where >1% max + ind = np.where(yf > 0.01*np.max(yf))[0] + lambdaMin, lambdaMax = xf[ind[0]], xf[ind[-1]] + norm = np.trapz(yf/xf, x=xf) # SDC: probably Cb + + # iz index on redshift + for iz in range(redshiftGrid.size): + opz = (redshiftGrid[iz] + 1) + xf_z = np.linspace(lambdaMin / opz, lambdaMax / opz, num=5000) + yf_z = interp1d(xf / opz, yf)(xf_z) + ysed = sed_interp(xf_z) + f_mod[iz, jf] = np.trapz(ysed * yf_z, x=xf_z) / norm + f_mod[iz, jf] *= opz**2. / DL(redshiftGrid[iz])**2. / (4*np.pi) + # for each SED, save the flux at each redshift (along row) for each + np.savetxt(dir_seds + '/' + sed_name + '_fluxredshiftmod.txt', f_mod) + + +#----------------------------------------------------------------------------------------- +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start processSEDs.py" + logger.info(msg) + logger.info("--- Process SEDs ---") + + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + processSEDs(sys.argv[1]) diff --git a/rail/estimation/algos/include_delightPZ/simulateWithSEDs.py b/rail/estimation/algos/include_delightPZ/simulateWithSEDs.py new file mode 100644 index 00000000..09d1b56a --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/simulateWithSEDs.py @@ -0,0 +1,146 @@ +####################################################################################################### +# +# script : simulateWithSED.py +# +# simulate mock data with those filters and SEDs +# produce files `galaxies-redshiftpdfs.txt` and `galaxies-redshiftpdfs2.txt` for training and target +# +######################################################################################################### + + +import sys +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +from delight.io import * +from delight.utils import * + + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + + +def simulateWithSEDs(configfilename): + """ + + :param configfilename: + :return: + """ + + + + + logger.info("--- Simulate with SED ---") + + params = parseParamFile(configfilename, verbose=False, catFilesNeeded=False) + dir_seds = params['templates_directory'] + sed_names = params['templates_names'] + + # redshift grid + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + + numZ = redshiftGrid.size + numT = len(sed_names) + numB = len(params['bandNames']) + numObjects = params['numObjects'] + noiseLevel = params['noiseLevel'] + + # f_mod : 2D-container of interpolation functions of flux over redshift: + # row sed, column bands + # one row per sed, one column per band + f_mod = np.zeros((numT, numB), dtype=object) + + # loop on SED + # read the fluxes file at different redshift in training data file + # in file sed_name + '_fluxredshiftmod.txt' + # to produce f_mod the interpolation function redshift --> flux for each band and sed template + for it, sed_name in enumerate(sed_names): + # data : redshifted fluxes (row vary with z, columns: filters) + data = np.loadtxt(dir_seds + '/' + sed_name + '_fluxredshiftmod.txt') + # build the interpolation of flux wrt redshift for each band + for jf in range(numB): + f_mod[it, jf] = interp1d(redshiftGrid, data[:, jf], kind='linear') + + # Generate training data + #------------------------- + # pick a set of redshift at random to be representative of training galaxies + redshifts = np.random.uniform(low=redshiftGrid[0],high=redshiftGrid[-1],size=numObjects) + #pick some SED type at random + types = np.random.randint(0, high=numT, size=numObjects) + + ell = 1e6 # I don't know why we have this value multiplicative constant + # it is to show that delightLearn can find this multiplicative number when calling + # utils:scalefree_flux_likelihood(returnedChi2=True) + #ell = 0.45e-4 # SDC may 14 2021 calibrate approximately to AB magnitude + + # what is fluxes and fluxes variance + fluxes, fluxesVar = np.zeros((numObjects, numB)), np.zeros((numObjects, numB)) + + # loop on objects to simulate for the training and save in output training file + for k in range(numObjects): + #loop on number of bands + for i in range(numB): + trueFlux = ell * f_mod[types[k], i](redshifts[k]) # noiseless flux at the random redshift + noise = trueFlux * noiseLevel + fluxes[k, i] = trueFlux + noise * np.random.randn() # noisy flux + fluxesVar[k, i] = noise**2. + + # container for training galaxies output + # at some redshift, provides the flux and its variance inside each band + data = np.zeros((numObjects, 1 + len(params['training_bandOrder']))) + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn,refBandColumn = readColumnPositions(params, prefix="training_") + + for ib, pf, pfv in zip(bandIndices, bandColumns, bandVarColumns): + data[:, pf] = fluxes[:, ib] + data[:, pfv] = fluxesVar[:, ib] + data[:, redshiftColumn] = redshifts + data[:, -1] = types + np.savetxt(params['trainingFile'], data) + + # Generate Target data : procedure similar to the training + #----------------------------------------------------------- + # pick set of redshift at random + redshifts = np.random.uniform(low=redshiftGrid[0],high=redshiftGrid[-1],size=numObjects) + types = np.random.randint(0, high=numT, size=numObjects) + + fluxes, fluxesVar = np.zeros((numObjects, numB)), np.zeros((numObjects, numB)) + + # loop on objects in target files + for k in range(numObjects): + # loop on bands + for i in range(numB): + # compute the flux in that band at the redshift + trueFlux = f_mod[types[k], i](redshifts[k]) + noise = trueFlux * noiseLevel + fluxes[k, i] = trueFlux + noise * np.random.randn() + fluxesVar[k, i] = noise**2. + + data = np.zeros((numObjects, 1 + len(params['target_bandOrder']))) + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn,refBandColumn = readColumnPositions(params, prefix="target_") + + for ib, pf, pfv in zip(bandIndices, bandColumns, bandVarColumns): + data[:, pf] = fluxes[:, ib] + data[:, pfv] = fluxesVar[:, ib] + data[:, redshiftColumn] = redshifts + data[:, -1] = types + np.savetxt(params['targetFile'], data) + + +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start simulateWithSEDs.py" + logger.info(msg) + logger.info("--- simulate with SED ---") + + + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + simulateWithSEDs(sys.argv[1]) diff --git a/rail/estimation/algos/include_delightPZ/templateFitting.py b/rail/estimation/algos/include_delightPZ/templateFitting.py new file mode 100644 index 00000000..a8acba8e --- /dev/null +++ b/rail/estimation/algos/include_delightPZ/templateFitting.py @@ -0,0 +1,205 @@ +######################################################################################## +# +# script : templateFitting.py +# +# Does the template fitting not calling gaussian processes +# +# output files : redshiftpdfFileTemp and metricsFileTemp +# +###################################################################################### +import sys +from mpi4py import MPI +import numpy as np +from scipy.interpolate import interp1d + +from delight.io import * +from delight.utils import * +from delight.photoz_gp import PhotozGP +from delight.photoz_kernels import Photoz_mean_function, Photoz_kernel + +from rail.estimation.algos.include_delightPZ.libPriorPZ import * + + + +import coloredlogs +import logging + + +logger = logging.getLogger(__name__) +coloredlogs.install(level='DEBUG', logger=logger,fmt='%(asctime)s,%(msecs)03d %(programname)s, %(name)s[%(process)d] %(levelname)s %(message)s') + +FLAG_NEW_PRIOR = True + +def templateFitting(configfilename): + """ + + :param configfilename: + :return: + """ + + comm = MPI.COMM_WORLD + threadNum = comm.Get_rank() + numThreads = comm.Get_size() + + if threadNum == 0: + logger.info("--- TEMPLATE FITTING ---") + + if FLAG_NEW_PRIOR: + logger.info("==> New Prior calculation from Benitez") + + # Parse parameters file + + paramFileName = configfilename + params = parseParamFile(paramFileName, verbose=False) + + if threadNum == 0: + msg = 'Thread number / number of threads: ' + str(threadNum+1) + " , " + str(numThreads) + logger.info(msg) + msg = 'Input parameter file:' + paramFileName + logger.info(msg) + + + + DL = approx_DL() + redshiftDistGrid, redshiftGrid, redshiftGridGP = createGrids(params) + numZ = redshiftGrid.size + + # Locate which columns of the catalog correspond to which bands. + + bandIndices, bandNames, bandColumns, bandVarColumns, redshiftColumn,refBandColumn = readColumnPositions(params, prefix="target_") + + dir_seds = params['templates_directory'] + dir_filters = params['bands_directory'] + lambdaRef = params['lambdaRef'] + sed_names = params['templates_names'] + + # f_mod : flux model in each band as a function of the sed and the band name + # axis 0 : redshifts + # axis 1 : sed names + # axis 2 : band names + + f_mod = np.zeros((redshiftGrid.size, len(sed_names),len(params['bandNames']))) + + # loop on SED to load the flux-redshift file from the training + # ture data or simulated by simulateWithSEDs.py + + for t, sed_name in enumerate(sed_names): + f_mod[:, t, :] = np.loadtxt(dir_seds + '/' + sed_name + '_fluxredshiftmod.txt') + + numObjectsTarget = np.sum(1 for line in open(params['target_catFile'])) + + firstLine = int(threadNum * numObjectsTarget / float(numThreads)) + lastLine = int(min(numObjectsTarget,(threadNum + 1) * numObjectsTarget / float(numThreads))) + numLines = lastLine - firstLine + + if threadNum == 0: + msg='Number of Target Objects ' + str(numObjectsTarget) + logger.info(msg) + + comm.Barrier() + + msg= 'Thread ' + str(threadNum) + ' , analyzes lines ' + str(firstLine) + ' , to ' + str(lastLine) + logger.info(msg) + + numMetrics = 7 + len(params['confidenceLevels']) + + # Create local files to store results + localPDFs = np.zeros((numLines, numZ)) + localMetrics = np.zeros((numLines, numMetrics)) + + # Now loop over each target galaxy (indexed bu loc index) to compute likelihood function + # with its flux in each bands + loc = - 1 + trainingDataIter = getDataFromFile(params, firstLine, lastLine,prefix="target_", getXY=False) + for z, normedRefFlux, bands, fluxes, fluxesVar,bCV, fCV, fvCV in trainingDataIter: + loc += 1 + # like_grid, _ = scalefree_flux_likelihood( + # fluxes, fluxesVar, + # f_mod[:, :, bands]) + # ell_hat_z = normedRefFlux * 4 * np.pi\ + # * params['fluxLuminosityNorm'] \ + # * (DL(redshiftGrid)**2. * (1+redshiftGrid))[:, None] + + # OLD way be keep it now + ell_hat_z = 1 + params['ellPriorSigma'] = 1e12 + + # Not working + #ell_hat_z=0.45e-4 + #params['ellPriorSigma'] = 1e12 + + # approximate flux likelihood, with scaling of both the mean and variance. + # This approximates the true likelihood with an iterative scheme. + # - data : fluxes, fluxesVar + # - model based on SED : f_mod + like_grid = approx_flux_likelihood(fluxes, fluxesVar, f_mod[:, :, bands],normalized=True, marginalizeEll=True,ell_hat=ell_hat_z, ell_var=(ell_hat_z*params['ellPriorSigma'])**2) + + if FLAG_NEW_PRIOR: + maglim=26 # M5 magnitude max + p_z = libPriorPZ(redshiftGrid,maglim=maglim) # return 2D template nz x nt, nt is 8 + + + else: + b_in = np.array(params['p_t'])[None, :] + beta2 = np.array(params['p_z_t'])**2.0 + + #compute prior on z + p_z = b_in * redshiftGrid[:, None] / beta2[None, :] *np.exp(-0.5 * redshiftGrid[:, None]**2 / beta2[None, :]) + + if loc < 0: + np.set_printoptions(threshold=20, edgeitems=10, linewidth=140,formatter=dict(float=lambda x: "%.3e" % x)) # float arrays %.3g + print(p_z) + + # Compute likelihood x prior + like_grid *= p_z + + localPDFs[loc, :] += like_grid.sum(axis=1) + + if localPDFs[loc, :].sum() > 0: + localMetrics[loc, :] = computeMetrics(z, redshiftGrid,localPDFs[loc, :],params['confidenceLevels']) + + comm.Barrier() + if threadNum == 0: + globalPDFs = np.zeros((numObjectsTarget, numZ)) + globalMetrics = np.zeros((numObjectsTarget, numMetrics)) + else: + globalPDFs = None + globalMetrics = None + + firstLines = [int(k*numObjectsTarget/numThreads) for k in range(numThreads)] + lastLines = [int(min(numObjectsTarget, (k+1)*numObjectsTarget/numThreads)) for k in range(numThreads)] + numLines = [lastLines[k] - firstLines[k] for k in range(numThreads)] + + sendcounts = tuple([numLines[k] * numZ for k in range(numThreads)]) + displacements = tuple([firstLines[k] * numZ for k in range(numThreads)]) + comm.Gatherv(localPDFs,[globalPDFs, sendcounts, displacements, MPI.DOUBLE]) + + sendcounts = tuple([numLines[k] * numMetrics for k in range(numThreads)]) + displacements = tuple([firstLines[k] * numMetrics for k in range(numThreads)]) + comm.Gatherv(localMetrics,[globalMetrics, sendcounts, displacements, MPI.DOUBLE]) + + comm.Barrier() + + if threadNum == 0: + fmt = '%.2e' + np.savetxt(params['redshiftpdfFileTemp'], globalPDFs, fmt=fmt) + if redshiftColumn >= 0: + np.savetxt(params['metricsFileTemp'], globalMetrics, fmt=fmt) + + + + +if __name__ == "__main__": + # execute only if run as a script + + + msg="Start templateFitting.py" + logger.info(msg) + logger.info("--- Template Fitting ---") + + + + if len(sys.argv) < 2: + raise Exception('Please provide a parameter file') + + templateFitting(sys.argv[1]) diff --git a/tests/data/parametersTest.cfg b/tests/data/parametersTest.cfg new file mode 100644 index 00000000..2568e666 --- /dev/null +++ b/tests/data/parametersTest.cfg @@ -0,0 +1,74 @@ + +# DELIGHT parameter file +# Syntactic rules: +# - You can set parameters with : or = +# - Lines starting with # or ; will be ignored +# - Multiple values (band names, band orders, confidence levels) +# must beb separated by spaces +# - The input files should contain numbers separated with spaces. +# - underscores mean unused column + +[Bands] +names: lsst_u lsst_g lsst_r lsst_i lsst_z lsst_y +directory: /Users/sylvie/MacOSX/GitHub/LSST/Delight/data/FILTERS +numCoefs: 15 +bands_verbose: True +bands_debug: True +bands_makeplots: True + +[Templates] +directory: /Users/sylvie/MacOSX/GitHub/LSST/Delight/data/CWW_SEDs +names: El_B2004a Sbc_B2004a Scd_B2004a SB3_B2004a SB2_B2004a Im_B2004a ssp_25Myr_z008 ssp_5Myr_z008 +p_t: 0.27 0.26 0.25 0.069 0.021 0.11 0.0061 0.0079 +p_z_t:0.23 0.39 0.33 0.31 1.1 0.34 1.2 0.14 +lambdaRef: 4.5e3 + +[Simulation] +numObjects: 1000 +noiseLevel: 0.03 +trainingFile: ./tmp/data/galaxies-fluxredshifts.txt +targetFile: ./tmp/galaxies-fluxredshifts2.txt + +[Training] +catFile: ./tmp/data/galaxies-fluxredshifts.txt +bandOrder: lsst_u lsst_u_var lsst_g lsst_g_var _ _ lsst_i lsst_i_var lsst_z lsst_z_var lsst_y lsst_y_var redshift +referenceBand: lsst_i +extraFracFluxError: 1e-4 +paramFile: ./tmp/galaxies-gpparams.txt +crossValidate: False +CVfile: data_lsst/galaxies-gpCV.txt +crossValidationBandOrder: _ _ _ _ lsst_r lsst_r_var _ _ _ _ _ _ +numChunks: 1 + +[Target] +catFile: ./tmp/data/galaxies-fluxredshifts2.txt +bandOrder: lsst_u lsst_u_var lsst_g lsst_g_var lsst_r lsst_r_var _ _ lsst_z lsst_z_var lsst_y lsst_y_var redshift +referenceBand: lsst_r +extraFracFluxError: 1e-4 +redshiftpdfFile: ./tmp/data/galaxies-redshiftpdfs.txt +redshiftpdfFileTemp: ./tmp/data/galaxies-redshiftpdfs-cww.txt +metricsFile: ./tmp/data/galaxies-redshiftmetrics.txt +metricsFileTemp: ./tmp/data/galaxies-redshiftmetrics-cww.txt +useCompression: False +Ncompress: 10 +compressIndicesFile: ./tmp/galaxies-compressionIndices.txt +compressMargLikFile: ./tmp/galaxies-compressionMargLikes.txt +redshiftpdfFileComp: ./tmp/galaxies-redshiftpdfs-comp.txt + +[Other] +rootDir: ./ +zPriorSigma: 0.2 +ellPriorSigma: 0.5 +fluxLuminosityNorm: 1.0 +alpha_C: 1.0e3 +V_C: 0.1 +alpha_L: 1.0e2 +V_L: 0.1 +lines_pos: 6500 5002.26 3732.22 +lines_width: 20.0 20.0 20.0 +redshiftMin: 0.1 +redshiftMax: 1.101 +redshiftNumBinsGPpred: 100 +redshiftBinSize: 0.001 +redshiftDisBinSize: 0.2 +confidenceLevels: 0.1 0.50 0.68 0.95