ACCarnall / bagpipes

Bagpipes is a state of the art code for generating realistic model galaxy spectra and fitting these to spectroscopic and photometric observations. Users should install with pip, not by cloning the repository.
http://bagpipes.readthedocs.io
GNU General Public License v3.0
71 stars 37 forks source link

Something about the SNR setting #37

Closed markliuchina closed 1 year ago

markliuchina commented 1 year ago

Dear Adams: Hello! I was trying to imitate the work of your paper (Carnall et al. 2020) when learning your code. I followed the instructions in your paper and prepared all the input files well. (I've checked the filter files and the photometry catalog and find no problems, as well as the double-power-law fitting model setting) At first, I set the function (load the photometry data) as:

def load_goodss(ID):
    cat = np.loadtxt('./photometry_catalog/staa1535_SED_input_catalog_goodss.txt', 
                 usecols = (3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35,   
                            4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36))
    row = int(ID) - 1
    fluxes = cat[row, :17]
    fluxerrs = cat[row, 17:]
    photometry = np.c_[fluxes, fluxerrs]

    for i in range(len(photometry)):
        if (photometry[i, 0] == 0.) or (photometry[i, 1] <= 0):
            photometry[i,:] = [0., 9.9*10**99.]
    for i in range(len(photometry)):
         #This SNR criteria is from the published paper 
        if i <= 12:
            max_snr = 20.
        elif i <= 14:
            max_snr = 10.
        else:
            max_snr = 5.

        if photometry[i, 0]/photometry[i, 1] > max_snr:
            photometry[i, 1] = photometry[i, 0]/max_snr
    return photometry

I changed the SNR criteria and now it is different from the SNR setting in your instruction example jupyter-notebook. And the fitting process crashed with the terminal showing me such infomation:

ValueError                                Traceback (most recent call last)
/tmp/ipykernel_140663/828412448.py in <module>
      1 fit = pipes.fit(galaxy, fit_instructions)
----> 2 fit.fit(verbose = False)
      3 fig = fit.plot_spectrum_posterior(save = False, show = True)

~/anaconda3/lib/python3.9/site-packages/bagpipes/fitting/fit.py in fit(self, verbose, n_live, use_MPI)
    134 
    135             # Load MultiNest outputs and save basic quantities to file.
--> 136             samples2d = np.loadtxt(self.fname + "post_equal_weights.dat")
    137             lnz_line = open(self.fname + "stats.dat").readline().split()
    138 

~/anaconda3/lib/python3.9/site-packages/numpy/lib/npyio.py in loadtxt(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin, encoding, max_rows, like)
   1146         # converting the data
   1147         X = None
-> 1148         for x in read_data(_loadtxt_chunksize):
   1149             if X is None:
   1150                 X = np.array(x, dtype)

~/anaconda3/lib/python3.9/site-packages/numpy/lib/npyio.py in read_data(chunk_size)
    997 
    998             # Convert each value according to its column and store
--> 999             items = [conv(val) for (conv, val) in zip(converters, vals)]
   1000 
   1001             # Then pack it according to the dtype's nesting

~/anaconda3/lib/python3.9/site-packages/numpy/lib/npyio.py in <listcomp>(.0)
    997 
    998             # Convert each value according to its column and store
--> 999             items = [conv(val) for (conv, val) in zip(converters, vals)]
   1000 
   1001             # Then pack it according to the dtype's nesting

~/anaconda3/lib/python3.9/site-packages/numpy/lib/npyio.py in floatconv(x)
    734         if '0x' in x:
    735             return float.fromhex(x)
--> 736         return float(x)
    737 
    738     typ = dtype.type

ValueError: could not convert string to float: '0.695279959934739884-309'

And also a red warning window from jupyter-notebook:

Exception ignored on calling ctypes callback function: <function run.<locals>.loglike at 0x7f06c171c8b0>
Traceback (most recent call last):
  File "/home/ubuntu/anaconda3/lib/python3.9/site-packages/pymultinest/run.py", line 220, in loglike
    Prior(cube, ndim, nparams)
  File "/home/ubuntu/anaconda3/lib/python3.9/site-packages/bagpipes/fitting/prior.py", line 75, in transform
    cube[i] = prior_function(cube[i], self.limits[i],
  File "/home/ubuntu/anaconda3/lib/python3.9/site-packages/bagpipes/fitting/prior.py", line 88, in log_10
    value = 10**((np.log10(limits[1]/limits[0]))*value
ZeroDivisionError: float division by zero 

I checked the filter files, model setting and the photometry data again but it still failed. Finally, when I recover back the SNR setting, following your example, like this:

def load_goodss(ID):
    cat = np.loadtxt('./photometry_catalog/staa1535_SED_input_catalog_goodss.txt', 
                 usecols = (3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35,   
                            4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36))
    row = int(ID) - 1
    fluxes = cat[row, :17]
    fluxerrs = cat[row, 17:]
    photometry = np.c_[fluxes, fluxerrs]

    for i in range(len(photometry)):
        if (photometry[i, 0] == 0.) or (photometry[i, 1] <= 0):
            photometry[i,:] = [0., 9.9*10**99.]
    for i in range(len(photometry)):

         for i in range(len(photometry)):
            if i < 10:
                max_snr = 20.
            else:
                max_snr = 10.
        if photometry[i, 0]/photometry[i, 1] > max_snr:
                photometry[i, 1] = photometry[i, 0]/max_snr
    return photometry

When I start running again, it succeed! I know about the concerns and the reasons of the SNR criteria and indeed enforced a SNR criteria as described in your paper when selecting the massive galaxies sample. I was thinking, in bagpipes fitting, maybe a simple SNR limiting (max_snr<=20 for all, max_snr<=10 for IRAC Channels) is just fine and will not influence the fitting results? (Because the SNR criteria described in your paper was indeed performed and there is no need to limit it so strictly?) Maybe I have a misunderstanding of it. (That precise SNR criteria is just for sample selection and has little impact on the fitting) Could you please tell me how should I set the SNR limit when fitting? Anything from you will be appreciated! Thanks you again for providing this wonderful code! Best wishes!

Mingfeng Liu

Department of Theoretical Physics Nanjing Normal University Wenyuan Road 1, Qixia District Nanjing City, Jiangsu Province, PRC 2023.3.3